man/lmmin.pod

=pod

=begin html

<link rel="stylesheet" href="podstyle.css" type="text/css" />

=end html

=head1 NAME

lmmin - Levenberg-Marquardt least-squares minimization


=head1 SYNOPSIS

B<#include <lmmin.h>>

B<void lmmin( const int> I<n_par>B<, double *>I<par>B<, const int> I<m_dat>B<,
            constS< >void *>I<data>B<,
            void *>I<evaluate>B<(
                 constS< >double *>I<par>B<, const int >I<m_dat>B<,
                 constS< >void *>I<data>B<, double *>I<fvec>B<, int *>I<userbreak>B<),
            constS< >lm_control_struct *>I<control>B<,
            lm_status_struct *>I<status>B< );>

B<extern const lm_control_struct lm_control_double;>

B<extern const lm_control_struct lm_control_float;>

B<extern const char *lm_infmsg[];>

B<extern const char *lm_shortmsg[];>

=head1 DESCRIPTION

B<lmmin()> determines a vector I<par> that minimizes the sum of squared elements of a vector I<fvec> that is computed by a user-supplied function I<evaluate>().
On success, I<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 B<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.

Function arguments:

=over

=item I<n_par>

Number of free variables.
Length of parameter vector I<par>.

=item I<par>

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

=item I<m_dat>

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

=item I<data>

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

=item I<evaluate>

Pointer to a user-supplied function that computes I<m_dat> elements of vector I<fvec> for a given parameter vector I<par>.
If I<evaluate> return with *I<userbreak> set to a negative value, B<lmmin()> will interrupt the fitting and terminate.

=item I<control>

Parameter collection for tuning the fit procedure.
In most cases, the default &I<lm_control_double> is adequate.
If I<f> is only computed with single-precision accuracy,
I<&lm_control_float> should be used.
See also below, NOTES on initializing parameter records.

I<control> has the following members (for more details, see the source file I<lmstruct.h>):

=over

=item B<double> I<control.ftol>

Relative error desired in the sum of squares.
Recommended setting: somewhat above machine precision; less if I<fvec> is computed with reduced accuracy.

=item B<double> I<control.xtol>

Relative error between last two approximations.
Recommended setting: as I<ftol>.

=item B<double> I<control.gtol>

A measure for degeneracy.
Recommended setting: as I<ftol>.

=item B<double> I<control.epsilon>

Step used to calculate the Jacobian.
Recommended setting: as I<ftol>, but definitely less than the accuracy of I<fvec>.

=item B<double> I<control.stepbound>

Initial bound to steps in the outer loop, generally between 0.01 and 100; recommended value is 100.

=item B<int> I<control.patience>

Used to set the maximum number of function evaluations to patience*n_par.

=item B<int> I<control.scale_diag>

Logical switch (0 or 1).
If 1, then scale parameters to their initial value.
This is the recommended setting.

=item B<FILE*> I<control.msgfile>

Progress messages will be written to this file.
Typically I<stdout> or I<stderr>.
The value I<NULL> will be interpreted as I<stdout>.

=item B<int> I<control.verbosity>

If nonzero, some progress information from within the LM algorithm
is written to control.stream.

=item B<int> I<control.n_maxpri>

-1, or maximum number of parameters to print.

=item B<int> I<control.m_maxpri>

-1, or maximum number of residuals to print.

=back

=item I<status>

A record used to return information about the minimization process:

=over

=item B<double> I<status.fnorm>

Norm of the vector I<fvec>;

=item B<int> I<status.nfev>

Actual number of iterations;

=item B<int> I<status.outcome>

Status of minimization;
for the corresponding text message, print I<lm_infmsg>B<[>I<status.outcome>B<]>;
for a short code, print I<lm_shortmsg>B<[>I<status.outcome>B<]>.

=item B<int> I<status.userbreak>

Set when termination has been forced by the user-supplied routine I<evaluate>.

=back

=back


=head1 NOTES

=head2 Initializing parameter records.

The parameter record I<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 I<lm_control_struct>.

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

=head2 More examples

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

=head1 COPYING

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

Software: FreeBSD License

Documentation: Creative Commons Attribution Share Alike


=head1 SEE ALSO

=begin html

<a href="http://apps.jcns.fz-juelich.de/man/lmcurve.html"><b>lmcurve</b>(3)</a>

=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 <j.wuttke@fz-juelich.de>.