ceph/src/boost/libs/math/doc/roots/minima.qbk

   1 [section:brent_minima Locating Function Minima using Brent's algorithm]
   2
   3 [import ../../example/brent_minimise_example.cpp]
   4
   5 [h4 Synopsis]
   6
   7 ``
   8 #include <boost/math/tools/minima.hpp>
   9
  10 ``
  11
  12    template <class F, class T>
  13    std::pair<T, T> brent_find_minima(F f, T min, T max, int bits);
  14
  15    template <class F, class T>
  16    std::pair<T, T> brent_find_minima(F f, T min, T max, int bits, boost::uintmax_t& max_iter);
  17
  18 [h4 Description]
  19
  20 These two functions locate the minima of the continuous function ['f] using
  21 [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method]: specifically it
  22 uses quadratic interpolation to locate the minima, or if that fails, falls back to
  23 a [@http://en.wikipedia.org/wiki/Golden_section_search golden-section search].
  24
  25 [*Parameters]
  26
  27 [variablelist
  28 [[f] [The function to minimise: a function object (functor) that should be smooth over the
  29           range ['\[min, max\]], with no maxima occurring in that interval.]]
  30 [[min] [The lower endpoint of the range in which to search  for the minima.]]
  31 [[max] [The upper endpoint of the range in which to search  for the minima.]]
  32 [[bits] [The number of bits precision to which the minima should be found.[br]
  33              Note that in principle, the minima can not be located to greater
  34              accuracy than the square root of machine epsilon (for 64-bit double, sqrt(1e-16)[cong]1e-8),
  35              therefore the value of ['bits] will be ignored if it's greater than half the number of bits
  36              in the mantissa of T.]]
  37 [[max_iter] [The maximum number of iterations to use
  38              in the algorithm, if not provided the algorithm will just
  39              keep on going until the minima is found.]]
  40 ] [/variablelist]
  41
  42 [*Returns:]
  43
  44 A `pair` of type T containing the value of the abscissa at the minima and the value
  45 of ['f(x)] at the minima.
  46
  47 [tip Defining BOOST_MATH_INSTRUMENT will show some parameters, for example:
  48 ``
  49         Type T is double
  50         bits = 24, maximum 26
  51         tolerance = 1.19209289550781e-007
  52         seeking minimum in range min-4 to 1.33333333333333
  53         maximum iterations 18446744073709551615
  54         10 iterations.
  55 ``
  56 ]
  57
  58 [h4:example Brent Minimisation Example]
  59
  60 As a demonstration, we replicate this  [@http://en.wikipedia.org/wiki/Brent%27s_method#Example Wikipedia example]
  61 minimising the function ['y= (x+3)(x-1)[super 2]].
  62
  63 It is obvious from the equation and the plot that there is a
  64 minimum at exactly one and the value of the function at one is exactly zero.
  65
  66 [tip This observation shows that an analytical or
  67 [@http://en.wikipedia.org/wiki/Closed-form_expression Closed-form expression]
  68 solution always beats brute-force  hands-down for both speed and precision.]
  69
  70 [graph brent_test_function_1]
  71
  72 First an include is needed:
  73
  74 [brent_minimise_include_1]
  75
  76 This function is encoded in C++ as function object (functor) using `double` precision thus:
  77
  78 [brent_minimise_double_functor]
  79
  80 The Brent function is conveniently accessed through a `using` statement (noting sub-namespace `::tools`).
  81
  82 The search minimum and maximum are chosen as -4 to 4/3 (as in the Wikipedia example).
  83
  84 [tip S A Stage (reference 6) reports that the Brent algorithm is ['slow to start, but fast to converge],
  85 so choosing a tight min-max range is good.]
  86
  87 For simplicity, we set the precision parameter `bits` to `std::numeric_limits<double>::digits`,
  88 which is effectively the maximum possible i.e. `std::numeric_limits<double>::digits`/2.
  89 Nor do we provide a maximum iterations parameter `max_iter`,
  90 (perhaps unwidely), so the function will iterate until it finds a minimum.
  91
  92 [brent_minimise_double_1]
  93
  94 The resulting [@http://en.cppreference.com/w/cpp/utility/pair std::pair]
  95 contains the minimum close to one and the minimum value close to zero.
  96
  97   x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018
  98
  99 The differences from the expected ['one] and ['zero] are less than the
 100 uncertainty (for `double`) 1.5e-008 calculated from
 101 `sqrt(std::numeric_limits<double>::digits) == 53`.
 102
 103 We can use it like this to check that the two values are close-enough to those expected,
 104
 105     using boost::math::fpc::is_close_to;
 106     using boost::math::fpc::is_small;
 107
 108    double uncertainty = sqrt(std::numeric_limits<double>::digits);
 109    is_close_to(1., r.first, uncertainty);
 110    is_small(r.second, uncertainty);
 111
 112    x == 1 (compared to uncertainty 0.00034527) is true
 113    f(x) == 0 (compared to uncertainty 0.00034527) is true
 114
 115 It is possible to make this comparison more generally with a templated function,
 116 returning `true` when this criterion is met, for example:
 117
 118 [brent_minimise_close]
 119
 120 In practical applications, we might want to know how many iterations,
 121 and maybe to limit iterations and
 122 perhaps to trade some loss of precision for speed, for example:
 123
 124 [brent_minimise_double_2]
 125
 126 limits to a maximum of 20 iterations
 127 (a reasonable estimate for this application, even for higher precision shown later).
 128
 129 The parameter `it` is updated to return the actual number of iterations
 130 (so it may be useful to also keep a record of the limit in `maxit`).
 131
 132 It is neat to avoid showing insignificant digits by computing the number of decimal digits to display.
 133
 134 [brent_minimise_double_3]
 135
 136   Showing 53 bits precision with 9 decimal digits from tolerance 1.49011611938477e-008
 137   x at minimum = 1, f(1) = 5.04852568e-018
 138
 139 We can also half the number of precision bits from 52 to 26.
 140
 141 [brent_minimise_double_4]
 142
 143 showing no change in the result and no change in the number of iterations, as expected.
 144
 145 It is only if we reduce the precision to a quarter, specifying only 13 precision bits
 146
 147 [brent_minimise_double_5]
 148
 149 that we reduce the number of iterations from 10 to 7 and the result significantly differing from ['one] and ['zero].
 150
 151   Showing 13 bits precision with 9 decimal digits from tolerance 0.015625
 152   x at minimum = 0.9999776, f(0.9999776) = 2.0069572e-009 after  7 iterations.
 153
 154 [h5:template Templating on floating-point type]
 155
 156 If we want to switch the floating-point type, then the functor must be revised.
 157 Since the functor is stateless, the easiest option is to simply make
 158 `operator()` a template member function:
 159
 160 [brent_minimise_T_functor]
 161
 162 The `brent_find_minima` function can now be used in template form.
 163
 164 [brent_minimise_template_1]
 165
 166 The form shown uses the floating-point type `long double` by deduction,
 167 but it is also possible to be more explicit, for example:
 168
 169   std::pair<long double, long double> r = brent_find_minima<func, long double>
 170   (func(), bracket_min, bracket_max, bits, it);
 171
 172 In order to show the use of multiprecision below, it may be convenient to write a templated function to use this.
 173
 174 [brent_minimise_T_show]
 175
 176 We can use this with all built-in floating-point types, for example
 177
 178 [brent_minimise_template_fd]
 179
 180 and, on platforms that provide it, a
 181 [@http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format 128-bit quad] type.
 182 (See [@boost:libs/multiprecision/doc/html/boost_multiprecision/tut/floats/float128.html float128]).
 183
 184 For this optional include, the build should define the macro BOOST_HAVE_QUADMATH:
 185
 186 [brent_minimise_mp_include_1]
 187
 188 or
 189
 190 [brent_minimise_template_quad]
 191
 192 [h5:multiprecision  Multiprecision]
 193
 194 If a higher precision than `double` (or `long double` if that is more precise) is required,
 195 then this is easily achieved using __multiprecision with some includes from
 196
 197 [brent_minimise_mp_include_0]
 198
 199 and some `typdef`s.
 200
 201 [brent_minimise_mp_typedefs]
 202 Using thus
 203
 204 [brent_minimise_mp_1]
 205
 206 and with our show function
 207
 208 [brent_minimise_mp_2]
 209
 210 [brent_minimise_mp_output_1]
 211
 212 [brent_minimise_mp_output_2]
 213
 214 [tip One can usually rely on template argument deduction
 215 to avoid specifying the verbose multiprecision types,
 216 but great care in needed with the ['type of the values] provided
 217 to avoid confusing the compiler.
 218 ]
 219
 220 [tip Using `std::cout.precision(std::numeric_limits<T>::digits10);`
 221 or `std::cout.precision(std::numeric_limits<T>::max_digits10);`
 222 during debugging may be wise because it gives some warning if construction of multiprecision values
 223 involves unintended conversion from `double` by showing trailing zero  or random  digits after
 224 [@http://en.cppreference.com/w/cpp/types/numeric_limits/max_digits10 max_digits10],
 225 that is 17 for `double`, digit 18... may be just noise.]
 226
 227 The complete example code is at [@../../example/brent_minimise_example.cpp brent_minimise_example.cpp].
 228
 229 [h4 Implementation]
 230
 231 This is a reasonably faithful implementation of Brent's algorithm.
 232
 233 [h4 References]
 234
 235 # Brent, R.P. 1973, Algorithms for Minimization without Derivatives,
 236 (Englewood Cliffs, NJ: Prentice-Hall), Chapter 5.
 237
 238 # Numerical Recipes in C, The Art of Scientific Computing,
 239 Second Edition, William H. Press, Saul A. Teukolsky,
 240 William T. Vetterling, and Brian P. Flannery.
 241 Cambridge University Press. 1988, 1992.
 242
 243 # An algorithm with guaranteed convergence for finding a zero
 244 of a function, R. P. Brent, The Computer Journal, Vol 44, 1971.
 245
 246 # [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method in Wikipedia.]
 247
 248 # Z. Zhang, An Improvement to the Brent's Method, IJEA, vol. 2, pp. 2 to 26, May 31, 2011.
 249 [@http://www.cscjournals.org/manuscript/Journals/IJEA/volume2/Issue1/IJEA-7.pdf ]
 250
 251 # Steven A. Stage, Comments on An Improvement to the Brent's Method
 252 (and comparison of various algorithms)
 253 [@http://www.cscjournals.org/manuscript/Journals/IJEA/volume4/Issue1/IJEA-33.pdf]
 254 Stage concludes that Brent's algorithm is slow to start, but fast to finish convergence, and has good accuracy.
 255
 256 [endsect] [/section:rebt_minima Locating Function Minima]
 257
 258 [/
 259   Copyright 2006, 2015 John Maddock and Paul A. Bristow.
 260   Distributed under the Boost Software License, Version 1.0.
 261   (See accompanying file LICENSE_1_0.txt or copy at
 262   http://www.boost.org/LICENSE_1_0.txt).
 263 ]
 264