lmmin2 - Man Page

Levenberg-Marquardt least-squares minimization with error estimation

Synopsis

#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[];

Description

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 legacy function lmmin(3) provides a simplified API without 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:

n_par

Number of free variables. Length of parameter vector par.

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||.

parerr

Parameter error vector, either of length n_par, or NULL. On output, unless it is NULL, it contains the parameter uncertainties, estimated from the diagonal elements of the covariance matrix.

covar

Covariance matrix, stored as vector of length n_par*n_par, or NULL. On output, unless it is NULL, it contains the covariance matrix. NOT YET IMPLEMENTED!

m_dat

Length of vector fvec. Must statisfy n_par <= m_dat.

y

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.

data

This pointer is ignored by the fit algorithm, except for appearing as an argument in all calls to the user-supplied routine evaluate.

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.

control

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):

double control.ftol: Relative error desired in the sum of squares. Recommended setting: somewhat above machine precision; less if fvec is computed with reduced accuracy.
double control.xtol: Relative error between last two approximations. Recommended setting: as ftol.
double control.gtol: A measure for degeneracy. Recommended setting: as ftol.
double control.epsilon: Step used to calculate the Jacobian. Recommended setting: as ftol, but definitely less than the accuracy of fvec.
double control.stepbound: Initial bound to steps in the outer loop, generally between 0.01 and 100; recommended value is 100.
int control.patience: Used to set the maximum number of function evaluations to patience*n_par.
int control.scale_diag: Logical switch (0 or 1). If 1, then scale parameters to their initial value. This is the recommended setting.
FILE* control.msgfile: Progress messages will be written to this file. Typically stdout or stderr. The value NULL will be interpreted as stdout.
int control.verbosity: If nonzero, some progress information from within the LM algorithm is written to control.stream.
int control.n_maxpri: -1, or maximum number of parameters to print.
int control.m_maxpri: -1, or maximum number of residuals to print.

status

A record used to return information about the minimization process:

double status.fnorm: Norm of the vector fvec;
int status.nfev: Actual number of iterations;
int status.outcome: Status of minimization; for the corresponding text message, print lm_infmsg[status.outcome]; for a short code, print lm_shortmsg[status.outcome].
int status.userbreak: Set when termination has been forced by the user-supplied routine evaluate.

Notes

Initializing parameter records.

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.

Examples

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 <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, (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;
    }

lmmin2 - Man Page

Synopsis

Description

Notes

Initializing parameter records.

Examples

Fitting a surface

More examples

Copying

See Also

Bugs

Referenced By