This chapter describes functions for finding minima of arbitrary one-dimensional functions.  The library provides low level components for a variety of iterative minimizers and convergence tests. These can be combined by the user to achieve the desired solution, with full access to the intermediate steps of the algorithms. Each class of methods uses the same framework, so that you can switch between minimizers at runtime without needing to recompile your program. Each instance of a minimizer keeps track of its own state, allowing the minimizers to be used in multi-threaded programs.

The header file gsl_min.h contains prototypes for the minimization functions and related declarations. To use the minimization algorithms to find the maximum of a function simply invert its sign.

Overview
========

The minimization algorithms begin with a bounded region known to contain a minimum. The region is described by an lower bound a and an upper bound b, with an estimate of the minimum x.

The value of the function at x must be less than the value of the function at the ends of the interval,

     f(a) > f(x) < f(b)

This condition guarantees that a minimum is contained somewhere within the interval. On each iteration a new point x' is selected using one of the available algorithms.  If the new point is a better estimate of the minimum, f(x') < f(x), then the current estimate of the minimum x is updated.  The new point also allows the size of the bounded interval to be reduced, by choosing the most compact set of points which satisfies the constraint f(a) > f(x) < f(b).  The interval is reduced until it encloses the true minimum to a desired tolerance. This provides a best estimate of the location of the minimum and a rigorous error estimate.

In GSL different bracketing algorithms are available in a similar framework. The user provides a high-level driver for the algorithms, and the library provides the individual functions necessary for each of the steps.  There are three main phases of the iteration. The steps are,

  * initialize minimizer state, S, for algorithm T

  * update S using the iteration T

  * test S for convergence, and repeat iteration if necessary

The state for the minimizers is held in a gsl_min_fminimizer struct. The updating procedure uses only function evaluations (not derivatives).

Caveats
=======

Note that minimization functions can only search for one minimum at a time.  When there are several minima in the search area, the first minimum to be found will be returned; however it is difficult to predict which of the minima this will be. In most cases, no error will be reported if you try to find a minimum in an area where there is more than one.

With all minimization algorithms it can be difficult to determine the location of the minimum to full numerical precision. The behavior of the function in the region of the minimum x^* can be approximated by a Taylor expansion,

     y = f(x^*) + (1/2) f''(x^*) (x - x^*)^2

and the second term of this expansion can be lost when added to the first term at finite precision. This magnifies the error in locating x^*, making it proportional to \sqrt \epsilon (where \epsilon is the relative accuracy of the floating point numbers). For functions with higher order minima, such as x^4, the magnification of the error is correspondingly worse. The best that can be achieved is to converge to the limit of numerical accuracy in the function values, rather than the location of the minimum itself.

Initializing the Minimizer
==========================

 -- Function: gsl_min_fminimizer * gsl_min_fminimizer_alloc (const gsl_min_fminimizer_type * T, gsl_function * F, double MINIMUM, gsl_interval X)
     This function returns a pointer to a a newly allocated instance of a minimizer of type T for the function F, with an initial bracket on the minimum of X containing a guess for the location of the minimum MINIMUM. For example, the following code creates an instance of a golden section minimizer,

          gsl_min_fminimizer * s =
              gsl_min_fminimizer_alloc (gsl_min_fminimizer_goldensection);

     If there is insufficient memory to create the minimizer then the function returns a null pointer and the error handler is invoked with an error code of GSL_ENOMEM.

 -- Function: gsl_min_fminimizer * gsl_min_fminimizer_alloc_with_values (const gsl_min_fminimizer_type * T, gsl_function * F, double MINIMUM, double F_MINIMUM, gsl_interval X, double F_LOWER, double F_UPPER)
     This function is equivalent to gsl_min_fminimizer_alloc but uses the values F_MINIMUM, F_LOWER and F_UPPER instead of computing f(minimum), f(x.lower) and f(x.upper).

 -- Function: int gsl_min_fminimizer_set (gsl_min_fminimizer * S, gsl_function * F, double MINIMUM, gsl_interval X)
     This function reinitializes an existing minimizer S to use the function F and the initial search interval X, with a guess for the location of the minimum MINIMUM.

 -- Function: int gsl_min_fminimizer_set_with_values (gsl_min_fminimizer * S, gsl_function * F, double MINIMUM, double F_MINIMUM, gsl_interval X, double F_LOWER, double F_UPPER)
     This function is equivalent to gsl_min_fminimizer_set but uses the values F_MINIMUM, F_LOWER and F_UPPER instead of computing f(minimum), f(x.lower) and f(x.upper).

 -- Function: void gsl_min_fminimizer_free (gsl_min_fminimizer * S)
     This function frees all the memory associated with the minimizer S.

 -- Function: const char * gsl_min_fminimizer_name (const gsl_min_fminimizer * S)
     This function returns a pointer to the name of the minimizer. For example,

          printf("s is a '%s' minimizer\n", gsl_min_fminimizer_name (s)) ;

     would print something like s is a 'brent' minimizer

