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 `gsl_multifit_nlin.h' contains prototypes for the
multidimensional nonlinear fitting functions and related declarations.
The problem of multidimensional nonlinear least-squares fitting requires
the minimization of the squared residuals of n functions,
f_i, in p parameters, x_i,
All algorithms proceed from an initial guess using the linearization,
\psi(p) = || F(x+p) || ~=~ || F(x) + J p ||
where x is the initial point, p is the proposed step
and J is the
Jacobian matrix
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 ||F|| on each
step or using a trust region to avoid steps which fall outside the linear
regime.
This function returns a pointer to a a newly allocated instance of a
solver of type T for n observations and p parameters.
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 GSL_ENOMEM.
This function returns a pointer to a a newly allocated instance of a
derivative solver of type T for n observations and p
parameters. For example, the following code creates an instance of a
Levenberg-Marquardt solver for 100 data points and 3 parameters,
const gsl_multifit_fdfsolver_type * T
= gsl_multifit_fdfsolver_lmder;
gsl_multifit_fdfsolver * s
= gsl_multifit_fdfsolver_alloc (T, 100, 3);
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 GSL_ENOMEM.
You must provide n functions of p variables for the minimization algorithms to operate on. In order to allow for general parameters the
functions are defined by the following data types:
Data Type:gsl_multifit_function
This data type defines a general system of functions with parameters.
this function should store the vector result
f(x,params) in f for argument x and parameters params,
returning an appropriate error code if the function cannot be computed.
size_t n
the number of functions, i.e. the number of components of the
vector f
size_t p
the number of independent variables, i.e. the number of components of
the vectors x
void * params
a pointer to the parameters of the function
Data Type:gsl_multifit_function_fdf
This data type defines a general system of functions with parameters and
the corresponding Jacobian matrix of derivatives,
this function should store the vector result
f(x,params) in f for argument x and parameters params,
returning an appropriate error code if the function cannot be computed.
this function should store the n-by-p matrix result
J_ij = d f_i(x,params) / d x_j in J for argument x
and parameters params, returning an appropriate error code if the
function cannot be computed.
This function should set the values of the f and J as above,
for arguments x and parameters params. This function provides
an optimization of the separate functions for f(x) and J(x) --
it is always faster to compute the function and its derivative at the
same time.
size_t n
the number of functions, i.e. the number of components of the
vector f
size_t p
the number of independent variables, i.e. the number of components of
the vectors x
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.
Function: int gsl_multifit_fsolver_iterate(gsl_multifit_fsolver * s)
Function: int gsl_multifit_fdfsolver_iterate(gsl_multifit_fdfsolver * s)
These functions perform a single iteration of the solver 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. This information can be accessed with the
following auxiliary functions,
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 been reached.
An error has occurred.
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.
This function tests for the convergence of the sequence by comparing the
last step dx with the absolute error epsabs and relative
error epsrel to the current position x. The test returns
GSL_SUCCESS if the following condition is achieved,
|dx_i| < epsabs + epsrel |x_i|
for each component of x and returns GSL_CONTINUE otherwise.
Function: int gsl_multifit_test_gradient(const gsl_vector * g, double epsabs)
This function tests the residual gradient g against the absolute
error bound epsabs. Mathematically, the gradient should be
exactly zero at the minimum. The test returns GSL_SUCCESS if the
following condition is achieved,
\sum_i |g_i| < epsabs
and returns GSL_CONTINUE otherwise. This criterion is suitable
for situations where the the precise location of the minimum, x,
is unimportant provided a value can be found where the gradient is small
enough.
The minimization algorithms described in this section make use of both
the function and its derivative. They require an initial guess for the
location of the minimum. There is no absolute guarantee of convergence
-- the function must be suitable for this technique and the initial
guess must be sufficiently close to the minimum for it to work.
Derivative Solver:gsl_multifit_fdfsolver_lmsder
This is a robust and efficient version of the Levenberg-Marquardt
algorithm as implemented in the scaled LMDER routine in
MINPACK. Minpack was written by Jorge J. Mor'e, Burton S. Garbow
and Kenneth E. Hillstrom.
The algorithm uses a generalized trust region to keep each step under
control. In order to be accepted a proposed new position x' must
satisfy the condition |D (x' - x)| < \delta, where D is a
diagonal scaling matrix and \delta is the size of the trust
region. The components of D are computed internally, using the
column norms of the Jacobian to estimate the sensitivity of the residual
to each component of x. This improves the behavior of the
algorithm for badly scaled functions.
On each iteration the algorithm attempts to minimize the linear system
|F + J p| subject to the constraint |D p| < \Delta. The
solution to this constrained linear system is found using the
Levenberg-Marquardt method.
The proposed step is now tested by evaluating the function at the
resulting point, x'. If the step reduces the norm of the
function sufficiently, and follows the predicted behavior of the
function within the trust region. then it is accepted and size of the
trust region is increased. If the proposed step fails to improve the
solution, or differs significantly from the expected behavior within
the trust region, then the size of the trust region is decreased and
another trial step is computed.
The algorithm also monitors the progress of the solution and returns an
error if the changes in the solution are smaller than the machine
precision. The possible error codes are,
GSL_ETOLF
the decrease in the function falls below machine precision
GSL_ETOLX
the change in the position vector falls below machine precision
GSL_ETOLG
the norm of the gradient, relative to the norm of the function, falls
below machine precision
These error codes indicate that further iterations will be unlikely to
change the solution from its current value.
Derivative Solver:gsl_multifit_fdfsolver_lmder
This is an unscaled version of the LMDER algorithm. The elements of the
diagonal scaling matrix D are set to 1. This algorithm may be
useful in circumstances where the scaled version of LMDER converges too
slowly, or the function is already scaled appropriately.
This function uses the Jacobian matrix J to compute the covariance
matrix of the best-fit parameters, covar. The parameter
epsrel is used to remove linear-dependent columns when J is
rank deficient.
The covariance matrix is given by,
covar = (J^T J)^{-1}
and is computed by QR decomposition of J with column-pivoting. Any
columns of R which satisfy
|R_{kk}| <= epsrel |R_{11}|
are considered linearly-dependent and are excluded from the covariance
matrix (the corresponding rows and columns of the covariance matrix are
set to zero).
The following example program fits a weighted exponential model with
background to experimental data, Y = A \exp(-\lambda t) + b. The
first part of the program sets up the functions expb_f and
expb_df to calculate the model and its Jacobian. The appropriate
fitting function is given by,
f_i = ((A \exp(-\lambda t_i) + b) - y_i)/\sigma_i
where we have chosen t_i = i. The Jacobian matrix J is
the derivative of these functions with respect to the three parameters
(A, \lambda, b). It is given by,
J_{ij} = d f_i / d x_j
where x_0 = A, x_1 = \lambda and x_2 = b.
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>
struct data {
size_t n;
double * y;
double * sigma;
};
int
expb_f (const gsl_vector * x, void *params,
gsl_vector * f)
{
size_t n = ((struct data *)params)->n;
double *y = ((struct data *)params)->y;
double *sigma = ((struct data *) params)->sigma;
double A = gsl_vector_get (x, 0);
double lambda = gsl_vector_get (x, 1);
double b = gsl_vector_get (x, 2);
size_t i;
for (i = 0; i < n; i++)
{
/* Model Yi = A * exp(-lambda * i) + b */
double t = i;
double Yi = A * exp (-lambda * t) + b;
gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
}
return GSL_SUCCESS;
}
int
expb_df (const gsl_vector * x, void *params,
gsl_matrix * J)
{
size_t n = ((struct data *)params)->n;
double *sigma = ((struct data *) params)->sigma;
double A = gsl_vector_get (x, 0);
double lambda = gsl_vector_get (x, 1);
size_t i;
for (i = 0; i < n; i++)
{
/* Jacobian matrix J(i,j) = dfi / dxj, */
/* where fi = (Yi - yi)/sigma[i], */
/* Yi = A * exp(-lambda * i) + b */
/* and the xj are the parameters (A,lambda,b) */
double t = i;
double s = sigma[i];
double e = exp(-lambda * t);
gsl_matrix_set (J, i, 0, e/s);
gsl_matrix_set (J, i, 1, -t * A * e/s);
gsl_matrix_set (J, i, 2, 1/s);
}
return GSL_SUCCESS;
}
int
expb_fdf (const gsl_vector * x, void *params,
gsl_vector * f, gsl_matrix * J)
{
expb_f (x, params, f);
expb_df (x, params, J);
return GSL_SUCCESS;
}
The main part of the program sets up a Levenberg-Marquardt solver and
some simulated random data. The data uses the known parameters
(1.0,5.0,0.1) combined with gaussian noise (standard deviation = 0.1)
over a range of 40 timesteps. The initial guess for the parameters is
chosen as (0.0, 1.0, 0.0).
#define N 40
int
main (void)
{
const gsl_multifit_fdfsolver_type *T;
gsl_multifit_fdfsolver *s;
int status;
size_t i, iter = 0;
const size_t n = N;
const size_t p = 3;
gsl_matrix *covar = gsl_matrix_alloc (p, p);
double y[N], sigma[N];
struct data d = { n, y, sigma};
gsl_multifit_function_fdf f;
double x_init[3] = { 1.0, 0.0, 0.0 };
gsl_vector_view x = gsl_vector_view_array (x_init, p);
const gsl_rng_type * type;
gsl_rng * r;
gsl_rng_env_setup();
type = gsl_rng_default;
r = gsl_rng_alloc (type);
f.f = &expb_f;
f.df = &expb_df;
f.fdf = &expb_fdf;
f.n = n;
f.p = p;
f.params = &d;
/* This is the data to be fitted */
for (i = 0; i < n; i++)
{
double t = i;
y[i] = 1.0 + 5 * exp (-0.1 * t)
+ gsl_ran_gaussian(r, 0.1);
sigma[i] = 0.1;
printf("data: %d %g %g\n", i, y[i], sigma[i]);
};
T = gsl_multifit_fdfsolver_lmsder;
s = gsl_multifit_fdfsolver_alloc (T, n, p);
gsl_multifit_fdfsolver_set (s, &f, &x.vector);
print_state (iter, s);
do
{
iter++;
status = gsl_multifit_fdfsolver_iterate (s);
printf ("status = %s\n", gsl_strerror (status));
print_state (iter, s);
if (status)
break;
status = gsl_multifit_test_delta (s->dx, s->x,
1e-4, 1e-4);
}
while (status == GSL_CONTINUE && iter < 500);
gsl_multifit_covar (s->J, 0.0, covar);
gsl_matrix_fprintf (stdout, covar, "%g");
#define FIT(i) gsl_vector_get(s->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
printf("A = %.5f +/- %.5f\n", FIT(0), ERR(0));
printf("lambda = %.5f +/- %.5f\n", FIT(1), ERR(1));
printf("b = %.5f +/- %.5f\n", FIT(2), ERR(2));
printf ("status = %s\n", gsl_strerror (status));
gsl_multifit_fdfsolver_free (s);
return 0;
}
int
print_state (size_t iter, gsl_multifit_fdfsolver * s)
{
printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
"|f(x)| = %g\n",
iter,
gsl_vector_get (s->x, 0),
gsl_vector_get (s->x, 1),
gsl_vector_get (s->x, 2),
gsl_blas_dnrm2 (s->f));
}
The iteration terminates when the change in x is smaller than 0.0001, as
both an absolute and relative change. Here are the results of running
the program,
iter: 0 x = 1.00000000 0.00000000 0.00000000 |f(x)| = 118.574
iter: 1 x = 1.64919392 0.01780040 0.64919392 |f(x)| = 77.2068
iter: 2 x = 2.86269020 0.08032198 1.45913464 |f(x)| = 38.0579
iter: 3 x = 4.97908864 0.11510525 1.06649948 |f(x)| = 10.1548
iter: 4 x = 5.03295496 0.09912462 1.00939075 |f(x)| = 6.4982
iter: 5 x = 5.05811477 0.10055914 0.99819876 |f(x)| = 6.33121
iter: 6 x = 5.05827645 0.10051697 0.99756444 |f(x)| = 6.33119
iter: 7 x = 5.05828006 0.10051819 0.99757710 |f(x)| = 6.33119
A = 5.05828 +/- 0.05983
lambda = 0.10052 +/- 0.00309
b = 0.99758 +/- 0.03944
status = success
The approximate values of the parameters are found correctly. The
errors on the parameters are given by the square roots of the diagonal
elements of the covariance matrix.