123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278 |
- [section:brent_minima Locating Function Minima using Brent's algorithm]
- [import ../../example/brent_minimise_example.cpp]
- [h4 Synopsis]
- ``
- #include <boost/math/tools/minima.hpp>
- ``
- template <class F, class T>
- std::pair<T, T> brent_find_minima(F f, T min, T max, int bits);
- template <class F, class T>
- std::pair<T, T> 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<double>::digits`,
- which is effectively the maximum possible `std::numeric_limits<double>::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<double>::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<long double, long double> r = brent_find_minima<func, long double>
- (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<T>::digits10);`
- or `std::cout.precision(std::numeric_limits<T>::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).
- ]
|