1 @cindex optimization, see minimization
2 @cindex maximization, see minimization
3 @cindex minimization, one-dimensional
5 @cindex nonlinear functions, minimization
7 This chapter describes routines for finding minima of arbitrary
8 one-dimensional functions. The library provides low level components
9 for a variety of iterative minimizers and convergence tests. These can be
10 combined by the user to achieve the desired solution, with full access
11 to the intermediate steps of the algorithms. Each class of methods uses
12 the same framework, so that you can switch between minimizers at runtime
13 without needing to recompile your program. Each instance of a minimizer
14 keeps track of its own state, allowing the minimizers to be used in
15 multi-threaded programs.
17 The header file @file{gsl_min.h} contains prototypes for the
18 minimization functions and related declarations. To use the minimization
19 algorithms to find the maximum of a function simply invert its sign.
22 * Minimization Overview::
23 * Minimization Caveats::
24 * Initializing the Minimizer::
25 * Providing the function to minimize::
26 * Minimization Iteration::
27 * Minimization Stopping Parameters::
28 * Minimization Algorithms::
29 * Minimization Examples::
30 * Minimization References and Further Reading::
33 @node Minimization Overview
35 @cindex minimization, overview
37 The minimization algorithms begin with a bounded region known to contain
38 a minimum. The region is described by a lower bound @math{a} and an
39 upper bound @math{b}, with an estimate of the location of the minimum
44 @center @image{min-interval,3.4in}
48 The value of the function at @math{x} must be less than the value of the
49 function at the ends of the interval,
51 $$f(a) > f(x) < f(b)$$
61 This condition guarantees that a minimum is contained somewhere within
62 the interval. On each iteration a new point @math{x'} is selected using
63 one of the available algorithms. If the new point is a better estimate
64 of the minimum, i.e.@: where @math{f(x') < f(x)}, then the current
65 estimate of the minimum @math{x} is updated. The new point also allows
66 the size of the bounded interval to be reduced, by choosing the most
67 compact set of points which satisfies the constraint @math{f(a) > f(x) <
68 f(b)}. The interval is reduced until it encloses the true minimum to a
69 desired tolerance. This provides a best estimate of the location of the
70 minimum and a rigorous error estimate.
72 Several bracketing algorithms are available within a single framework.
73 The user provides a high-level driver for the algorithm, and the
74 library provides the individual functions necessary for each of the
75 steps. There are three main phases of the iteration. The steps are,
79 initialize minimizer state, @var{s}, for algorithm @var{T}
82 update @var{s} using the iteration @var{T}
85 test @var{s} for convergence, and repeat iteration if necessary
89 The state for the minimizers is held in a @code{gsl_min_fminimizer}
90 struct. The updating procedure uses only function evaluations (not
93 @node Minimization Caveats
95 @cindex minimization, caveats
97 Note that minimization functions can only search for one minimum at a
98 time. When there are several minima in the search area, the first
99 minimum to be found will be returned; however it is difficult to predict
100 which of the minima this will be. @emph{In most cases, no error will be
101 reported if you try to find a minimum in an area where there is more
104 With all minimization algorithms it can be difficult to determine the
105 location of the minimum to full numerical precision. The behavior of the
106 function in the region of the minimum @math{x^*} can be approximated by
110 y = f(x^*) + {1 \over 2} f''(x^*) (x - x^*)^2
116 y = f(x^*) + (1/2) f''(x^*) (x - x^*)^2
121 and the second term of this expansion can be lost when added to the
122 first term at finite precision. This magnifies the error in locating
123 @math{x^*}, making it proportional to @math{\sqrt \epsilon} (where
124 @math{\epsilon} is the relative accuracy of the floating point numbers).
125 For functions with higher order minima, such as @math{x^4}, the
126 magnification of the error is correspondingly worse. The best that can
127 be achieved is to converge to the limit of numerical accuracy in the
128 function values, rather than the location of the minimum itself.
130 @node Initializing the Minimizer
131 @section Initializing the Minimizer
133 @deftypefun {gsl_min_fminimizer *} gsl_min_fminimizer_alloc (const gsl_min_fminimizer_type * @var{T})
134 This function returns a pointer to a newly allocated instance of a
135 minimizer of type @var{T}. For example, the following code
136 creates an instance of a golden section minimizer,
139 const gsl_min_fminimizer_type * T
140 = gsl_min_fminimizer_goldensection;
141 gsl_min_fminimizer * s
142 = gsl_min_fminimizer_alloc (T);
145 If there is insufficient memory to create the minimizer then the function
146 returns a null pointer and the error handler is invoked with an error
147 code of @code{GSL_ENOMEM}.
150 @deftypefun int gsl_min_fminimizer_set (gsl_min_fminimizer * @var{s}, gsl_function * @var{f}, double @var{x_minimum}, double @var{x_lower}, double @var{x_upper})
151 This function sets, or resets, an existing minimizer @var{s} to use the
152 function @var{f} and the initial search interval [@var{x_lower},
153 @var{x_upper}], with a guess for the location of the minimum
156 If the interval given does not contain a minimum, then the function
157 returns an error code of @code{GSL_EINVAL}.
160 @deftypefun int gsl_min_fminimizer_set_with_values (gsl_min_fminimizer * @var{s}, gsl_function * @var{f}, double @var{x_minimum}, double @var{f_minimum}, double @var{x_lower}, double @var{f_lower}, double @var{x_upper}, double @var{f_upper})
161 This function is equivalent to @code{gsl_min_fminimizer_set} but uses
162 the values @var{f_minimum}, @var{f_lower} and @var{f_upper} instead of
163 computing @code{f(x_minimum)}, @code{f(x_lower)} and @code{f(x_upper)}.
167 @deftypefun void gsl_min_fminimizer_free (gsl_min_fminimizer * @var{s})
168 This function frees all the memory associated with the minimizer
172 @deftypefun {const char *} gsl_min_fminimizer_name (const gsl_min_fminimizer * @var{s})
173 This function returns a pointer to the name of the minimizer. For example,
176 printf ("s is a '%s' minimizer\n",
177 gsl_min_fminimizer_name (s));
181 would print something like @code{s is a 'brent' minimizer}.
184 @node Providing the function to minimize
185 @section Providing the function to minimize
186 @cindex minimization, providing a function to minimize
188 You must provide a continuous function of one variable for the
189 minimizers to operate on. In order to allow for general parameters the
190 functions are defined by a @code{gsl_function} data type
191 (@pxref{Providing the function to solve}).
193 @node Minimization Iteration
196 The following functions drive the iteration of each algorithm. Each
197 function performs one iteration to update the state of any minimizer of the
198 corresponding type. The same functions work for all minimizers so that
199 different methods can be substituted at runtime without modifications to
202 @deftypefun int gsl_min_fminimizer_iterate (gsl_min_fminimizer * @var{s})
203 This function performs a single iteration of the minimizer @var{s}. If the
204 iteration encounters an unexpected problem then an error code will be
209 the iteration encountered a singular point where the function evaluated
210 to @code{Inf} or @code{NaN}.
213 the algorithm could not improve the current best approximation or
218 The minimizer maintains a current best estimate of the position of the
219 minimum at all times, and the current interval bounding the minimum.
220 This information can be accessed with the following auxiliary functions,
222 @deftypefun double gsl_min_fminimizer_x_minimum (const gsl_min_fminimizer * @var{s})
223 This function returns the current estimate of the position of the
224 minimum for the minimizer @var{s}.
227 @deftypefun double gsl_min_fminimizer_x_upper (const gsl_min_fminimizer * @var{s})
228 @deftypefunx double gsl_min_fminimizer_x_lower (const gsl_min_fminimizer * @var{s})
229 These functions return the current upper and lower bound of the interval
230 for the minimizer @var{s}.
233 @deftypefun double gsl_min_fminimizer_f_minimum (const gsl_min_fminimizer * @var{s})
234 @deftypefunx double gsl_min_fminimizer_f_upper (const gsl_min_fminimizer * @var{s})
235 @deftypefunx double gsl_min_fminimizer_f_lower (const gsl_min_fminimizer * @var{s})
236 These functions return the value of the function at the current estimate
237 of the minimum and at the upper and lower bounds of the interval for the
241 @node Minimization Stopping Parameters
242 @section Stopping Parameters
243 @cindex minimization, stopping parameters
245 A minimization procedure should stop when one of the following
250 A minimum has been found to within the user-specified precision.
253 A user-specified maximum number of iterations has been reached.
256 An error has occurred.
260 The handling of these conditions is under user control. The function
261 below allows the user to test the precision of the current result.
263 @deftypefun int gsl_min_test_interval (double @var{x_lower}, double @var{x_upper}, double @var{epsabs}, double @var{epsrel})
264 This function tests for the convergence of the interval [@var{x_lower},
265 @var{x_upper}] with absolute error @var{epsabs} and relative error
266 @var{epsrel}. The test returns @code{GSL_SUCCESS} if the following
267 condition is achieved,
271 |a - b| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, \min(|a|,|b|)
278 |a - b| < epsabs + epsrel min(|a|,|b|)
283 when the interval @math{x = [a,b]} does not include the origin. If the
284 interval includes the origin then @math{\min(|a|,|b|)} is replaced by
285 zero (which is the minimum value of @math{|x|} over the interval). This
286 ensures that the relative error is accurately estimated for minima close
289 This condition on the interval also implies that any estimate of the
290 minimum @math{x_m} in the interval satisfies the same condition with respect
291 to the true minimum @math{x_m^*},
295 |x_m - x_m^*| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, x_m^*
302 |x_m - x_m^*| < epsabs + epsrel x_m^*
307 assuming that the true minimum @math{x_m^*} is contained within the interval.
310 @comment ============================================================
312 @node Minimization Algorithms
313 @section Minimization Algorithms
315 The minimization algorithms described in this section require an initial
316 interval which is guaranteed to contain a minimum---if @math{a} and
317 @math{b} are the endpoints of the interval and @math{x} is an estimate
318 of the minimum then @math{f(a) > f(x) < f(b)}. This ensures that the
319 function has at least one minimum somewhere in the interval. If a valid
320 initial interval is used then these algorithm cannot fail, provided the
321 function is well-behaved.
323 @deffn {Minimizer} gsl_min_fminimizer_goldensection
325 @cindex golden section algorithm for finding minima
326 @cindex minimum finding, golden section algorithm
328 The @dfn{golden section algorithm} is the simplest method of bracketing
329 the minimum of a function. It is the slowest algorithm provided by the
330 library, with linear convergence.
332 On each iteration, the algorithm first compares the subintervals from
333 the endpoints to the current minimum. The larger subinterval is divided
334 in a golden section (using the famous ratio @math{(3-\sqrt 5)/2 =
335 0.3189660}@dots{}) and the value of the function at this new point is
336 calculated. The new value is used with the constraint @math{f(a') >
337 f(x') < f(b')} to a select new interval containing the minimum, by
338 discarding the least useful point. This procedure can be continued
339 indefinitely until the interval is sufficiently small. Choosing the
340 golden section as the bisection ratio can be shown to provide the
341 fastest convergence for this type of algorithm.
345 @comment ============================================================
347 @deffn {Minimizer} gsl_min_fminimizer_brent
348 @cindex Brent's method for finding minima
349 @cindex minimum finding, Brent's method
351 The @dfn{Brent minimization algorithm} combines a parabolic
352 interpolation with the golden section algorithm. This produces a fast
353 algorithm which is still robust.
355 The outline of the algorithm can be summarized as follows: on each
356 iteration Brent's method approximates the function using an
357 interpolating parabola through three existing points. The minimum of the
358 parabola is taken as a guess for the minimum. If it lies within the
359 bounds of the current interval then the interpolating point is accepted,
360 and used to generate a smaller interval. If the interpolating point is
361 not accepted then the algorithm falls back to an ordinary golden section
362 step. The full details of Brent's method include some additional checks
363 to improve convergence.
366 @comment ============================================================
368 @node Minimization Examples
371 The following program uses the Brent algorithm to find the minimum of
372 the function @math{f(x) = \cos(x) + 1}, which occurs at @math{x = \pi}.
373 The starting interval is @math{(0,6)}, with an initial guess for the
377 @verbatiminclude examples/min.c
381 Here are the results of the minimization procedure.
385 @verbatiminclude examples/min.out
388 @node Minimization References and Further Reading
389 @section References and Further Reading
391 Further information on Brent's algorithm is available in the following
396 Richard Brent, @cite{Algorithms for minimization without derivatives},
397 Prentice-Hall (1973), republished by Dover in paperback (2002), ISBN