Providing the function to minimize
==================================

You must provide a continuous function of one variable for the minimizers to operate on. In order to allow for general parameters the functions are defined by a gsl_function data type (see Providing the function to solve).

Iteration
=========

The following functions drive the iteration of each algorithm. Each function performs one iteration to update the state of any minimizer of the corresponding type. The same functions work for all minimizers so that different methods can be substituted at runtime without modifications to the code.

 -- Function: int gsl_min_fminimizer_iterate (gsl_min_fminimizer * S)
     This function performs a single iteration of the minimizer S. If the iteration encounters an unexpected problem then an error code will be returned,

     GSL_EBADFUNC
          the iteration encountered a singular point where the function evaluated to Inf or NaN.

     GSL_FAILURE
          the algorithm could not improve the current best approximation or bounding interval.

The minimizer maintains a current best estimate of the minimum at all times, and the current interval bounding the minimum. This information can be accessed with the following auxiliary functions,

 -- Function: double gsl_min_fminimizer_minimum (const gsl_min_fminimizer * S)
     This function returns the current estimate of the minimum for the minimizer S.

 -- Function: gsl_interval gsl_min_fminimizer_interval (const gsl_min_fminimizer * S)
     This function returns the current bounding interval for the minimizer S.

Stopping Parameters
===================

A minimization procedure should stop when one of the following conditions is true:

  * A minimum has been found to within the user-specified precision.

  * A user-specified maximum number of iterations has executed.

  * An error has occurred.

In the minimization framework of GSL the handling of these conditions is under user control.  The function below allows the user to test the precision of the current result.

 -- Function: int gsl_min_test_interval (gsl_interval X, double EPSREL, double EPSABS)
     This function tests for the convergence of the interval X with absolute error EPSABS and relative error EPSREL. The test returns GSL_SUCCESS if the following condition is achieved,

          |a - b| < epsabs + epsrel min(|a|,|b|)

     when the interval x = [a,b] does not include the origin.  If the interval includes the origin then \min(|a|,|b|) is replaced by zero (which is the minimum value of |x| over the interval). This ensures that the relative error is accurately estimated for minima close to the origin.

     This condition on the interval also implies that any estimate of the minimum x_m in the interval satisfies the same condition with respect to the true minimum x_m^*,

          |x_m - x_m^*| < epsabs + epsrel x_m^*

     assuming that the true minimum x_m^* is contained within the interval.

Minimization Algorithms
=======================

