blob: cd3fd8bab9e360ed1d6f162720ee2ae534d2f7c5 [file] [log] [blame]
 =pod =begin html =end html =head1 NAME lmmin - Levenberg-Marquardt least-squares minimization =head1 SYNOPSIS B<#include > B IB<, double *>IB<, const int> IB<, constS< >void *>IB<, void *>IB<( constS< >double *>IB<, const int >IB<, constS< >void *>IB<, double *>IB<, int *>IB<), constS< >lm_control_struct *>IB<, lm_status_struct *>IB< );> B B B B =head1 DESCRIPTION B determines a vector I that minimizes the sum of squared elements of a vector I that is computed by a user-supplied function I(). On success, I represents a local minimum, not necessarily a global one; it may depend on its starting value. For applications in curve fitting, the wrapper function B 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. Function arguments: =over =item I Number of free variables. Length of parameter vector I. =item I Parameter vector. On input, it must contain a reasonable guess. On output, it contains the solution found to minimize ||I||. =item I Length of vector I. Must statisfy I <= I. =item I This pointer is ignored by the fit algorithm, except for appearing as an argument in all calls to the user-supplied routine I. =item I Pointer to a user-supplied function that computes I elements of vector I for a given parameter vector I. If I return with *I set to a negative value, B will interrupt the fitting and terminate. =item I Parameter collection for tuning the fit procedure. In most cases, the default &I is adequate. If I is only computed with single-precision accuracy, I<&lm_control_float> should be used. See also below, NOTES on initializing parameter records. I has the following members (for more details, see the source file I): =over =item B I Relative error desired in the sum of squares. Recommended setting: somewhat above machine precision; less if I is computed with reduced accuracy. =item B I Relative error between last two approximations. Recommended setting: as I. =item B I A measure for degeneracy. Recommended setting: as I. =item B I Step used to calculate the Jacobian. Recommended setting: as I, but definitely less than the accuracy of I. =item B I Initial bound to steps in the outer loop, generally between 0.01 and 100; recommended value is 100. =item B I Used to set the maximum number of function evaluations to patience*n_par. =item B I Logical switch (0 or 1). If 1, then scale parameters to their initial value. This is the recommended setting. =item B I Progress messages will be written to this file. Typically I or I. The value I will be interpreted as I. =item B I If nonzero, some progress information from within the LM algorithm is written to control.stream. =item B I -1, or maximum number of parameters to print. =item B I -1, or maximum number of residuals to print. =back =item I A record used to return information about the minimization process: =over =item B I Norm of the vector I; =item B I Actual number of iterations; =item B I Status of minimization; for the corresponding text message, print IB<[>IB<]>; for a short code, print IB<[>IB<]>. =item B I Set when termination has been forced by the user-supplied routine I. =back =back =head1 NOTES =head2 Initializing parameter records. The parameter record I 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 I. Conversely, addition of parameters is not considered an API change; it may happen without increment of the major version number. =head1 EXAMPLES =head2 Fitting a surface Fit a data set y(t) by a function f(t;p) where t is a two-dimensional vector: #include "lmmin.h" #include /* fit model: a plane p0 + p1*tx + p2*tz */ double f( double tx, double tz, const double *p ) { return p + p*tx + p*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 = { -1, 0, 1 }; /* arbitrary starting value */ /* data points */ int m_dat = 4; double tx = { -1, -1, 1, 1 }; double tz = { -1, 1, -1, 1 }; double y = { 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; ilmcurve(3) =end html =begin man \fBlmcurve\fR(3) .PP =end man Homepage: http://apps.jcns.fz-juelich.de/lmfit =head1 BUGS Please send bug reports and suggestions to the author .