1 [section:brent_minima Locating Function Minima using Brent's algorithm]
3 [import ../../example/brent_minimise_example.cpp]
8 #include <boost/math/tools/minima.hpp>
12 template <class F, class T>
13 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits);
15 template <class F, class T>
16 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits, boost::uintmax_t& max_iter);
20 These two functions locate the minima of the continuous function ['f] using
21 [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method]: specifically it
22 uses quadratic interpolation to locate the minima, or if that fails, falls back to
23 a [@http://en.wikipedia.org/wiki/Golden_section_search golden-section search].
28 [[f] [The function to minimise: a function object (functor) that should be smooth over the
29 range ['\[min, max\]], with no maxima occurring in that interval.]]
30 [[min] [The lower endpoint of the range in which to search for the minima.]]
31 [[max] [The upper endpoint of the range in which to search for the minima.]]
32 [[bits] [The number of bits precision to which the minima should be found.[br]
33 Note that in principle, the minima can not be located to greater
34 accuracy than the square root of machine epsilon (for 64-bit double, sqrt(1e-16)[cong]1e-8),
35 therefore the value of ['bits] will be ignored if it's greater than half the number of bits
36 in the mantissa of T.]]
37 [[max_iter] [The maximum number of iterations to use
38 in the algorithm, if not provided the algorithm will just
39 keep on going until the minima is found.]]
44 A `pair` of type T containing the value of the abscissa at the minima and the value
45 of ['f(x)] at the minima.
47 [tip Defining BOOST_MATH_INSTRUMENT will show some parameters, for example:
51 tolerance = 1.19209289550781e-007
52 seeking minimum in range min-4 to 1.33333333333333
53 maximum iterations 18446744073709551615
58 [h4:example Brent Minimisation Example]
60 As a demonstration, we replicate this [@http://en.wikipedia.org/wiki/Brent%27s_method#Example Wikipedia example]
61 minimising the function ['y= (x+3)(x-1)[super 2]].
63 It is obvious from the equation and the plot that there is a
64 minimum at exactly one and the value of the function at one is exactly zero.
66 [tip This observation shows that an analytical or
67 [@http://en.wikipedia.org/wiki/Closed-form_expression Closed-form expression]
68 solution always beats brute-force hands-down for both speed and precision.]
70 [graph brent_test_function_1]
72 First an include is needed:
74 [brent_minimise_include_1]
76 This function is encoded in C++ as function object (functor) using `double` precision thus:
78 [brent_minimise_double_functor]
80 The Brent function is conveniently accessed through a `using` statement (noting sub-namespace `::tools`).
82 The search minimum and maximum are chosen as -4 to 4/3 (as in the Wikipedia example).
84 [tip S A Stage (reference 6) reports that the Brent algorithm is ['slow to start, but fast to converge],
85 so choosing a tight min-max range is good.]
87 For simplicity, we set the precision parameter `bits` to `std::numeric_limits<double>::digits`,
88 which is effectively the maximum possible i.e. `std::numeric_limits<double>::digits`/2.
89 Nor do we provide a maximum iterations parameter `max_iter`,
90 (perhaps unwidely), so the function will iterate until it finds a minimum.
92 [brent_minimise_double_1]
94 The resulting [@http://en.cppreference.com/w/cpp/utility/pair std::pair]
95 contains the minimum close to one and the minimum value close to zero.
97 x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018
99 The differences from the expected ['one] and ['zero] are less than the
100 uncertainty (for `double`) 1.5e-008 calculated from
101 `sqrt(std::numeric_limits<double>::digits) == 53`.
103 We can use it like this to check that the two values are close-enough to those expected,
105 using boost::math::fpc::is_close_to;
106 using boost::math::fpc::is_small;
108 double uncertainty = sqrt(std::numeric_limits<double>::digits);
109 is_close_to(1., r.first, uncertainty);
110 is_small(r.second, uncertainty);
112 x == 1 (compared to uncertainty 0.00034527) is true
113 f(x) == 0 (compared to uncertainty 0.00034527) is true
115 It is possible to make this comparison more generally with a templated function,
116 returning `true` when this criterion is met, for example:
118 [brent_minimise_close]
120 In practical applications, we might want to know how many iterations,
121 and maybe to limit iterations and
122 perhaps to trade some loss of precision for speed, for example:
124 [brent_minimise_double_2]
126 limits to a maximum of 20 iterations
127 (a reasonable estimate for this application, even for higher precision shown later).
129 The parameter `it` is updated to return the actual number of iterations
130 (so it may be useful to also keep a record of the limit in `maxit`).
132 It is neat to avoid showing insignificant digits by computing the number of decimal digits to display.
134 [brent_minimise_double_3]
136 Showing 53 bits precision with 9 decimal digits from tolerance 1.49011611938477e-008
137 x at minimum = 1, f(1) = 5.04852568e-018
139 We can also half the number of precision bits from 52 to 26.
141 [brent_minimise_double_4]
143 showing no change in the result and no change in the number of iterations, as expected.
145 It is only if we reduce the precision to a quarter, specifying only 13 precision bits
147 [brent_minimise_double_5]
149 that we reduce the number of iterations from 10 to 7 and the result significantly differing from ['one] and ['zero].
151 Showing 13 bits precision with 9 decimal digits from tolerance 0.015625
152 x at minimum = 0.9999776, f(0.9999776) = 2.0069572e-009 after 7 iterations.
154 [h5:template Templating on floating-point type]
156 If we want to switch the floating-point type, then the functor must be revised.
157 Since the functor is stateless, the easiest option is to simply make
158 `operator()` a template member function:
160 [brent_minimise_T_functor]
162 The `brent_find_minima` function can now be used in template form.
164 [brent_minimise_template_1]
166 The form shown uses the floating-point type `long double` by deduction,
167 but it is also possible to be more explicit, for example:
169 std::pair<long double, long double> r = brent_find_minima<func, long double>
170 (func(), bracket_min, bracket_max, bits, it);
172 In order to show the use of multiprecision below, it may be convenient to write a templated function to use this.
174 [brent_minimise_T_show]
176 We can use this with all built-in floating-point types, for example
178 [brent_minimise_template_fd]
180 and, on platforms that provide it, a
181 [@http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format 128-bit quad] type.
182 (See [@boost:libs/multiprecision/doc/html/boost_multiprecision/tut/floats/float128.html float128]).
184 For this optional include, the build should define the macro BOOST_HAVE_QUADMATH:
186 [brent_minimise_mp_include_1]
190 [brent_minimise_template_quad]
192 [h5:multiprecision Multiprecision]
194 If a higher precision than `double` (or `long double` if that is more precise) is required,
195 then this is easily achieved using __multiprecision with some includes from
197 [brent_minimise_mp_include_0]
201 [brent_minimise_mp_typedefs]
204 [brent_minimise_mp_1]
206 and with our show function
208 [brent_minimise_mp_2]
210 [brent_minimise_mp_output_1]
212 [brent_minimise_mp_output_2]
214 [tip One can usually rely on template argument deduction
215 to avoid specifying the verbose multiprecision types,
216 but great care in needed with the ['type of the values] provided
217 to avoid confusing the compiler.
220 [tip Using `std::cout.precision(std::numeric_limits<T>::digits10);`
221 or `std::cout.precision(std::numeric_limits<T>::max_digits10);`
222 during debugging may be wise because it gives some warning if construction of multiprecision values
223 involves unintended conversion from `double` by showing trailing zero or random digits after
224 [@http://en.cppreference.com/w/cpp/types/numeric_limits/max_digits10 max_digits10],
225 that is 17 for `double`, digit 18... may be just noise.]
227 The complete example code is at [@../../example/brent_minimise_example.cpp brent_minimise_example.cpp].
231 This is a reasonably faithful implementation of Brent's algorithm.
235 # Brent, R.P. 1973, Algorithms for Minimization without Derivatives,
236 (Englewood Cliffs, NJ: Prentice-Hall), Chapter 5.
238 # Numerical Recipes in C, The Art of Scientific Computing,
239 Second Edition, William H. Press, Saul A. Teukolsky,
240 William T. Vetterling, and Brian P. Flannery.
241 Cambridge University Press. 1988, 1992.
243 # An algorithm with guaranteed convergence for finding a zero
244 of a function, R. P. Brent, The Computer Journal, Vol 44, 1971.
246 # [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method in Wikipedia.]
248 # Z. Zhang, An Improvement to the Brent's Method, IJEA, vol. 2, pp. 2 to 26, May 31, 2011.
249 [@http://www.cscjournals.org/manuscript/Journals/IJEA/volume2/Issue1/IJEA-7.pdf ]
251 # Steven A. Stage, Comments on An Improvement to the Brent's Method
252 (and comparison of various algorithms)
253 [@http://www.cscjournals.org/manuscript/Journals/IJEA/volume4/Issue1/IJEA-33.pdf]
254 Stage concludes that Brent's algorithm is slow to start, but fast to finish convergence, and has good accuracy.
256 [endsect] [/section:rebt_minima Locating Function Minima]
259 Copyright 2006, 2015 John Maddock and Paul A. Bristow.
260 Distributed under the Boost Software License, Version 1.0.
261 (See accompanying file LICENSE_1_0.txt or copy at
262 http://www.boost.org/LICENSE_1_0.txt).