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 pair of type T containing the value of the abscissa 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 one and the value of the function at one is exactly zero.

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

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 i.e. std::numeric_limits<double>::digits/2. Nor do we provide a maximum iterations parameter max_iter, (perhaps unwidely), so the function will iterate until it finds a minimum.

int bits = std::numeric_limits<double>::digits;

std::pair<double, double> r = brent_find_minima(funcdouble(), -4., 4. / 3, bits);

std::cout.precision(std::numeric_limits<double>::digits10);
std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << std::endl;
// x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018

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>::digits) == 53.

We can use it like this to check that the two values are close-enough to those expected,

 using boost::math::fpc::is_close_to;
 using boost::math::fpc::is_small;

double uncertainty = sqrt(std::numeric_limits<double>::digits);
is_close_to(1., r.first, uncertainty);
is_small(r.second, uncertainty);

x == 1 (compared to uncertainty 0.00034527) is true
f(x) == 0 (compared to uncertainty 0.00034527) is true

It is possible to make this comparison more generally with a templated function, returning true when this criterion is met, for example:

//
template <class T = double>
bool close(T expect, T got, T tolerance)
{
  using boost::math::fpc::is_close_to;
  using boost::math::fpc::is_small;

  if (is_small<T>(expect, tolerance))
  {
    return is_small<T>(got, tolerance);
  }
  else
  {
    return is_close_to<T>(expect, got, tolerance);
  }
}

In practical applications, we might want to know how many iterations, and maybe to limit iterations and perhaps to trade some loss of precision for speed, for example:

const boost::uintmax_t maxit = 20;
boost::uintmax_t it = maxit;
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 application, even for 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 in 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(bits));  // Number of significant decimal digits.
std::cout << "Showing " << bits << " bits precision with " << prec
  << " decimal digits from tolerance " << sqrt(std::numeric_limits<double>::epsilon())
  << std::endl;
std::streamsize precision = std::cout.precision(prec); // Save.

std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second
  << " after " << it << " iterations. " << std::endl;

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.

bits /= 2; // Half digits precision (effective maximum).
double epsilon_2 = boost::math::pow<-(std::numeric_limits<double>::digits/2 - 1), double>(2);

std::cout << "Showing " << bits << " bits precision with " << prec
  << " decimal digits from tolerance " << sqrt(epsilon_2)
  << std::endl;
std::streamsize precision = std::cout.precision(prec); // Save.

boost::uintmax_t it = maxit;
r = brent_find_minima(funcdouble(), -4., 4. / 3, bits, it);
std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << std::endl;
std::cout << it << " iterations. " << std::endl;

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

bits /= 2; // Quarter precision.
double epsilon_4 = boost::math::pow<-(std::numeric_limits<double>::digits / 4 - 1), double>(2);

std::cout << "Showing " << bits << " bits precision with " << prec
  << " decimal digits from tolerance " << sqrt(epsilon_4)
  << std::endl;
std::streamsize precision = std::cout.precision(prec); // Save.

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

that we reduce the number of iterations from 10 to 7 and the result significantly differing 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); //
  }
};

The brent_find_minima function can now be used in template form.

std::cout.precision(std::numeric_limits<long double>::digits10);
long double bracket_min = -4.;
long double bracket_max = 4. / 3;
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;

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, it may be convenient to write a templated function to use this.

template <class T>
void show_minima()
{
  using boost::math::tools::brent_find_minima;
  try
  { // Always use try'n'catch blocks with Boost.Math to get any error messages.

    int bits = std::numeric_limits<T>::digits/2; // Maximum is digits/2;
    std::streamsize prec = static_cast<int>(2 + sqrt(bits));  // Number of significant decimal digits.
    std::streamsize precision = std::cout.precision(prec); // Save.

    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 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 " << close(static_cast<T>(1), r.first, uncertainty) << std::endl;
    //std::cout << std::boolalpha << "f(x) == (0 compared to uncertainty " << uncertainty << ") 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()

We can use this with all built-in floating-point types, for example

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

and, on platforms that provide it, a 128-bit quad type. (See float128).

For this optional include, the build should define the macro BOOST_HAVE_QUADMATH:

#ifdef BOOST_HAVE_QUADMATH  // Define only if GCC or Intel and have quadmath.lib or .dll library available.
  using boost::multiprecision::float128;
#endif

or

// #ifndef _MSC_VER
#ifdef BOOST_HAVE_QUADMATH  // Define 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 from

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

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;

using boost::multiprecision::cpp_dec_float_50; // decimal.

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;

Using thus

std::cout.precision(std::numeric_limits<cpp_bin_float_50>::digits10);

cpp_bin_float_50 fpv("-1.2345");
cpp_bin_float_50 absv;

absv = fpv < static_cast<cpp_bin_float_50>(0) ? -fpv : fpv;
std::cout << fpv << ' ' << absv << std::endl;


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 << bracket_min << " " << bracket_max << std::endl;
const boost::uintmax_t maxit = 20;
boost::uintmax_t it = maxit;
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 << ", f(" << r.first << ") = " << r.second
// x at minimum = 1, f(1) = 5.04853e-018
  << ", after " << it << " iterations. " << std::endl;

close(static_cast<cpp_bin_float_50>(1), r.first, sqrt(std::numeric_limits<cpp_bin_float_50>::epsilon()));

and with our show 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.

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