Go to the first, previous, next, last section, table of contents.
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,
\Phi(x) = (1/2) \sum_{i=1}^{n} f_i(x_1, ..., x_p)^2
= (1/2) || F(x) ||^2
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.
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.
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.
printf ("s is a '%s' solver\n",
gsl_multifit_fdfsolver_name (s));
would print something like s is a 'lmder' solver.
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:
int (* f) (const gsl_vector * x, void * params, gsl_vector * f)
size_t n
size_t p
void * params
int (* f) (const gsl_vector * x, void * params, gsl_vector * f)
int (* df) (const gsl_vector * x, void * params, gsl_matrix * J)
int (* fdf) (const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix * J)
size_t n
size_t p
void * params
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.
A minimization procedure should stop when one of the following conditions is true:
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.
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 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.
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
GSL_ETOLX
GSL_ETOLG
These error codes indicate that further iterations will be unlikely to change the solution from its current value.
There are no algorithms implemented in this section at the moment.
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));
{
double chi = gsl_blas_dnrm2(s->f);
printf("chisq/dof = %g\n", pow(chi, 2.0)/ (n - p));
}
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 chisq/dof = 1.08335 status = success
The approximate values of the parameters are found correctly, and the chi-squared value indicates a good fit (the chi-squared per degree of freedom is approximately 1). In this case the errors on the parameters can be estimated from the square roots of the diagonal elements of the covariance matrix. If the chi-squared value indicates a poor fit then error estimates obtained from the covariance matrix are not valid, since the Gaussian approximation would not apply.
The MINPACK algorithm is described in the following article,
The following paper is also relevant to the algorithms described in this section,
Go to the first, previous, next, last section, table of contents.