279 lines
10 KiB
Plaintext
279 lines
10 KiB
Plaintext
[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).
|
|
]
|
|
|