3
<!-- This HTML file has been created by texi2html 1.54+ (gsl)
4
from ../gsl-ref.texi -->
6
<TITLE>GNU Scientific Library -- Reference Manual - One dimensional Minimization</TITLE>
7
<!-- <LINK rel="stylesheet" title="Default Style Sheet" href="/css/texinfo.css" type="text/css"> -->
8
<link href="gsl-ref_34.html" rel=Next>
9
<link href="gsl-ref_32.html" rel=Previous>
10
<link href="gsl-ref_toc.html" rel=ToC>
14
<p>Go to the <A HREF="gsl-ref_1.html">first</A>, <A HREF="gsl-ref_32.html">previous</A>, <A HREF="gsl-ref_34.html">next</A>, <A HREF="gsl-ref_50.html">last</A> section, <A HREF="gsl-ref_toc.html">table of contents</A>.
18
<H1><A NAME="SEC440" HREF="gsl-ref_toc.html#TOC440">One dimensional Minimization</A></H1>
20
<A NAME="IDX2113"></A>
21
<A NAME="IDX2114"></A>
22
<A NAME="IDX2115"></A>
23
<A NAME="IDX2116"></A>
24
<A NAME="IDX2117"></A>
28
This chapter describes routines for finding minima of arbitrary
29
one-dimensional functions. The library provides low level components
30
for a variety of iterative minimizers and convergence tests. These can be
31
combined by the user to achieve the desired solution, with full access
32
to the intermediate steps of the algorithms. Each class of methods uses
33
the same framework, so that you can switch between minimizers at runtime
34
without needing to recompile your program. Each instance of a minimizer
35
keeps track of its own state, allowing the minimizers to be used in
36
multi-threaded programs.
40
The header file <TT>'gsl_min.h'</TT> contains prototypes for the
41
minimization functions and related declarations. To use the minimization
42
algorithms to find the maximum of a function simply invert its sign.
48
<H2><A NAME="SEC441" HREF="gsl-ref_toc.html#TOC441">Overview</A></H2>
50
<A NAME="IDX2118"></A>
54
The minimization algorithms begin with a bounded region known to contain
55
a minimum. The region is described by a lower bound a and an
56
upper bound b, with an estimate of the location of the minimum
62
The value of the function at x must be less than the value of the
63
function at the ends of the interval,
68
f(a) > f(x) < f(b)
72
This condition guarantees that a minimum is contained somewhere within
73
the interval. On each iteration a new point x' is selected using
74
one of the available algorithms. If the new point is a better estimate
75
of the minimum, f(x') < f(x), then the current estimate of the
76
minimum x is updated. The new point also allows the size of the
77
bounded interval to be reduced, by choosing the most compact set of
78
points which satisfies the constraint f(a) > f(x) < f(b). The
79
interval is reduced until it encloses the true minimum to a desired
80
tolerance. This provides a best estimate of the location of the minimum
81
and a rigorous error estimate.
85
Several bracketing algorithms are available within a single framework.
86
The user provides a high-level driver for the algorithm, and the
87
library provides the individual functions necessary for each of the
88
steps. There are three main phases of the iteration. The steps are,
95
initialize minimizer state, <VAR>s</VAR>, for algorithm <VAR>T</VAR>
99
update <VAR>s</VAR> using the iteration <VAR>T</VAR>
103
test <VAR>s</VAR> for convergence, and repeat iteration if necessary
107
The state for the minimizers is held in a <CODE>gsl_min_fminimizer</CODE>
108
struct. The updating procedure uses only function evaluations (not
114
<H2><A NAME="SEC442" HREF="gsl-ref_toc.html#TOC442">Caveats</A></H2>
116
<A NAME="IDX2119"></A>
120
Note that minimization functions can only search for one minimum at a
121
time. When there are several minima in the search area, the first
122
minimum to be found will be returned; however it is difficult to predict
123
which of the minima this will be. <EM>In most cases, no error will be
124
reported if you try to find a minimum in an area where there is more
129
With all minimization algorithms it can be difficult to determine the
130
location of the minimum to full numerical precision. The behavior of the
131
function in the region of the minimum x^* can be approximated by
136
<PRE class="example">
137
y = f(x^*) + (1/2) f"(x^*) (x - x^*)^2
141
and the second term of this expansion can be lost when added to the
142
first term at finite precision. This magnifies the error in locating
143
x^*, making it proportional to \sqrt \epsilon (where
144
\epsilon is the relative accuracy of the floating point numbers).
145
For functions with higher order minima, such as x^4, the
146
magnification of the error is correspondingly worse. The best that can
147
be achieved is to converge to the limit of numerical accuracy in the
148
function values, rather than the location of the minimum itself.
153
<H2><A NAME="SEC443" HREF="gsl-ref_toc.html#TOC443">Initializing the Minimizer</A></H2>
157
<DT><U>Function:</U> gsl_min_fminimizer * <B>gsl_min_fminimizer_alloc</B> <I>(const gsl_min_fminimizer_type * <VAR>T</VAR>)</I>
158
<DD><A NAME="IDX2120"></A>
159
This function returns a pointer to a newly allocated instance of a
160
minimizer of type <VAR>T</VAR>. For example, the following code
161
creates an instance of a golden section minimizer,
165
<PRE class="example">
166
const gsl_min_fminimizer_type * T
167
= gsl_min_fminimizer_goldensection;
168
gsl_min_fminimizer * s
169
= gsl_min_fminimizer_alloc (T);
173
If there is insufficient memory to create the minimizer then the function
174
returns a null pointer and the error handler is invoked with an error
175
code of <CODE>GSL_ENOMEM</CODE>.
181
<DT><U>Function:</U> int <B>gsl_min_fminimizer_set</B> <I>(gsl_min_fminimizer * <VAR>s</VAR>, gsl_function * <VAR>f</VAR>, double <VAR>x_minimum</VAR>, double <VAR>x_lower</VAR>, double <VAR>x_upper</VAR>)</I>
182
<DD><A NAME="IDX2121"></A>
183
This function sets, or resets, an existing minimizer <VAR>s</VAR> to use the
184
function <VAR>f</VAR> and the initial search interval [<VAR>x_lower</VAR>,
185
<VAR>x_upper</VAR>], with a guess for the location of the minimum
186
<VAR>x_minimum</VAR>.
190
If the interval given does not contain a minimum, then the function
191
returns an error code of <CODE>GSL_FAILURE</CODE>.
197
<DT><U>Function:</U> int <B>gsl_min_fminimizer_set_with_values</B> <I>(gsl_min_fminimizer * <VAR>s</VAR>, gsl_function * <VAR>f</VAR>, double <VAR>x_minimum</VAR>, double <VAR>f_minimum</VAR>, double <VAR>x_lower</VAR>, double <VAR>f_lower</VAR>, double <VAR>x_upper</VAR>, double <VAR>f_upper</VAR>)</I>
198
<DD><A NAME="IDX2122"></A>
199
This function is equivalent to <CODE>gsl_min_fminimizer_set</CODE> but uses
200
the values <VAR>f_minimum</VAR>, <VAR>f_lower</VAR> and <VAR>f_upper</VAR> instead of
201
computing <CODE>f(x_minimum)</CODE>, <CODE>f(x_lower)</CODE> and <CODE>f(x_upper)</CODE>.
208
<DT><U>Function:</U> void <B>gsl_min_fminimizer_free</B> <I>(gsl_min_fminimizer * <VAR>s</VAR>)</I>
209
<DD><A NAME="IDX2123"></A>
210
This function frees all the memory associated with the minimizer
217
<DT><U>Function:</U> const char * <B>gsl_min_fminimizer_name</B> <I>(const gsl_min_fminimizer * <VAR>s</VAR>)</I>
218
<DD><A NAME="IDX2124"></A>
219
This function returns a pointer to the name of the minimizer. For example,
223
<PRE class="example">
224
printf ("s is a '%s' minimizer\n",
225
gsl_min_fminimizer_name (s));
229
would print something like <CODE>s is a 'brent' minimizer</CODE>.
235
<H2><A NAME="SEC444" HREF="gsl-ref_toc.html#TOC444">Providing the function to minimize</A></H2>
237
<A NAME="IDX2125"></A>
241
You must provide a continuous function of one variable for the
242
minimizers to operate on. In order to allow for general parameters the
243
functions are defined by a <CODE>gsl_function</CODE> data type
244
(see section <A HREF="gsl-ref_32.html#SEC432">Providing the function to solve</A>).
249
<H2><A NAME="SEC445" HREF="gsl-ref_toc.html#TOC445">Iteration</A></H2>
252
The following functions drive the iteration of each algorithm. Each
253
function performs one iteration to update the state of any minimizer of the
254
corresponding type. The same functions work for all minimizers so that
255
different methods can be substituted at runtime without modifications to
261
<DT><U>Function:</U> int <B>gsl_min_fminimizer_iterate</B> <I>(gsl_min_fminimizer * <VAR>s</VAR>)</I>
262
<DD><A NAME="IDX2126"></A>
263
This function performs a single iteration of the minimizer <VAR>s</VAR>. If the
264
iteration encounters an unexpected problem then an error code will be
270
<DT><CODE>GSL_EBADFUNC</CODE>
272
the iteration encountered a singular point where the function evaluated
273
to <CODE>Inf</CODE> or <CODE>NaN</CODE>.
275
<DT><CODE>GSL_FAILURE</CODE>
277
the algorithm could not improve the current best approximation or
283
The minimizer maintains a current best estimate of the position of the
284
minimum at all times, and the current interval bounding the minimum.
285
This information can be accessed with the following auxiliary functions,
290
<DT><U>Function:</U> double <B>gsl_min_fminimizer_x_minimum</B> <I>(const gsl_min_fminimizer * <VAR>s</VAR>)</I>
291
<DD><A NAME="IDX2127"></A>
292
This function returns the current estimate of the position of the
293
minimum for the minimizer <VAR>s</VAR>.
299
<DT><U>Function:</U> double <B>gsl_min_fminimizer_x_upper</B> <I>(const gsl_min_fminimizer * <VAR>s</VAR>)</I>
300
<DD><A NAME="IDX2128"></A>
301
<DT><U>Function:</U> double <B>gsl_min_fminimizer_x_lower</B> <I>(const gsl_min_fminimizer * <VAR>s</VAR>)</I>
302
<DD><A NAME="IDX2129"></A>
303
These functions return the current upper and lower bound of the interval
304
for the minimizer <VAR>s</VAR>.
310
<DT><U>Function:</U> double <B>gsl_min_fminimizer_f_minimum</B> <I>(const gsl_min_fminimizer *<VAR>s</VAR>)</I>
311
<DD><A NAME="IDX2130"></A>
312
<DT><U>Function:</U> double <B>gsl_min_fminimizer_f_upper</B> <I>(const gsl_min_fminimizer *<VAR>s</VAR>)</I>
313
<DD><A NAME="IDX2131"></A>
314
<DT><U>Function:</U> double <B>gsl_min_fminimizer_f_lower</B> <I>(const gsl_min_fminimizer *<VAR>s</VAR>)</I>
315
<DD><A NAME="IDX2132"></A>
316
These functions return the value of the function at the current estimate
317
of the minimum and at the upper and lower bounds of interval for the
318
minimizer <VAR>s</VAR>.
324
<H2><A NAME="SEC446" HREF="gsl-ref_toc.html#TOC446">Stopping Parameters</A></H2>
326
<A NAME="IDX2133"></A>
330
A minimization procedure should stop when one of the following
338
A minimum has been found to within the user-specified precision.
342
A user-specified maximum number of iterations has been reached.
346
An error has occurred.
350
The handling of these conditions is under user control. The function
351
below allows the user to test the precision of the current result.
356
<DT><U>Function:</U> int <B>gsl_min_test_interval</B> <I>(double <VAR>x_lower</VAR>, double <VAR>x_upper</VAR>, double <VAR>epsabs</VAR>, double <VAR>epsrel</VAR>)</I>
357
<DD><A NAME="IDX2134"></A>
358
This function tests for the convergence of the interval [<VAR>x_lower</VAR>,
359
<VAR>x_upper</VAR>] with absolute error <VAR>epsabs</VAR> and relative error
360
<VAR>epsrel</VAR>. The test returns <CODE>GSL_SUCCESS</CODE> if the following
361
condition is achieved,
365
<PRE class="example">
366
|a - b| < epsabs + epsrel min(|a|,|b|)
370
when the interval x = [a,b] does not include the origin. If the
371
interval includes the origin then \min(|a|,|b|) is replaced by
372
zero (which is the minimum value of |x| over the interval). This
373
ensures that the relative error is accurately estimated for minima close
378
This condition on the interval also implies that any estimate of the
379
minimum x_m in the interval satisfies the same condition with respect
380
to the true minimum x_m^*,
384
<PRE class="example">
385
|x_m - x_m^*| < epsabs + epsrel x_m^*
389
assuming that the true minimum x_m^* is contained within the interval.
396
<H2><A NAME="SEC447" HREF="gsl-ref_toc.html#TOC447">Minimization Algorithms</A></H2>
399
The minimization algorithms described in this section require an initial
400
interval which is guaranteed to contain a minimum -- if a and
401
b are the endpoints of the interval and x is an estimate
402
of the minimum then f(a) > f(x) < f(b). This ensures that the
403
function has at least one minimum somewhere in the interval. If a valid
404
initial interval is used then these algorithm cannot fail, provided the
405
function is well-behaved.
410
<DT><U>Minimizer:</U> <B>gsl_min_fminimizer_goldensection</B>
411
<DD><A NAME="IDX2135"></A>
415
<A NAME="IDX2136"></A>
416
<A NAME="IDX2137"></A>
420
The <I>golden section algorithm</I> is the simplest method of bracketing
421
the minimum of a function. It is the slowest algorithm provided by the
422
library, with linear convergence.
426
On each iteration, the algorithm first compares the subintervals from
427
the endpoints to the current minimum. The larger subinterval is divided
428
in a golden section (using the famous ratio (3-\sqrt 5)/2 =
429
0.3189660...) and the value of the function at this new point is
430
calculated. The new value is used with the constraint f(a') >
431
f(x') < f(b') to a select new interval containing the minimum, by
432
discarding the least useful point. This procedure can be continued
433
indefinitely until the interval is sufficiently small. Choosing the
434
golden section as the bisection ratio can be shown to provide the
435
fastest convergence for this type of algorithm.
442
<DT><U>Minimizer:</U> <B>gsl_min_fminimizer_brent</B>
443
<DD><A NAME="IDX2138"></A>
444
<A NAME="IDX2139"></A>
445
<A NAME="IDX2140"></A>
449
The <I>Brent minimization algorithm</I> combines a parabolic
450
interpolation with the golden section algorithm. This produces a fast
451
algorithm which is still robust.
455
The outline of the algorithm can be summarized as follows: on each
456
iteration Brent's method approximates the function using an
457
interpolating parabola through three existing points. The minimum of the
458
parabola is taken as a guess for the minimum. If it lies within the
459
bounds of the current interval then the interpolating point is accepted,
460
and used to generate a smaller interval. If the interpolating point is
461
not accepted then the algorithm falls back to an ordinary golden section
462
step. The full details of Brent's method include some additional checks
463
to improve convergence.
470
<H2><A NAME="SEC448" HREF="gsl-ref_toc.html#TOC448">Examples</A></H2>
473
The following program uses the Brent algorithm to find the minimum of
474
the function f(x) = \cos(x) + 1, which occurs at x = \pi.
475
The starting interval is (0,6), with an initial guess for the
480
<PRE class="example">
481
#include <stdio.h>
482
#include <gsl/gsl_errno.h>
483
#include <gsl/gsl_math.h>
484
#include <gsl/gsl_min.h>
486
double fn1 (double x, void * params)
495
int iter = 0, max_iter = 100;
496
const gsl_min_fminimizer_type *T;
497
gsl_min_fminimizer *s;
498
double m = 2.0, m_expected = M_PI;
499
double a = 0.0, b = 6.0;
502
F.function = &fn1;
505
T = gsl_min_fminimizer_brent;
506
s = gsl_min_fminimizer_alloc (T);
507
gsl_min_fminimizer_set (s, &F, m, a, b);
509
printf ("using %s method\n",
510
gsl_min_fminimizer_name (s));
512
printf ("%5s [%9s, %9s] %9s %10s %9s\n",
513
"iter", "lower", "upper", "min",
516
printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
518
m, m - m_expected, b - a);
523
status = gsl_min_fminimizer_iterate (s);
525
m = gsl_min_fminimizer_x_minimum (s);
526
a = gsl_min_fminimizer_x_lower (s);
527
b = gsl_min_fminimizer_x_upper (s);
530
= gsl_min_test_interval (a, b, 0.001, 0.0);
532
if (status == GSL_SUCCESS)
533
printf ("Converged:\n");
535
printf ("%5d [%.7f, %.7f] "
536
"%.7f %.7f %+.7f %.7f\n",
538
m, m_expected, m - m_expected, b - a);
540
while (status == GSL_CONTINUE && iter < max_iter);
547
Here are the results of the minimization procedure.
551
<PRE class="smallexample">
553
0 [0.0000000, 6.0000000] 2.0000000 -1.1415927 6.0000000
554
1 [2.0000000, 6.0000000] 3.2758640 +0.1342713 4.0000000
555
2 [2.0000000, 3.2831929] 3.2758640 +0.1342713 1.2831929
556
3 [2.8689068, 3.2831929] 3.2758640 +0.1342713 0.4142862
557
4 [2.8689068, 3.2831929] 3.2758640 +0.1342713 0.4142862
558
5 [2.8689068, 3.2758640] 3.1460585 +0.0044658 0.4069572
559
6 [3.1346075, 3.2758640] 3.1460585 +0.0044658 0.1412565
560
7 [3.1346075, 3.1874620] 3.1460585 +0.0044658 0.0528545
561
8 [3.1346075, 3.1460585] 3.1460585 +0.0044658 0.0114510
562
9 [3.1346075, 3.1460585] 3.1424060 +0.0008133 0.0114510
563
10 [3.1346075, 3.1424060] 3.1415885 -0.0000041 0.0077985
565
11 [3.1415885, 3.1424060] 3.1415927 -0.0000000 0.0008175
570
<H2><A NAME="SEC449" HREF="gsl-ref_toc.html#TOC449">References and Further Reading</A></H2>
573
Further information on Brent's algorithm is available in the following
581
Richard Brent, <CITE>Algorithms for minimization without derivatives</CITE>,
582
Prentice-Hall (1973), republished by Dover in paperback (2002), ISBN
587
<p>Go to the <A HREF="gsl-ref_1.html">first</A>, <A HREF="gsl-ref_32.html">previous</A>, <A HREF="gsl-ref_34.html">next</A>, <A HREF="gsl-ref_50.html">last</A> section, <A HREF="gsl-ref_toc.html">table of contents</A>.
added
removed
gsl-ref_33.html