The minimization algorithms described in this section require an initial interval which is guaranteed to contain a minimum --- if a and b are the endpoints of the interval and x is an estimate of the minimum then f(a) > f(x) < f(b). This ensures that the function has at least one minimum somewhere in the interval. If a valid initial interval is used then these algorithm cannot fail, provided the function is well-behaved.

 -- Minimizer: gsl_min_fminimizer_goldensection

     The golden section algorithm is the simplest method of bracketing the minimum of a function.  It is the slowest algorithm provided by the library, with linear convergence.

     On each iteration, the algorithm first compares the subintervals from the endpoints to the current minimum. The larger subinterval is divided in a golden section (using the famous ratio (3-\sqrt 5)/2) and the value of the function at this new point is calculated. The new value is used with the constraint f(a') > f(x') < f(b') to a select new interval containing the minimum, by discarding the least useful point. This procedure can be continued indefinitely until the interval is sufficiently small.  Choosing the golden section as the bisection ratio can be shown to provide the fastest convergence for this type of algorithm.

 -- Minimizer: gsl_min_fminimizer_brent

     The Brent minimization algorithm combines a parabolic interpolation with the golden section algorithm. This produces a fast algorithm which is still robust.

     The outline of the algorithm can be summarized as follows: on each iteration Brent's method approximates the function using an interpolating parabola through three existing points. The minimum of the parabola is taken as a guess for the minimum.  If it lies within the bounds of the current interval then the interpolating point is accepted, and used to generate a smaller interval. If the interpolating point is not accepted then the algorithm falls back to an ordinary golden section step. The full details of Brent's method include some additional checks to improve convergence.

Examples
========

The following program uses the Brent algorithm to find the minimum of the function f(x) = \cos(x) + 1, which occurs at x = \pi. The starting interval is (0,6), with an initial guess for the minimum of 2.

     #include <stdio.h>
     #include <gsl_errno.h>
     #include <gsl_math.h>
     #include <gsl_min.h>

     double fn1 (double x, void * params);

     double fn1 (double x, void * params)
     {
       return cos(x) + 1.0 ;
     }

     int
     main ()
     {
       int status;
       int iterations = 0, max_iterations = 100;
       gsl_min_fminimizer *s;
       double m = 2.0, m_expected = M_PI;
       gsl_interval x =  {0, 6.0};
       gsl_function F;

       F.function = &fn1;
       F.params = 0;

       s = gsl_min_fminimizer_alloc (gsl_min_fminimizer_brent, &F, m, x);

       printf ("using %s method\n", gsl_min_fminimizer_name (s));

       printf ("%5s [%9s, %9s] %9s %9s %10s %9s\n",
     	  "iter", "lower", "upper", "min", "actual", "err", "err(est)");

       printf ("%5d [%.7f, %.7f] %.7f %.7f %+.7f %.7f\n",
               iterations, x.lower, x.upper,
               m, m_expected, m - m_expected, x.upper - x.lower);

       do
         {
           iterations++;
           status = gsl_min_fminimizer_iterate (s);

           m = gsl_min_fminimizer_minimum (s);
           x = gsl_min_fminimizer_interval (s);

           status = gsl_min_test_interval (x, 0.001, 0.0);

           if (status == GSL_SUCCESS)
     	printf ("Converged:\n");

           printf ("%5d [%.7f, %.7f] %.7f %.7f %+.7f %.7f\n",
     	      iterations, x.lower, x.upper,
     	      m, m_expected, m - m_expected, x.upper - x.lower);
         }
       while (status == GSL_CONTINUE && iterations < max_iterations);

     }

Here are the results of the minimization procedure.

     bash$ ./a.out
         0 [0.0000000, 6.0000000] 2.0000000 3.1415927 -1.1415927 6.0000000
         1 [2.0000000, 6.0000000] 3.2758640 3.1415927 +0.1342713 4.0000000
         2 [2.0000000, 3.2831929] 3.2758640 3.1415927 +0.1342713 1.2831929
         3 [2.8689068, 3.2831929] 3.2758640 3.1415927 +0.1342713 0.4142862
         4 [2.8689068, 3.2831929] 3.2758640 3.1415927 +0.1342713 0.4142862
         5 [2.8689068, 3.2758640] 3.1460585 3.1415927 +0.0044658 0.4069572
         6 [3.1346075, 3.2758640] 3.1460585 3.1415927 +0.0044658 0.1412565
         7 [3.1346075, 3.1874620] 3.1460585 3.1415927 +0.0044658 0.0528545
         8 [3.1346075, 3.1460585] 3.1460585 3.1415927 +0.0044658 0.0114510
         9 [3.1346075, 3.1460585] 3.1424060 3.1415927 +0.0008133 0.0114510
        10 [3.1346075, 3.1424060] 3.1415885 3.1415927 -0.0000041 0.0077985
     Converged:
        11 [3.1415885, 3.1424060] 3.1415927 3.1415927 -0.0000000 0.0008175

References and Further Reading
==============================

Further information on Brent's algorithm is available in the following book,

  Richard Brent, Algorithms for minimization without derivatives, Prentice-Hall (1973)