multifit.texi

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY without ev

@cindex nonlinear least squares fitting
@cindex least squares fitting, nonlinear

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

The header file @file{gsl_multifit_nlin.h} contains prototypes for the
multidimensional nonlinear fitting functions and related declarations.

@menu
* Overview of Nonlinear Least-Squares Fitting::  
* Initializing the Nonlinear Least-Squares Solver::  
* Providing the Function to be Minimized::  
* Iteration of the Minimization Algorithm::  
* Search Stopping Parameters for Minimization Algorithms::  
* Minimization Algorithms using Derivatives::  
* Minimization Algorithms without Derivatives::  
* Computing the covariance matrix of best fit parameters::  
* Example programs for Nonlinear Least-Squares Fitting::  
* References and Further Reading for Nonlinear Least-Squares Fitting::  
@end menu

@node Overview of Nonlinear Least-Squares Fitting
@section Overview
@cindex nonlinear least squares fitting, overview

The problem of multidimensional nonlinear least-squares fitting requires
the minimization of the squared residuals of @math{n} functions,
@math{f_i}, in @math{p} parameters, @math{x_i},
@tex
\beforedisplay
$$
\Phi(x)  = {1 \over 2} || F(x) ||^2
         = {1 \over 2} \sum_{i=1}^{n} f_i (x_1, \dots, x_p)^2
$$
\afterdisplay
@end tex
@ifinfo

@example
\Phi(x) = (1/2) || F(x) ||^2
        = (1/2) \sum_@{i=1@}^@{n@} f_i(x_1, ..., x_p)^2 
@end example

@end ifinfo
@noindent
All algorithms proceed from an initial guess using the linearization,
@tex
\beforedisplay
$$
\psi(p) = || F(x+p) || \approx || F(x) + J p\, ||
$$
\afterdisplay
@end tex
@ifinfo

@example
\psi(p) = || F(x+p) || ~=~ || F(x) + J p ||
@end example

@end ifinfo
@noindent
where @math{x} is the initial point, @math{p} is the proposed step
and @math{J} is the
Jacobian matrix @c{$J_{ij} = \partial f_i / \partial x_j$}
@math{J_@{ij@} = d f_i / d x_j}.  
Additional strategies are used to enlarge the region of convergence.
These include requiring a decrease in the norm @math{||F||} on each
step or using a trust region to avoid steps which fall outside the linear 
regime.

To perform a weighted least-squares fit of a nonlinear model
@math{Y(x,t)} to data (@math{t_i}, @math{y_i}) with independent gaussian
errors @math{\sigma_i}, use function components of the following form,
@tex
\beforedisplay
$$
f_i = {(Y(x, t_i) - y_i) \over \sigma_i}
$$
\afterdisplay
@end tex
@ifinfo

@example
f_i = (Y(x, t_i) - y_i) / \sigma_i
@end example

@end ifinfo
@noindent
Note that the model parameters are denoted by @math{x} in this chapter
since the non-linear least-squares algorithms are described
geometrically (i.e. finding the minimum of a surface).  The
independent variable of any data to be fitted is denoted by @math{t}.

With the definition above the Jacobian is 
@c{$J_{ij} = (1 / \sigma_i) \partial Y_i / \partial x_j$}
@math{J_@{ij@} =(1 / \sigma_i)  d Y_i / d x_j}, where @math{Y_i = Y(x,t_i)}.  


@node Initializing the Nonlinear Least-Squares Solver
@section Initializing the Solver

@deftypefun {gsl_multifit_fsolver *} gsl_multifit_fsolver_alloc (const gsl_multifit_fsolver_type * @var{T}, size_t @var{n}, size_t @var{p})
This function returns a pointer to a newly allocated instance of a
solver of type @var{T} for @var{n} observations and @var{p} parameters.
The number of observations @var{n} must be greater than or equal to
parameters @var{p}. 

If there is insufficient memory to create the solver then the function
returns a null pointer and the error handler is invoked with an error
code of @code{GSL_ENOMEM}.
@end deftypefun

@deftypefun {gsl_multifit_fdfsolver *} gsl_multifit_fdfsolver_alloc (const gsl_multifit_fdfsolver_type * @var{T}, size_t @var{n}, size_t @var{p})
This function returns a pointer to a newly allocated instance of a
derivative solver of type @var{T} for @var{n} observations and @var{p}
parameters.  For example, the following code creates an instance of a
Levenberg-Marquardt solver for 100 data points and 3 parameters,

@example
const gsl_multifit_fdfsolver_type * T 
    = gsl_multifit_fdfsolver_lmder;
gsl_multifit_fdfsolver * s 
    = gsl_multifit_fdfsolver_alloc (T, 100, 3);
@end example

@noindent
The number of observations @var{n} must be greater than or equal to
parameters @var{p}.

If there is insufficient memory to create the solver then the function
returns a null pointer and the error handler is invoked with an error
code of @code{GSL_ENOMEM}.
@end deftypefun

@deftypefun int gsl_multifit_fsolver_set (gsl_multifit_fsolver * @var{s}, gsl_multifit_function * @var{f}, gsl_vector * @var{x})
This function initializes, or reinitializes, an existing solver @var{s}
to use the function @var{f} and the initial guess @var{x}.
@end deftypefun

@deftypefun int gsl_multifit_fdfsolver_set (gsl_multifit_fdfsolver * @var{s}, gsl_function_fdf * @var{fdf}, gsl_vector * @var{x})
This function initializes, or reinitializes, an existing solver @var{s}
to use the function and derivative @var{fdf} and the initial guess
@var{x}.
@end deftypefun

@deftypefun void gsl_multifit_fsolver_free (gsl_multifit_fsolver * @var{s})
@deftypefunx void gsl_multifit_fdfsolver_free (gsl_multifit_fdfsolver * @var{s})
These functions free all the memory associated with the solver @var{s}.
@end deftypefun

