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);

Description

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

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.

Brent Minimisation Example

As a demonstration, we replicate this Wikipedia example minimising the function y= (x+3)(x-1)².

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 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:

//! 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 boost::uintmax_t maxit = 20;
boost::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 boost::uintmax_t maxit = 20;
boost::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 boost::uintmax_t maxit = 20;

boost::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.

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:

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 boost::uintmax_t maxit = 20;
boost::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 boost::uintmax_t maxit = 20;
    boost::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 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

show_minima<float>();
show_minima<double>();
show_minima<long double>();

Quad 128-bit precision

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

Multiprecision

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 typdefs.

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 boost::uintmax_t maxit = 20;
boost::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 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 max_digits10, that is 17 for double, digit 18... may be just noise.

The complete example code is at brent_minimise_example.cpp.

Implementation

This is a reasonably faithful implementation of Brent's algorithm.

References