lmmin3 - Levenberg-Marquardt least-squares minimization with error estimates
#include <lmmin.h>
void lmmin3( const int n_par, double *par, double *parerr, double *covar, const int scale_cov, 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[];
lmmin3() 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. For fitting a parametric curve to a data set, use lmcurve3(3) or lmcurve(3).
lmmin3() replaces the obsolescent lmmin2(3), which returned its two error estimates in mutually inconsistent conventions; see ERROR ESTIMATES.
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 uncertainty vector, either of length n_par, or NULL. On output, unless it is NULL, it contains the square roots of the diagonal elements of covar,
parerr[j] = sqrt( covar[j*n_par+j] )in the convention selected by scale_cov. This identity holds in both conventions.
Covariance matrix, stored as vector of length n_par*n_par, or NULL. On output, unless it is NULL, it contains
covar = (J^T J)^-1 if scale_cov == 0
covar = (J^T J)^-1 * status.s2 if scale_cov != 0where J is the Jacobian at the solution, computed by forward differences. 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.
Selects the convention in which parerr and covar are expressed, i.e. the assumption about the errors of the data relative to which the uncertainties are defined. Use 0 if the residuals computed by evaluate are already divided by the standard deviations of the data, nonzero if the scale of the data errors is to be inferred from the residuals themselves. There is no defensible default; that is why this is an argument and not a member of control. See ERROR ESTIMATES for the formulas and for guidance.
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, lmmin3() 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 returns with *userbreak set to a nonzero value, lmmin3() will interrupt the fitting and terminate.
Parameter collection for tuning the fit procedure. In most cases, the default &lm_control_double is adequate. If fvec 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. It is reported irrespective of scale_cov, so that a caller can convert from either convention to the other without recomputing anything. 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.
covar = (J^T J)^-1
parerr[j] = sqrt( covar[j*n_par+j] )Correct if the residuals returned by evaluate are already divided by the standard deviations of the data, as lmcurve3(3) divides them by dy. This is the convention of MINPACK's covar.f, from which this code descends, and of GSL's gsl_multifit_covar(). Choose it if you trust the error bars of your data; the uncertainties then do not depend on how well the model happens to fit.
covar = (J^T J)^-1 * s2
parerr[j] = sqrt( covar[j*n_par+j] )Correct if the scale of the data errors is unknown and is to be deduced from the residuals themselves, under the assumption that the model is correct; s2 is then 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. Choose it if you trust only the relative magnitudes of your error bars, or have none at all: with unweighted residuals, i.e. plain y-f as in the example below, this is the setting you want.
The two differ by a factor that is 1 only by accident: covar by s2, parerr by sqrt(s2). Since status.s2 is reported in either case, a caller can convert after the fact,
trusted from inferred: parerr[j] / sqrt(status.s2)
inferred from trusted: parerr[j] * sqrt(status.s2)so that a program that must report both need not fit twice.
A caveat on the inferred convention: it takes the model to be correct and the residuals to be a fair sample of the data errors. If the model is wrong, or the data contain outliers, s2 exceeds 1 and the uncertainties are inflated accordingly; that is a feature only insofar as the inflation is a warning. With m_dat==n_par there is no degree of freedom left, s2 is infinite, and no uncertainty can be inferred.
lmmin2(3) predates this distinction and mixes the two conventions: it returns covar in the trusted-sigma convention but parerr in the inferred-sigma one, so that parerr[j] != sqrt(covar[j*n_par+j]). No call of lmmin3() reproduces both of its outputs at once: scale_cov=1 reproduces its parerr, scale_cov=0 its covar.
lmcurve2(3) is consistent, and equals lmcurve3(3) with scale_cov=0.
Up to lmfit 8, the Jacobian was divided by fnorm internally, so that lmmin2()'s covar was scaled and its 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. Preventing a repetition of that is why the convention is now an argument that has to be written out at every call.
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 compile 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. Beware, however, that a new parameter arrives with a default value, and that this default may not be the one your application would have chosen. This is why scale_cov is an argument of lmmin3() and not a member of control: a choice that changes the numbers you publish has to be made where it can be seen.
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 arrays 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 *info )
{
/* 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(void)
{
/* 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[3*3];
/* data points */
int m_dat = 5;
double tx[5] = { -1, -1, 1, 1, 0 };
double tz[5] = { -1, 1, -1, 1, 0 };
double y[5] = { 0, 1, 1, 2, 0.4 };
data_struct data = { tx, tz, y, f };
/* auxiliary parameters */
lm_status_struct status;
lm_control_struct control = lm_control_double;
control.verbosity = 1;
/* perform the fit; the residues above are plain y-f, hence unweighted:
the scale of the data errors is unknown, and has to be inferred from
the residues themselves. That is what scale_cov=1 requests. */
printf( "Fitting:\n" );
lmmin3( n_par, par, parerr, covar, 1, 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("residual variance s2:\n %12g\n", status.s2 );
printf("covariance matrix:\n");
int j;
for ( i = 0; i < n_par; ++i ) {
printf(" ");
for ( j = 0; j < n_par; ++j )
printf("%12g ", covar[i*n_par+j]);
printf("\n");
}
/* parerr is the square root of the diagonal of covar, in either
convention; with scale_cov=0 both would be smaller by the factor
s2 resp. sqrt(s2) */
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-2026 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>.