lmmin - Levenberg-Marquardt least-squares minimization
#include <lmmin.h>
void lmmin( int n_par, double *par, int m_dat, const void *data, void *evaluate( const double *par, 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[];
lmmin() determines a vector par that minimizes the sum of squared elements of a vector fvec that is computed by a user-supplied function evaluate(). On success, par represents a local minimum, not necessarily a global one; it may depend on its starting value.
For applications in curve fitting, the wrapper function lmcurve(3) offers a simplified API.
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. For more background, see Moré 1978, Madsen et al. 2004, and comments in the source code.
Function arguments:
Number of free variables. Dimension of parameter vector par.
Parameter vector. On input, it must contain a reasonable guess; on output, it contains the solution found to minimize ||fvec||.
Dimension of vector fvec. Must statisfy n_par <= m_dat.
This pointer is ignored by the fit algorithm, except for appearing as an argument in all calls to the user-supplied routines evaluate and printout.
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.
A function that prints messages about the fit progress. Legitimate values of this pointer argument are NULL (to indicate that no messages are wanted), &lm_printout_std (pointer to a default implementation), or a pointer to a user-supplied function.
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.
Pointer to the stream the output shall be writen to. Typically &stdout or &stderr.
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 forded by the user supplied routine evaluate.
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 */
/* 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" ); lmmin( n_par, par, m_dat, (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\n", i, par[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-2013 Joachim Wuttke, Forschungszentrum Juelich GmbH
Software: FreeBSD License
Documentation: Creative Commons Attribution Share Alike
Homepage: http://apps.jcns.fz-juelich.de/lmfit
Please report bugs to the author <j.wuttke@fz-juelich.de>