[section:brent_minima Locating Function Minima using Brent's algorithm] [import ../../example/brent_minimise_example.cpp] [h4 Synopsis] `` #include `` template std::pair brent_find_minima(F f, T min, T max, int bits); template std::pair brent_find_minima(F f, T min, T max, int bits, boost::uintmax_t& max_iter); [h4 Description] These two functions locate the minima of the continuous function ['f] using [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method]: specifically it uses quadratic interpolation to locate the minima, or if that fails, falls back to a [@http://en.wikipedia.org/wiki/Golden_section_search golden-section search]. [*Parameters] [variablelist [[f] [The function to minimise: a function object (or C++ lambda) that should be smooth over the range ['\[min, max\]], with no maxima occurring in that interval.]] [[min] [The lower endpoint of the range in which to search for the minima.]] [[max] [The upper endpoint of the range in which to search for the minima.]] [[bits] [The number of bits precision to which the minima should be found.[br] Note that in principle, the minima can not be located to greater accuracy than the square root of machine epsilon (for 64-bit double, sqrt(1e-16)[cong]1e-8), therefore the value of ['bits] will be ignored if it's greater than half the number of bits in the mantissa of T.]] [[max_iter] [The maximum number of iterations to use in the algorithm, if not provided the algorithm will just keep on going until the minima is found.]] ] [/variablelist] [*Returns:] A `std::pair` of type T containing the value of the abscissa (x) at the minima and the value of ['f(x)] at the minima. [tip Defining BOOST_MATH_INSTRUMENT will show some parameters, for example: `` Type T is double bits = 24, maximum 26 tolerance = 1.19209289550781e-007 seeking minimum in range min-4 to 1.33333333333333 maximum iterations 18446744073709551615 10 iterations. `` ] [h4:example Brent Minimisation Example] As a demonstration, we replicate this [@http://en.wikipedia.org/wiki/Brent%27s_method#Example Wikipedia example] minimising the function ['y= (x+3)(x-1)[super 2]]. It is obvious from the equation and the plot that there is a minimum at exactly unity (x = 1) and the value of the function at one is exactly zero (y = 0). [tip This observation shows that an analytical or [@http://en.wikipedia.org/wiki/Closed-form_expression Closed-form expression] solution always beats brute-force hands-down for both speed and precision.] [graph brent_test_function_1] First an include is needed: [brent_minimise_include_1] This function is encoded in C++ as function object (functor) using `double` precision thus: [brent_minimise_double_functor] The Brent function is conveniently accessed through a `using` statement (noting sub-namespace `::tools`). using boost::math::tools::brent_find_minima; The search minimum and maximum are chosen as -4 to 4/3 (as in the Wikipedia example). [tip S A Stage (reference 6) reports that the Brent algorithm is ['slow to start, but fast to converge], so choosing a tight min-max range is good.] For simplicity, we set the precision parameter `bits` to `std::numeric_limits::digits`, which is effectively the maximum possible `std::numeric_limits::digits/2`. Nor do we provide a value for maximum iterations parameter `max_iter`, (probably unwisely), so the function will iterate until it finds a minimum. [brent_minimise_double_1] The resulting [@http://en.cppreference.com/w/cpp/utility/pair std::pair] contains the minimum close to one, and the minimum value close to zero. x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018 The differences from the expected ['one] and ['zero] are less than the uncertainty, for `double` 1.5e-008, calculated from `sqrt(std::numeric_limits::epsilon())`. We can use this uncertainty to check that the two values are close-enough to those expected, [brent_minimise_double_1a] that outputs x == 1 (compared to uncertainty 1.49011611938477e-08) is true f(x) == 0 (compared to uncertainty 1.49011611938477e-08) is true [note The function `close_at_tolerance` is ineffective for testing if a value is small or zero (because it may divide by small or zero and cause overflow). Function `is_small` does this job.] It is possible to make this comparison more generally with a templated function, `is_close`, first checking `is_small` before checking `close_at_tolerance`, returning `true` when small or close, for example: [brent_minimise_close] In practical applications, we might want to know how many iterations, and maybe to limit iterations (in case the function proves ill-behaved), and/or perhaps to trade some loss of precision for speed, for example: [brent_minimise_double_2] limits to a maximum of 20 iterations (a reasonable estimate for this example function, even for much higher precision shown later). The parameter `it` is updated to return the actual number of iterations (so it may be useful to also keep a record of the limit set in `const maxit`). It is neat to avoid showing insignificant digits by computing the number of decimal digits to display. [brent_minimise_double_3] Showing 53 bits precision with 9 decimal digits from tolerance 1.49011611938477e-008 x at minimum = 1, f(1) = 5.04852568e-018 We can also half the number of precision bits from 52 to 26: [brent_minimise_double_4] showing no change in the result and no change in the number of iterations, as expected. It is only if we reduce the precision to a quarter, specifying only 13 precision bits [brent_minimise_double_5] that we reduce the number of iterations from 10 to 7 that the result slightly differs from ['one] and ['zero]. Showing 13 bits precision with 9 decimal digits from tolerance 0.015625 x at minimum = 0.9999776, f(0.9999776) = 2.0069572e-009 after 7 iterations. [h5:template Templating on floating-point type] If we want to switch the floating-point type, then the functor must be revised. Since the functor is stateless, the easiest option is to simply make `operator()` a template member function: [brent_minimise_T_functor] The `brent_find_minima` function can now be used in template form, for example using `long double`: [brent_minimise_template_1] The form shown uses the floating-point type `long double` by deduction, but it is also possible to be more explicit, for example: std::pair r = brent_find_minima (func(), bracket_min, bracket_max, bits, it); In order to show the use of multiprecision below (or other user-defined types), it may be convenient to write a templated function to use this: [brent_minimise_T_show] [note the prudent addition of `try'n'catch` blocks within the function. This ensures that any malfunction will provide a Boost.Math error message rather than uncommunicatively calling `std::abort`.] We can use this with all built-in floating-point types, for example [brent_minimise_template_fd] [h5:quad_precision Quad 128-bit precision] On platforms that provide it, a [@http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format 128-bit quad] type can be used. (See [@boost:libs/multiprecision/doc/html/boost_multiprecision/tut/floats/float128.html float128]). If this type is available, the build should define the macro BOOST_HAVE_QUADMATH: [brent_minimise_mp_include_1] and can be (conditionally) used this: [brent_minimise_template_quad] [h5:multiprecision Multiprecision] If a higher precision than `double` (or `long double` if that is more precise) is required, then this is easily achieved using __multiprecision with some includes, for example: [brent_minimise_mp_include_0] and some `typdef`s. [brent_minimise_mp_typedefs] Used thus [brent_minimise_mp_1] and with our `show_minima` function [brent_minimise_mp_2] [brent_minimise_mp_output_1] [brent_minimise_mp_output_2] [tip One can usually rely on template argument deduction to avoid specifying the verbose multiprecision types, but great care in needed with the ['type of the values] provided to avoid confusing the compiler. ] [tip Using `std::cout.precision(std::numeric_limits::digits10);` or `std::cout.precision(std::numeric_limits::max_digits10);` during debugging may be wise because it gives some warning if construction of multiprecision values involves unintended conversion from `double` by showing trailing zero or random digits after [@http://en.cppreference.com/w/cpp/types/numeric_limits/max_digits10 max_digits10], that is 17 for `double`, digit 18... may be just noise.] The complete example code is at [@../../example/brent_minimise_example.cpp brent_minimise_example.cpp]. [h4 Implementation] This is a reasonably faithful implementation of Brent's algorithm. [h4 References] # Brent, R.P. 1973, Algorithms for Minimization without Derivatives, (Englewood Cliffs, NJ: Prentice-Hall), Chapter 5. # Numerical Recipes in C, The Art of Scientific Computing, Second Edition, William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Cambridge University Press. 1988, 1992. # An algorithm with guaranteed convergence for finding a zero of a function, R. P. Brent, The Computer Journal, Vol 44, 1971. # [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method in Wikipedia.] # Z. Zhang, An Improvement to the Brent's Method, IJEA, vol. 2, pp. 2 to 26, May 31, 2011. [@http://www.cscjournals.org/manuscript/Journals/IJEA/volume2/Issue1/IJEA-7.pdf ] # Steven A. Stage, Comments on An Improvement to the Brent's Method (and comparison of various algorithms) [@http://www.cscjournals.org/manuscript/Journals/IJEA/volume4/Issue1/IJEA-33.pdf] Stage concludes that Brent's algorithm is slow to start, but fast to finish convergence, and has good accuracy. [endsect] [/section:rebt_minima Locating Function Minima] [/ Copyright 2006, 2015, 2018 John Maddock and Paul A. Bristow. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt). ]