NAME

lmmin - Levenberg-Marquardt least-squares minimization


SYNOPSIS

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


DESCRIPTION

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:

n_par

Number of free variables. Dimension of parameter vector par.

par

Parameter vector. On input, it must contain a reasonable guess; on output, it contains the solution found to minimize ||fvec||.

m_dat

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

data

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.

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.

printout

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.

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

Pointer to the stream the output shall be writen to. Typically &stdout or &stderr.

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 forded 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 */
        /* 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;
    }

More examples

For more examples, see the homepage and directories demo/ and test/ in the source distribution.


COPYING

Copyright (C): 1980-1999 University of Chicago 2004-2013 Joachim Wuttke, Forschungszentrum Juelich GmbH

Software: FreeBSD License

Documentation: Creative Commons Attribution Share Alike


SEE ALSO

lmcurve(3)

Homepage: http://apps.jcns.fz-juelich.de/lmfit

Please report bugs to the author <j.wuttke@fz-juelich.de>