...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#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, std::uintmax_t& max_iter);
These two functions locate the minima of the continuous function f using Brent's method: specifically it uses quadratic interpolation to locate the minima, or if that fails, falls back to a golden-section search.
Parameters
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.
The lower endpoint of the range in which to search for the minima.
The upper endpoint of the range in which to search for the minima.
The number of bits precision to which the minima should be found.
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)≅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.
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.
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. |
As a demonstration, we replicate this Wikipedia example minimising the function y= (x+3)(x-1)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 Closed-form expression solution always beats brute-force hands-down for both speed and precision. |
First an include is needed:
#include <boost/math/tools/minima.hpp>
This function is encoded in C++ as function object (functor) using double
precision thus:
struct funcdouble { double operator()(double const& x) { return (x + 3) * (x - 1) * (x - 1); // (x + 3)(x - 1)^2 } };
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.
const int double_bits = std::numeric_limits<double>::digits; std::pair<double, double> r = brent_find_minima(funcdouble(), -4., 4. / 3, double_bits); std::streamsize precision_1 = std::cout.precision(std::numeric_limits<double>::digits10); // Show all double precision decimal digits and trailing zeros. std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << std::endl;
The resulting 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,
using boost::math::fpc::close_at_tolerance; using boost::math::fpc::is_small; std::cout << "x = " << r.first << ", f(x) = " << r.second << std::endl; std::cout << std::boolalpha << "x == 1 (compared to uncertainty " << uncertainty << ") is " << is_close(1., r.first, uncertainty) << std::endl; // true std::cout << std::boolalpha << "f(x) == 0 (compared to uncertainty " << uncertainty << ") is " << is_close(0., r.second, uncertainty) << std::endl; // true
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 |
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:
//! Compare if value got is close to expected, //! checking first if expected is very small //! (to avoid divide by tiny or zero during comparison) //! before comparing expect with value got. template <class T> bool is_close(T expect, T got, T tolerance) { using boost::math::fpc::close_at_tolerance; using boost::math::fpc::is_small; using boost::math::fpc::FPC_STRONG; if (is_small<T>(expect, tolerance)) { return is_small<T>(got, tolerance); } return close_at_tolerance<T>(tolerance, FPC_STRONG) (expect, got); } // bool is_close(T expect, T got, T tolerance)
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:
const std::uintmax_t maxit = 20; std::uintmax_t it = maxit; std::pair<double, double> r = brent_find_minima(funcdouble(), -4., 4. / 3, bits, it); std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << " after " << it << " iterations. " << std::endl;
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.
std::streamsize prec = static_cast<int>(2 + sqrt((double)bits)); // Number of significant decimal digits. std::streamsize precision_3 = std::cout.precision(prec); // Save and set new precision. std::cout << "Showing " << bits << " bits " "precision with " << prec << " decimal digits from tolerance " << sqrt(std::numeric_limits<double>::epsilon()) << std::endl; std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << " after " << it << " iterations. " << std::endl; std::cout.precision(precision_3); // Restore.
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:
const int bits_div_2 = std::numeric_limits<double>::digits / 2; // Half digits precision (effective maximum). double epsilon_2 = boost::math::pow<-(std::numeric_limits<double>::digits/2 - 1), double>(2); std::streamsize prec = static_cast<int>(2 + sqrt((double)bits_div_2)); // Number of significant decimal digits. std::cout << "Showing " << bits_div_2 << " bits precision with " << prec << " decimal digits from tolerance " << sqrt(epsilon_2) << std::endl; std::streamsize precision_4 = std::cout.precision(prec); // Save. const std::uintmax_t maxit = 20; std::uintmax_t it_4 = maxit; std::pair<double, double> r = brent_find_minima(funcdouble(), -4., 4. / 3, bits_div_2, it_4); std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << std::endl; std::cout << it_4 << " iterations. " << std::endl; std::cout.precision(precision_4); // Restore.
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
const int bits_div_4 = std::numeric_limits<double>::digits / 4; // Quarter precision. double epsilon_4 = boost::math::pow<-(std::numeric_limits<double>::digits / 4 - 1), double>(2); std::streamsize prec = static_cast<int>(2 + sqrt((double)bits_div_4)); // Number of significant decimal digits. std::cout << "Showing " << bits_div_4 << " bits precision with " << prec << " decimal digits from tolerance " << sqrt(epsilon_4) << std::endl; std::streamsize precision_5 = std::cout.precision(prec); // Save & set. const std::uintmax_t maxit = 20; std::uintmax_t it_5 = maxit; std::pair<double, double> r = brent_find_minima(funcdouble(), -4., 4. / 3, bits_div_4, it_5); std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << ", after " << it_5 << " iterations. " << std::endl; std::cout.precision(precision_5); // Restore.
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.
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:
struct func { template <class T> T operator()(T const& x) { return (x + 3) * (x - 1) * (x - 1); // (x + 3)(x - 1)^2 } };
The brent_find_minima
function
can now be used in template form, for example using long
double
:
std::streamsize precision_t1 = std::cout.precision(std::numeric_limits<long double>::digits10); // Save & set. long double bracket_min = -4.; long double bracket_max = 4. / 3; const int bits = std::numeric_limits<long double>::digits; const std::uintmax_t maxit = 20; std::uintmax_t it = maxit; std::pair<long double, long double> r = brent_find_minima(func(), bracket_min, bracket_max, bits, it); std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << ", after " << it << " iterations. " << std::endl; std::cout.precision(precision_t1); // Restore.
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:
//! Example template function to find and show minima. //! \tparam T floating-point or fixed_point type. template <class T> void show_minima() { using boost::math::tools::brent_find_minima; using std::sqrt; try { // Always use try'n'catch blocks with Boost.Math to ensure you get any error messages. int bits = std::numeric_limits<T>::digits/2; // Maximum is digits/2; std::streamsize prec = static_cast<int>(2 + sqrt((double)bits)); // Number of significant decimal digits. std::streamsize precision = std::cout.precision(prec); // Save and set. std::cout << "\n\nFor type: " << typeid(T).name() << ",\n epsilon = " << std::numeric_limits<T>::epsilon() // << ", precision of " << bits << " bits" << ",\n the maximum theoretical precision from Brent's minimization is " << sqrt(std::numeric_limits<T>::epsilon()) << "\n Displaying to std::numeric_limits<T>::digits10 " << prec << ", significant decimal digits." << std::endl; const std::uintmax_t maxit = 20; std::uintmax_t it = maxit; // Construct using string, not double, avoids loss of precision. //T bracket_min = static_cast<T>("-4"); //T bracket_max = static_cast<T>("1.3333333333333333333333333333333333333333333333333"); // Construction from double may cause loss of precision for multiprecision types like cpp_bin_float, // but brackets values are good enough for using Brent minimization. T bracket_min = static_cast<T>(-4); T bracket_max = static_cast<T>(1.3333333333333333333333333333333333333333333333333); std::pair<T, T> r = brent_find_minima<func, T>(func(), bracket_min, bracket_max, bits, it); std::cout << " x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second; if (it < maxit) { std::cout << ",\n met " << bits << " bits precision" << ", after " << it << " iterations." << std::endl; } else { std::cout << ",\n did NOT meet " << bits << " bits precision" << " after " << it << " iterations!" << std::endl; } // Check that result is that expected (compared to theoretical uncertainty). T uncertainty = sqrt(std::numeric_limits<T>::epsilon()); std::cout << std::boolalpha << "x == 1 (compared to uncertainty " << uncertainty << ") is " << is_close(static_cast<T>(1), r.first, uncertainty) << std::endl; std::cout << std::boolalpha << "f(x) == (0 compared to uncertainty " << uncertainty << ") is " << is_close(static_cast<T>(0), r.second, uncertainty) << std::endl; // Problems with this using multiprecision with expression template on? std::cout.precision(precision); // Restore. } catch (const std::exception& e) { // Always useful to include try & catch blocks because default policies // are to throw exceptions on arguments that cause errors like underflow, overflow. // Lacking try & catch blocks, the program will abort without a message below, // which may give some helpful clues as to the cause of the exception. std::cout << "\n""Message from thrown exception was:\n " << e.what() << std::endl; } } // void show_minima()
Note | |
---|---|
the prudent addition of |
We can use this with all built-in floating-point types, for example
show_minima<float>(); show_minima<double>(); show_minima<long double>();
On platforms that provide it, a 128-bit quad type can be used. (See float128).
If this type is available, the build should define the macro BOOST_HAVE_QUADMATH:
#ifdef BOOST_HAVE_QUADMATH // Defined only if GCC or Intel and have quadmath.lib or .dll library available. using boost::multiprecision::float128; #endif
and can be (conditionally) used this:
#ifdef BOOST_HAVE_QUADMATH // Defined only if GCC or Intel and have quadmath.lib or .dll library available. show_minima<float128>(); // Needs quadmath_snprintf, sqrtQ, fabsq that are in in quadmath library. #endif
If a higher precision than double
(or long double
if that is more precise) is required, then this is easily achieved using Boost.Multiprecision
with some includes, for example:
#include <boost/multiprecision/cpp_dec_float.hpp> // For decimal boost::multiprecision::cpp_dec_float_50. #include <boost/multiprecision/cpp_bin_float.hpp> // For binary boost::multiprecision::cpp_bin_float_50;
and some typdef
s.
using boost::multiprecision::cpp_bin_float_50; // binary multiprecision typedef. using boost::multiprecision::cpp_dec_float_50; // decimal multiprecision typedef. // One might also need typedefs like these to switch expression templates off and on (default is on). typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<50>, boost::multiprecision::et_on> cpp_bin_float_50_et_on; // et_on is default so is same as cpp_bin_float_50. typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<50>, boost::multiprecision::et_off> cpp_bin_float_50_et_off; typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<50>, boost::multiprecision::et_on> // et_on is default so is same as cpp_dec_float_50. cpp_dec_float_50_et_on; typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<50>, boost::multiprecision::et_off> cpp_dec_float_50_et_off;
Used thus
std::cout.precision(std::numeric_limits<cpp_bin_float_50>::digits10); int bits = std::numeric_limits<cpp_bin_float_50>::digits / 2 - 2; cpp_bin_float_50 bracket_min = static_cast<cpp_bin_float_50>("-4"); cpp_bin_float_50 bracket_max = static_cast<cpp_bin_float_50>("1.3333333333333333333333333333333333333333333333333"); std::cout << "Bracketing " << bracket_min << " to " << bracket_max << std::endl; const std::uintmax_t maxit = 20; std::uintmax_t it = maxit; // Will be updated with actual iteration count. std::pair<cpp_bin_float_50, cpp_bin_float_50> r = brent_find_minima(func(), bracket_min, bracket_max, bits, it); std::cout << "x at minimum = " << r.first << ",\n f(" << r.first << ") = " << r.second // x at minimum = 1, f(1) = 5.04853e-018 << ", after " << it << " iterations. " << std::endl; is_close_to(static_cast<cpp_bin_float_50>("1"), r.first, sqrt(std::numeric_limits<cpp_bin_float_50>::epsilon())); is_close_to(static_cast<cpp_bin_float_50>("0"), r.second, sqrt(std::numeric_limits<cpp_bin_float_50>::epsilon()));
and with our show_minima
function
show_minima<cpp_bin_float_50_et_on>(); //
For type class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,1>, epsilon = 5.3455294202e-51, the maximum theoretical precision from Brent minimization is 7.311312755e-26 Displaying to std::numeric_limits<T>::digits10 11 significant decimal digits. x at minimum = 1, f(1) = 5.6273022713e-58, met 84 bits precision, after 14 iterations. x == 1 (compared to uncertainty 7.311312755e-26) is true f(x) == (0 compared to uncertainty 7.311312755e-26) is true -4 1.3333333333333333333333333333333333333333333333333 x at minimum = 0.99999999999999999999999999998813903221565569205253, f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58 14 iterations
For type class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50, 10, void, int, 0, 0>, 1>,
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 |
The complete example code is at brent_minimise_example.cpp.
This is a reasonably faithful implementation of Brent's algorithm.