lmmin2 - Levenberg-Marquardt least-squares minimization with error estimation (obsolescent)
lmmin2() is obsolescent, and may be withdrawn in lmfit 13, because the semantics of its error estimates are unclear: it returns covar and parerr in two different conventions, so that parerr is not the square root of the diagonal of covar, and the pair as a whole is correct under no single assumption about the errors of the data. For the full account, see ERROR ESTIMATES below.
The recommended replacement is lmmin3(3), which takes one additional argument, scale_cov, and returns a consistent pair either way. No migration reproduces both outputs of lmmin2() at once: scale_cov=1 reproduces this parerr, scale_cov=0 this covar.
lmmin2() itself is unchanged, in signature and in results, and remains available throughout lmfit 11 and 12.
#include <lmmin.h>
void lmmin2( const int n_par, double *par, double *parerr, double *covar, const int m_dat, const void *y, const void *data, void *evaluate( const double *par, const int m_dat, const void *data, double *fvec, int *userbreak), const lm_control_struct *control, lm_status_struct *status );
extern const lm_control_struct lm_control_double;
extern const lm_control_struct lm_control_float;
extern const char *lm_infmsg[];
extern const char *lm_shortmsg[];
lmmin2() determines a vector par that minimizes the sum of squared elements of fvec-y. The vector fvec is computed by a user-supplied function evaluate(); the vector y contains user-provided values. On success, par represents a local minimum, not necessarily a global one; it may depend on its starting value.
The function lmmin(3) provides a simplified API without error estimates; it is not obsolescent and will be maintained.
The Levenberg-Marquardt minimization starts with a steepest-descent exploration of the parameter space, and achieves rapid convergence by crossing over into the Newton-Gauss method.
Function arguments:
Number of free variables. Length of parameter vector par.
Parameter vector, of length n_par. On input, it must contain a reasonable guess. On output, it contains the solution found to minimize ||fvec||.
Parameter error vector, either of length n_par, or NULL. On output, unless it is NULL, it contains
parerr[j] = sqrt( covar[j*n_par+j] * status.s2 )i.e. the uncertainties under the assumption that the scale of the data errors is not known and is to be inferred from the residuals. Note that this is not the square root of the diagonal of covar as returned by the same call; see ERROR ESTIMATES.
Covariance matrix, stored as vector of length n_par*n_par, or NULL. On output, unless it is NULL, it contains the inverse curvature matrix
covar = (J^T J)^-1where J is the Jacobian at the solution, computed by forward differences. This is the covariance matrix under the assumption that the residuals returned by evaluate are already divided by the standard deviations of the data -- the opposite assumption to the one built into parerr; see ERROR ESTIMATES.
Both parerr and covar are set to zero if no error estimate can be obtained, which is the case if the sum of squares underflows, if evaluate fails during the finite differencing, or if the curvature matrix is singular.
Length of vector fvec. Must statisfy n_par <= m_dat.
Input vector of length m_dat. May also be the null pointer; in this case, lmmin() minimizes the squared sum of fvec instead of fvec-y.
This pointer is ignored by the fit algorithm, except for appearing as an argument in all calls to the user-supplied routine evaluate.
Pointer to a user-supplied function that computes m_dat elements of vector fvec for a given parameter vector par. If evaluate return with *userbreak set to a negative value, lmmin() will interrupt the fitting and terminate.
Parameter collection for tuning the fit procedure. In most cases, the default &lm_control_double is adequate. If f is only computed with single-precision accuracy, &lm_control_float should be used. See also below, NOTES on initializing parameter records.
control has the following members (for more details, see the source file lmstruct.h):
Relative error desired in the sum of squares. Recommended setting: somewhat above machine precision; less if fvec is computed with reduced accuracy.
Relative error between last two approximations. Recommended setting: as ftol.
A measure for degeneracy. Recommended setting: as ftol.
Step used to calculate the Jacobian. Recommended setting: as ftol, but definitely less than the accuracy of fvec.
Initial bound to steps in the outer loop, generally between 0.01 and 100; recommended value is 100.
Used to set the maximum number of function evaluations to patience*n_par.
Logical switch (0 or 1). If 1, then scale parameters to their initial value. This is the recommended setting.
Progress messages will be written to this file. Typically stdout or stderr. The value NULL will be interpreted as stdout.
If nonzero, some progress information from within the LM algorithm is written to control.stream.
-1, or maximum number of parameters to print.
-1, or maximum number of residuals to print.
A record used to return information about the minimization process:
Norm of the vector fvec;
Actual number of iterations;
Status of minimization; for the corresponding text message, print lm_infmsg[status.outcome]; for a short code, print lm_shortmsg[status.outcome].
Set when termination has been forced by the user-supplied routine evaluate.
Residual variance status.fnorm^2/(m_dat-n_par), the factor by which the two conventions of ERROR ESTIMATES differ, and the factor by which parerr and covar of this function are inconsistent with each other. Infinite if m_dat==n_par, where no degree of freedom is left to estimate it from.
Parameter uncertainties are not defined absolutely. They are defined relative to an assumption about the errors of the data, and there are two such assumptions in common use. Write J for the Jacobian at the solution and s2 = status.s2 = fnorm^2/(m_dat-n_par) for the residual variance, m_dat-n_par being the number of degrees of freedom.
The residuals returned by evaluate are already divided by the standard deviations of the data, as lmcurve2(3) and lmcurve3(3) divide them by dy. Then the covariance matrix is (J^T J)^-1, and the uncertainty of parameter j is sqrt of its j-th diagonal element. This is the convention of MINPACK's covar.f, from which this code descends, and of GSL's gsl_multifit_covar(). Appropriate if you trust the error bars of your data.
The scale of the data errors is unknown, and is deduced from the residuals themselves under the assumption that the model is correct. Then the covariance matrix is (J^T J)^-1 * s2, and the uncertainty of parameter j is sqrt of its j-th diagonal element; s2 is the reduced chi-squared if the residuals are weighted. This is what scipy.optimize.curve_fit does by default (absolute_sigma=False), and what gnuplot's fit and the summaries of R's nls report. Appropriate if you trust only the relative magnitudes of your error bars.
lmmin2() predates this distinction and mixes the two: covar follows the first convention, parerr the second. Consequently
parerr[j] != sqrt( covar[j*n_par+j] )for the same call, the two outputs differing by the factor sqrt(s2). Three ways in which this misleads have been observed:
A user reads parerr while working with weighted residuals, believing it to be the uncertainty of the weighted model. It is that uncertainty multiplied by sqrt(s2), i.e. by the square root of the reduced chi-squared, which for a good fit is near 1 and hence easy to overlook, and for a poor fit is not.
A user reads covar while working with unweighted residuals, believing it to be a covariance matrix in the units of the parameters. It is in units of residual^2 and must be multiplied by s2.
A user compares lmmin2() with lmcurve2(3) on the same fit and finds that the two disagree in parerr by sqrt(s2), because lmcurve2() applies the first convention to both of its outputs. That is what led to this documentation being written.
Up to lmfit 8, the Jacobian was divided by fnorm internally, so that covar was scaled and parerr was sqrt of its diagonal. Commit 81fcee5, released in 9.0 as "Correct parameter error estimates", corrected parerr and, as a side effect, changed the meaning of covar from scaled to raw; lmcurve2() was not adapted. Programs that moved from lmfit 8 to 9 and read covar therefore changed their reported uncertainties by a factor fnorm without any warning.
In lmmin3(3), the choice between the two conventions is an explicit argument, and parerr is the square root of the diagonal of covar in both of them.
The parameter record control should always be initialized from supplied default records:
lm_control_struct control = lm_control_double; /* or _float */After this, parameters may be overwritten:
control.patience = 500; /* allow more iterations */
control.verbosity = 15; /* for verbose monitoring */An application written this way is guaranteed to work even if new parameters are added to lm_control_struct.
Conversely, addition of parameters is not considered an API change; it may happen without increment of the major version number.
Fit a data set y(t) by a function f(t;p) where t is a two-dimensional vector:
#include "lmmin.h"
#include <stdio.h>
/* fit model: a plane p0 + p1*tx + p2*tz */
double f( double tx, double tz, const double *p )
{
return p[0] + p[1]*tx + p[2]*tz;
}
/* data structure to transmit data arays and fit model */
typedef struct {
double *tx, *tz;
double *y;
double (*f)( double tx, double tz, const double *p );
} data_struct;
/* function evaluation, determination of residues */
void evaluate_surface( const double *par, int m_dat,
const void *data, double *fvec, int *userbreak )
{
/* for readability, explicit type conversion */
data_struct *D;
D = (data_struct*)data;
int i;
for ( i = 0; i < m_dat; i++ )
fvec[i] = D->y[i] - D->f( D->tx[i], D->tz[i], par );
}
int main()
{
/* parameter vector */
int n_par = 3; /* number of parameters in model function f */
double par[3] = { -1, 0, 1 }; /* arbitrary starting value */
double parerr[3];
double covar[9];
/* data points */
int m_dat = 4;
double tx[4] = { -1, -1, 1, 1 };
double tz[4] = { -1, 1, -1, 1 };
double y[4] = { 0, 1, 1, 2 };
data_struct data = { tx, tz, y, f };
/* auxiliary parameters */
lm_status_struct status;
lm_control_struct control = lm_control_double;
control.verbosity = 3;
/* perform the fit */
printf( "Fitting:\n" );
lmmin2( n_par, par, parerr, covar, m_dat, NULL,
(const void*) &data, evaluate_surface, &control, &status );
/* print results */
printf( "\nResults:\n" );
printf( "status after %d function evaluations:\n %s\n",
status.nfev, lm_infmsg[status.outcome] );
printf("obtained parameters:\n");
int i;
for ( i=0; i<n_par; ++i )
printf(" par[%i] = %12g +- %12g\n", i, par[i], parerr[i]);
printf("obtained norm:\n %12g\n", status.fnorm );
printf("fitting data as follows:\n");
double ff;
for ( i=0; i<m_dat; ++i ){
ff = f(tx[i], tz[i], par);
printf( " t[%2d]=%12g,%12g y=%12g fit=%12g residue=%12g\n",
i, tx[i], tz[i], y[i], ff, y[i] - ff );
}
return 0;
}For more examples, see the homepage and directories demo/ and test/ in the source distribution.
Copyright (C): 1980-1999 University of Chicago 2004-2018 Joachim Wuttke, Forschungszentrum Juelich GmbH
Software: FreeBSD License
Documentation: Creative Commons Attribution Share Alike
Homepage: https://jugit.fz-juelich.de/mlz/lmfit
Please send bug reports and suggestions to the author <j.wuttke@fz-juelich.de>.