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 20These 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 22uses quadratic interpolation to locate the minima, or if that fails, falls back to 23a [@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 (or C++ lambda) 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 44A `std::pair` of type T containing the value of the abscissa (x) at the minima and the value 45of ['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 60As a demonstration, we replicate this [@http://en.wikipedia.org/wiki/Brent%27s_method#Example Wikipedia example] 61minimising the function ['y= (x+3)(x-1)[super 2]]. 62 63It is obvious from the equation and the plot that there is a 64minimum at exactly unity (x = 1) and the value of the function at one is exactly zero (y = 0). 65 66[tip This observation shows that an analytical or 67[@http://en.wikipedia.org/wiki/Closed-form_expression Closed-form expression] 68solution always beats brute-force hands-down for both speed and precision.] 69 70[graph brent_test_function_1] 71 72First an include is needed: 73 74[brent_minimise_include_1] 75 76This function is encoded in C++ as function object (functor) using `double` precision thus: 77 78[brent_minimise_double_functor] 79 80The Brent function is conveniently accessed through a `using` statement (noting sub-namespace `::tools`). 81 82 using boost::math::tools::brent_find_minima; 83 84The search minimum and maximum are chosen as -4 to 4/3 (as in the Wikipedia example). 85 86[tip S A Stage (reference 6) reports that the Brent algorithm is ['slow to start, but fast to converge], 87so choosing a tight min-max range is good.] 88 89For simplicity, we set the precision parameter `bits` to `std::numeric_limits<double>::digits`, 90which is effectively the maximum possible `std::numeric_limits<double>::digits/2`. 91 92Nor do we provide a value for maximum iterations parameter `max_iter`, 93(probably unwisely), so the function will iterate until it finds a minimum. 94 95[brent_minimise_double_1] 96 97The resulting [@http://en.cppreference.com/w/cpp/utility/pair std::pair] 98contains the minimum close to one, and the minimum value close to zero. 99 100 x at minimum = 1.00000000112345, 101 f(1.00000000112345) = 5.04852568272458e-018 102 103The differences from the expected ['one] and ['zero] are less than the 104uncertainty, for `double` 1.5e-008, calculated from 105`sqrt(std::numeric_limits<double>::epsilon())`. 106 107We can use this uncertainty to check that the two values are close-enough to those expected, 108 109[brent_minimise_double_1a] 110 111that outputs 112 113 x == 1 (compared to uncertainty 1.49011611938477e-08) is true 114 f(x) == 0 (compared to uncertainty 1.49011611938477e-08) is true 115 116[note The function `close_at_tolerance` is ineffective for testing if a value is small or zero 117(because it may divide by small or zero and cause overflow). 118Function `is_small` does this job.] 119 120It is possible to make this comparison more generally with a templated function, 121`is_close`, first checking `is_small` before checking `close_at_tolerance`, 122returning `true` when small or close, for example: 123 124[brent_minimise_close] 125 126In practical applications, we might want to know how many iterations, 127and maybe to limit iterations (in case the function proves ill-behaved), 128and/or perhaps to trade some loss of precision for speed, for example: 129 130[brent_minimise_double_2] 131 132limits to a maximum of 20 iterations 133(a reasonable estimate for this example function, even for much higher precision shown later). 134 135The parameter `it` is updated to return the actual number of iterations 136(so it may be useful to also keep a record of the limit set in `const maxit`). 137 138It is neat to avoid showing insignificant digits by computing the number of decimal digits to display. 139 140[brent_minimise_double_3] 141 142 Showing 53 bits precision with 9 decimal digits from tolerance 1.49011611938477e-008 143 x at minimum = 1, f(1) = 5.04852568e-018 144 145We can also half the number of precision bits from 52 to 26: 146 147[brent_minimise_double_4] 148 149showing no change in the result and no change in the number of iterations, as expected. 150 151It is only if we reduce the precision to a quarter, specifying only 13 precision bits 152 153[brent_minimise_double_5] 154 155that we reduce the number of iterations from 10 to 7 that the result slightly differs from ['one] and ['zero]. 156 157 Showing 13 bits precision with 9 decimal digits from tolerance 0.015625 158 x at minimum = 0.9999776, f(0.9999776) = 2.0069572e-009 after 7 iterations. 159 160[h5:template Templating on floating-point type] 161 162If we want to switch the floating-point type, then the functor must be revised. 163Since the functor is stateless, the easiest option is to simply make 164`operator()` a template member function: 165 166[brent_minimise_T_functor] 167 168The `brent_find_minima` function can now be used in template form, for example using `long double`: 169 170[brent_minimise_template_1] 171 172The form shown uses the floating-point type `long double` by deduction, 173but it is also possible to be more explicit, for example: 174 175 std::pair<long double, long double> r = brent_find_minima<func, long double> 176 (func(), bracket_min, bracket_max, bits, it); 177 178In order to show the use of multiprecision below (or other user-defined types), 179it may be convenient to write a templated function to use this: 180 181[brent_minimise_T_show] 182 183[note the prudent addition of `try'n'catch` blocks within the function. 184This ensures that any malfunction will provide a Boost.Math error message 185rather than uncommunicatively calling `std::abort`.] 186 187We can use this with all built-in floating-point types, for example 188 189[brent_minimise_template_fd] 190 191[h5:quad_precision Quad 128-bit precision] 192 193On platforms that provide it, a 194[@http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format 128-bit quad] type can be used. 195(See [@boost:libs/multiprecision/doc/html/boost_multiprecision/tut/floats/float128.html float128]). 196 197If this type is available, the build should define the macro BOOST_HAVE_QUADMATH: 198 199[brent_minimise_mp_include_1] 200 201and can be (conditionally) used this: 202 203[brent_minimise_template_quad] 204 205[h5:multiprecision Multiprecision] 206 207If a higher precision than `double` (or `long double` if that is more precise) is required, 208then this is easily achieved using __multiprecision with some includes, for example: 209 210[brent_minimise_mp_include_0] 211 212and some `typdef`s. 213 214[brent_minimise_mp_typedefs] 215 216Used thus 217 218[brent_minimise_mp_1] 219 220and with our `show_minima` function 221 222[brent_minimise_mp_2] 223 224[brent_minimise_mp_output_1] 225 226[brent_minimise_mp_output_2] 227 228[tip One can usually rely on template argument deduction 229to avoid specifying the verbose multiprecision types, 230but great care in needed with the ['type of the values] provided 231to avoid confusing the compiler. 232] 233 234[tip Using `std::cout.precision(std::numeric_limits<T>::digits10);` 235or `std::cout.precision(std::numeric_limits<T>::max_digits10);` 236during debugging may be wise because it gives some warning if construction of multiprecision values 237involves unintended conversion from `double` by showing trailing zero or random digits after 238[@http://en.cppreference.com/w/cpp/types/numeric_limits/max_digits10 max_digits10], 239that is 17 for `double`, digit 18... may be just noise.] 240 241The complete example code is at [@../../example/brent_minimise_example.cpp brent_minimise_example.cpp]. 242 243[h4 Implementation] 244 245This is a reasonably faithful implementation of Brent's algorithm. 246 247[h4 References] 248 249# Brent, R.P. 1973, Algorithms for Minimization without Derivatives, 250(Englewood Cliffs, NJ: Prentice-Hall), Chapter 5. 251 252# Numerical Recipes in C, The Art of Scientific Computing, 253Second Edition, William H. Press, Saul A. Teukolsky, 254William T. Vetterling, and Brian P. Flannery. 255Cambridge University Press. 1988, 1992. 256 257# An algorithm with guaranteed convergence for finding a zero 258of a function, R. P. Brent, The Computer Journal, Vol 44, 1971. 259 260# [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method in Wikipedia.] 261 262# Z. Zhang, An Improvement to the Brent's Method, IJEA, vol. 2, pp. 2 to 26, May 31, 2011. 263[@http://www.cscjournals.org/manuscript/Journals/IJEA/volume2/Issue1/IJEA-7.pdf ] 264 265# Steven A. Stage, Comments on An Improvement to the Brent's Method 266(and comparison of various algorithms) 267[@http://www.cscjournals.org/manuscript/Journals/IJEA/volume4/Issue1/IJEA-33.pdf] 268Stage concludes that Brent's algorithm is slow to start, but fast to finish convergence, and has good accuracy. 269 270[endsect] [/section:rebt_minima Locating Function Minima] 271 272[/ 273 Copyright 2006, 2015, 2018 John Maddock and Paul A. Bristow. 274 Distributed under the Boost Software License, Version 1.0. 275 (See accompanying file LICENSE_1_0.txt or copy at 276 http://www.boost.org/LICENSE_1_0.txt). 277] 278 279