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

Эти две функции обнаруживают минимум непрерывной функцииfс помощью методаБрента: В частности, он использует квадратичную интерполяцию, чтобы найти минимум, или, если это не удается, возвращается к поиску золотого сечения.

Параметры

Возвращение:

A<pair>типа T, содержащего значение абсциссы при минимуме и значениеf(x)при минимуме.

	Tip
	Обсуждение Boost_MATH_INSTRUMENT Покажут некоторые параметры, например: 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

В качестве демонстрации мы воспроизводим этотпример Википедии, минимизируя функциюy= (x+3) (x-1)².

Из уравнения и графика очевидно, что существует минимум в одном, а значение функции в одном равно нулю.

	Tip
	Это наблюдение показывает, что аналитический илиРешение для выражения в закрытой формевсегда превосходит грубую силу как для скорости, так и для точности.

Во-первых, необходимо включить:

#include <boost/math/tools/minima.hpp>

Эта функция кодируется на C++ в качестве функционального объекта (функтора) с использованием<double>точности таким образом:

struct funcdouble
{
  double operator()(double const& x)
  { //
    return (x + 3) * (x - 1) * (x - 1); // (x + 3)(x - 1)^2
  }
};

Функция Brent легко доступна через утверждение<using>(с указанием подпространства имен<::tools>).

Минимальный и максимальный поисковые запросы выбираются как -4 до 4/3 (как в примере Википедии).

	Tip
	S A Stage (ссылка 6) сообщает, что алгоритм Brentмедленно запускается, но быстро сходится, поэтому выбор узкого диапазона min-max хорош.

Для простоты мы устанавливаем прецизионный параметр<bits>на<std::numeric_limits<double>::digits>, что фактически является максимально возможным, т.е.<std::numeric_limits<double>::digits>/2. Также мы не предоставляем максимальный параметр итераций<max_iter>, (возможно, нешироко), поэтому функция будет повторяться, пока не найдет минимум.

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

Полученныйstd::pairсодержит минимум, близкий к одному, и минимальное значение, близкое к нулю.

x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018

Отличия от ожидаемогоодногоинуляменьше неопределенности (для<double>1,5e-008, рассчитанной из<sqrt(std::numeric_limits<double>::digits)==53>.

Мы можем использовать его таким образом, чтобы проверить, что эти два значения достаточно близки к ожидаемым.

 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

Можно сделать это сравнение более общим с шаблонной функцией, возвращаясь<true>при выполнении этого критерия, например:

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

В практическом применении мы могли бы захотеть узнать, сколько итераций, и, возможно, ограничить итерации и, возможно, обменять некоторую потерю точности на скорость, например:

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;

предел до 20 итераций (разумная оценка для этого приложения, даже для более высокой точности, показанной позже).

Параметр<it>обновляется, чтобы вернуть фактическое количество итераций (так что может быть полезно также вести запись предела в<maxit>).

Опрятно избегать показа незначительных цифр, вычисляя число десятичных цифр для отображения.

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

Мы также можем вдвое уменьшить количество прецизионных битов с 52 до 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;

не показывает никаких изменений в результате и никаких изменений в количестве итераций, как ожидалось.

Только если мы уменьшим точность до четверти, указывая только 13 прецизионных битов.

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;

Мы уменьшаем число итераций с 10 до 7, и результат значительно отличается отодногоинуля.

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

Если мы хотим переключить тип с плавающей точкой, то функтор должен быть пересмотрен. Поскольку функтор не имеет состояния, проще всего сделать<operator()>функцию члена шаблона:

struct func
{
  template <class T>
  T operator()(T const& x)
  { //
    return (x + 3) * (x - 1) * (x - 1); //
  }
};

Функция<brent_find_minima>теперь может использоваться в форме шаблона.

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;

Показываемая форма использует тип плавающей точки<long double>путем вычета, но также можно быть более явным, например:

std::pair<long double, long double> r = brent_find_minima<func, long double>
(func(), bracket_min, bracket_max, bits, it);

Чтобы показать использование многоточности ниже, может быть удобно написать шаблонную функцию для использования этого.

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

Мы можем использовать это со всеми встроенными типами плавающих точек.

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

и, на платформах, которые обеспечивают его,128-битный квадтип (см.float128).

Для этого опционального включения сборка должна определять макрос 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

или

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

Если требуется более высокая точность, чем<double>(или<longdouble>, если это более точно), то это легко достигается с помощьюBoost.Multiprecision.

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

и некоторые из них<typdef>.

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;

Используя таким образом

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

и с нашей функцией шоу

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
	Обычно можно полагаться на вычисление аргументов шаблона, чтобы избежать указания многоточных типов глаголов, но при этом следует проявлять большую осторожность с типом значений, чтобы не путать компилятор.

Tip

Использование<std::cout.precision(std::numeric_limits<T>::digits10);>или<std::cout.precision(std::numeric_limits<T>::max_digits10);>во время отладки может быть разумным, поскольку оно дает некоторое предупреждение, если построение многоточных значений включает непреднамеренное преобразование из<double>, показывая нулевые или случайные цифры послеmax_digits10, то есть 17 для<double>, цифра 18 ... может быть просто шумом.

Полный пример кода находится по адресуbrent_minimise_example.cpp.

Implementation

Это достаточно надежная реализация алгоритма Брента.

References