@deftypefun {const char *} gsl_multifit_fsolver_name (const gsl_multifit_fdfsolver * @var{s})
@deftypefunx {const char *} gsl_multifit_fdfsolver_name (const gsl_multifit_fdfsolver * @var{s})
These functions return a pointer to the name of the solver.  For example,

@example
printf ("s is a '%s' solver\n", 
        gsl_multifit_fdfsolver_name (s));
@end example

@noindent
would print something like @code{s is a 'lmder' solver}.
@end deftypefun

@node Providing the Function to be Minimized
@section Providing the Function to be Minimized

You must provide @math{n} functions of @math{p} variables for the
minimization algorithms to operate on.  In order to allow for
arbitrary parameters the functions are defined by the following data
types:

@deftp {Data Type} gsl_multifit_function 
This data type defines a general system of functions with arbitrary parameters.  

@table @code
@item int (* f) (const gsl_vector * @var{x}, void * @var{params}, gsl_vector * @var{f})
this function should store the vector result
@c{$f(x,\hbox{\it params})$}
@math{f(x,params)} in @var{f} for argument @var{x} and arbitrary parameters @var{params},
returning an appropriate error code if the function cannot be computed.

@item size_t n
the number of functions, i.e. the number of components of the
vector @var{f}.

@item size_t p
the number of independent variables, i.e. the number of components of
the vector @var{x}.

@item void * params
a pointer to the arbitrary parameters of the function.
@end table
@end deftp

@deftp {Data Type} gsl_multifit_function_fdf
This data type defines a general system of functions with arbitrary parameters and
the corresponding Jacobian matrix of derivatives,

@table @code
@item int (* f) (const gsl_vector * @var{x}, void * @var{params}, gsl_vector * @var{f})
this function should store the vector result
@c{$f(x,\hbox{\it params})$}
@math{f(x,params)} in @var{f} for argument @var{x} and arbitrary parameters @var{params},
returning an appropriate error code if the function cannot be computed.

@item int (* df) (const gsl_vector * @var{x}, void * @var{params}, gsl_matrix * @var{J})
this function should store the @var{n}-by-@var{p} matrix result
@c{$J_{ij} = \partial f_i(x,\hbox{\it params}) / \partial x_j$}
@math{J_ij = d f_i(x,params) / d x_j} in @var{J} for argument @var{x} 
and arbitrary parameters @var{params}, returning an appropriate error code if the
function cannot be computed.

@item int (* fdf) (const gsl_vector * @var{x}, void * @var{params}, gsl_vector * @var{f}, gsl_matrix * @var{J})
This function should set the values of the @var{f} and @var{J} as above,
for arguments @var{x} and arbitrary parameters @var{params}.  This function
provides an optimization of the separate functions for @math{f(x)} and
@math{J(x)}---it is always faster to compute the function and its
derivative at the same time.

@item size_t n
the number of functions, i.e. the number of components of the
vector @var{f}.

@item size_t p
the number of independent variables, i.e. the number of components of
the vector @var{x}.

@item void * params
a pointer to the arbitrary parameters of the function.
@end table
@end deftp

Note that when fitting a non-linear model against experimental data,
the data is passed to the functions above using the
@var{params} argument and the trial best-fit parameters through the
@var{x} argument.

@node Iteration of the Minimization Algorithm
@section Iteration

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

@deftypefun int gsl_multifit_fsolver_iterate (gsl_multifit_fsolver * @var{s})
@deftypefunx int gsl_multifit_fdfsolver_iterate (gsl_multifit_fdfsolver * @var{s})
These functions perform a single iteration of the solver @var{s}.  If
the iteration encounters an unexpected problem then an error code will
be returned.  The solver maintains a current estimate of the best-fit
parameters at all times. 
@end deftypefun

The solver struct @var{s} contains the following entries, which can
be used to track the progress of the solution:

@table @code
@item gsl_vector * x
The current position.

@item gsl_vector * f
The function value at the current position.

@item gsl_vector * dx
The difference between the current position and the previous position,
i.e. the last step, taken as a vector.

@item gsl_matrix * J
The Jacobian matrix at the current position (for the
@code{gsl_multifit_fdfsolver} struct only)
@end table

The best-fit information also can be accessed with the following
auxiliary functions,

@deftypefun {gsl_vector *} gsl_multifit_fsolver_position (const gsl_multifit_fsolver * @var{s})
@deftypefunx {gsl_vector *} gsl_multifit_fdfsolver_position (const gsl_multifit_fdfsolver * @var{s})
These functions return the current position (i.e. best-fit parameters)
@code{s->x} of the solver @var{s}.
@end deftypefun

@node Search Stopping Parameters for Minimization Algorithms
@section Search Stopping Parameters
@cindex nonlinear fitting, stopping parameters

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

@itemize @bullet
@item
A minimum has been found to within the user-specified precision.

@item
A user-specified maximum number of iterations has been reached.

@item
An error has occurred.
@end itemize

@noindent
The handling of these conditions is under user control.  The functions
below allow the user to test the current estimate of the best-fit
parameters in several standard ways.

@deftypefun int gsl_multifit_test_delta (const gsl_vector * @var{dx}, const gsl_vector * @var{x}, double @var{epsabs}, double @var{epsrel})

This function tests for the convergence of the sequence by comparing the
last step @var{dx} with the absolute error @var{epsabs} and relative
error @var{epsrel} to the current position @var{x}.  The test returns
@code{GSL_SUCCESS} if the following condition is achieved,
@tex
\beforedisplay
$$
|dx_i| < \hbox{\it epsabs} + \hbox{\it epsrel\/}\, |x_i|