]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | [section:brent_minima Locating Function Minima using Brent's algorithm] |
2 | ||
3 | [import ../../example/brent_minimise_example.cpp] | |
4 | ||
5 | [h4 Synopsis] | |
6 | ||
7 | `` | |
8 | #include <boost/math/tools/minima.hpp> | |
9 | ||
10 | `` | |
11 | ||
12 | template <class F, class T> | |
13 | std::pair<T, T> brent_find_minima(F f, T min, T max, int bits); | |
14 | ||
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); | |
17 | ||
18 | [h4 Description] | |
19 | ||
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]. | |
24 | ||
25 | [*Parameters] | |
26 | ||
27 | [variablelist | |
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.]] | |
40 | ] [/variablelist] | |
41 | ||
42 | [*Returns:] | |
43 | ||
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. | |
46 | ||
47 | [tip Defining BOOST_MATH_INSTRUMENT will show some parameters, for example: | |
48 | `` | |
49 | Type T is double | |
50 | bits = 24, maximum 26 | |
51 | tolerance = 1.19209289550781e-007 | |
52 | seeking minimum in range min-4 to 1.33333333333333 | |
53 | maximum iterations 18446744073709551615 | |
54 | 10 iterations. | |
55 | `` | |
56 | ] | |
57 | ||
58 | [h4:example Brent Minimisation Example] | |
59 | ||
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]]. | |
62 | ||
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. | |
65 | ||
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.] | |
69 | ||
70 | [graph brent_test_function_1] | |
71 | ||
72 | First an include is needed: | |
73 | ||
74 | [brent_minimise_include_1] | |
75 | ||
76 | This function is encoded in C++ as function object (functor) using `double` precision thus: | |
77 | ||
78 | [brent_minimise_double_functor] | |
79 | ||
80 | The Brent function is conveniently accessed through a `using` statement (noting sub-namespace `::tools`). | |
81 | ||
82 | The search minimum and maximum are chosen as -4 to 4/3 (as in the Wikipedia example). | |
83 | ||
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.] | |
86 | ||
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. | |
91 | ||
92 | [brent_minimise_double_1] | |
93 | ||
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. | |
96 | ||
97 | x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018 | |
98 | ||
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`. | |
102 | ||
103 | We can use it like this to check that the two values are close-enough to those expected, | |
104 | ||
105 | using boost::math::fpc::is_close_to; | |
106 | using boost::math::fpc::is_small; | |
107 | ||
108 | double uncertainty = sqrt(std::numeric_limits<double>::digits); | |
109 | is_close_to(1., r.first, uncertainty); | |
110 | is_small(r.second, uncertainty); | |
111 | ||
112 | x == 1 (compared to uncertainty 0.00034527) is true | |
113 | f(x) == 0 (compared to uncertainty 0.00034527) is true | |
114 | ||
115 | It is possible to make this comparison more generally with a templated function, | |
116 | returning `true` when this criterion is met, for example: | |
117 | ||
118 | [brent_minimise_close] | |
119 | ||
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: | |
123 | ||
124 | [brent_minimise_double_2] | |
125 | ||
126 | limits to a maximum of 20 iterations | |
127 | (a reasonable estimate for this application, even for higher precision shown later). | |
128 | ||
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`). | |
131 | ||
132 | It is neat to avoid showing insignificant digits by computing the number of decimal digits to display. | |
133 | ||
134 | [brent_minimise_double_3] | |
135 | ||
136 | Showing 53 bits precision with 9 decimal digits from tolerance 1.49011611938477e-008 | |
137 | x at minimum = 1, f(1) = 5.04852568e-018 | |
138 | ||
139 | We can also half the number of precision bits from 52 to 26. | |
140 | ||
141 | [brent_minimise_double_4] | |
142 | ||
143 | showing no change in the result and no change in the number of iterations, as expected. | |
144 | ||
145 | It is only if we reduce the precision to a quarter, specifying only 13 precision bits | |
146 | ||
147 | [brent_minimise_double_5] | |
148 | ||
149 | that we reduce the number of iterations from 10 to 7 and the result significantly differing from ['one] and ['zero]. | |
150 | ||
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. | |
153 | ||
154 | [h5:template Templating on floating-point type] | |
155 | ||
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: | |
159 | ||
160 | [brent_minimise_T_functor] | |
161 | ||
162 | The `brent_find_minima` function can now be used in template form. | |
163 | ||
164 | [brent_minimise_template_1] | |
165 | ||
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: | |
168 | ||
169 | std::pair<long double, long double> r = brent_find_minima<func, long double> | |
170 | (func(), bracket_min, bracket_max, bits, it); | |
171 | ||
172 | In order to show the use of multiprecision below, it may be convenient to write a templated function to use this. | |
173 | ||
174 | [brent_minimise_T_show] | |
175 | ||
176 | We can use this with all built-in floating-point types, for example | |
177 | ||
178 | [brent_minimise_template_fd] | |
179 | ||
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]). | |
183 | ||
184 | For this optional include, the build should define the macro BOOST_HAVE_QUADMATH: | |
185 | ||
186 | [brent_minimise_mp_include_1] | |
187 | ||
188 | or | |
189 | ||
190 | [brent_minimise_template_quad] | |
191 | ||
192 | [h5:multiprecision Multiprecision] | |
193 | ||
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 | |
196 | ||
197 | [brent_minimise_mp_include_0] | |
198 | ||
199 | and some `typdef`s. | |
200 | ||
201 | [brent_minimise_mp_typedefs] | |
202 | Using thus | |
203 | ||
204 | [brent_minimise_mp_1] | |
205 | ||
206 | and with our show function | |
207 | ||
208 | [brent_minimise_mp_2] | |
209 | ||
210 | [brent_minimise_mp_output_1] | |
211 | ||
212 | [brent_minimise_mp_output_2] | |
213 | ||
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. | |
218 | ] | |
219 | ||
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.] | |
226 | ||
227 | The complete example code is at [@../../example/brent_minimise_example.cpp brent_minimise_example.cpp]. | |
228 | ||
229 | [h4 Implementation] | |
230 | ||
231 | This is a reasonably faithful implementation of Brent's algorithm. | |
232 | ||
233 | [h4 References] | |
234 | ||
235 | # Brent, R.P. 1973, Algorithms for Minimization without Derivatives, | |
236 | (Englewood Cliffs, NJ: Prentice-Hall), Chapter 5. | |
237 | ||
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. | |
242 | ||
243 | # An algorithm with guaranteed convergence for finding a zero | |
244 | of a function, R. P. Brent, The Computer Journal, Vol 44, 1971. | |
245 | ||
246 | # [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method in Wikipedia.] | |
247 | ||
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 ] | |
250 | ||
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. | |
255 | ||
256 | [endsect] [/section:rebt_minima Locating Function Minima] | |
257 | ||
258 | [/ | |
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). | |
263 | ] | |
264 |