EQ2EMu/EQ2/source/depends/boost_1_72_0/libs/math/doc/roots/minima.qbk

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