1=pod
2
3=begin html
4
5<link rel="stylesheet" href="podstyle.css" type="text/css" />
6
7=end html
8
9=head1 NAME
10
11lmmin - Levenberg-Marquardt least-squares minimization
12
13
14=head1 SYNOPSIS
15
16B<#include <lmmin.h>>
17
18B<void lmmin( const int> I<n_par>B<, double *>I<par>B<, const int> I<m_dat>B<,
19            constS< >void *>I<data>B<,
20            void *>I<evaluate>B<(
21                 constS< >double *>I<par>B<, const int >I<m_dat>B<,
22                 constS< >void *>I<data>B<, double *>I<fvec>B<, int *>I<userbreak>B<),
23            constS< >lm_control_struct *>I<control>B<,
24            lm_status_struct *>I<status>B< );>
25
26B<extern const lm_control_struct lm_control_double;>
27
28B<extern const lm_control_struct lm_control_float;>
29
30B<extern const char *lm_infmsg[];>
31
32B<extern const char *lm_shortmsg[];>
33
34=head1 DESCRIPTION
35
36B<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>().
37On success, I<par> represents a local minimum, not necessarily a global one; it may depend on its starting value.
38
39For applications in curve fitting, the wrapper function B<lmcurve(3)> offers a simplified API.
40
41The 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.
42
43Function arguments:
44
45=over
46
47=item I<n_par>
48
49Number of free variables.
50Length of parameter vector I<par>.
51
52=item I<par>
53
54Parameter vector.
55On input, it must contain a reasonable guess.
56On output, it contains the solution found to minimize ||I<fvec>||.
57
58=item I<m_dat>
59
60Length of vector I<fvec>.
61Must statisfy I<n_par> <= I<m_dat>.
62
63=item I<data>
64
65This pointer is ignored by the fit algorithm,
66except for appearing as an argument in all calls to the user-supplied
67routine I<evaluate>.
68
69=item I<evaluate>
70
71Pointer to a user-supplied function that computes I<m_dat> elements of vector I<fvec> for a given parameter vector I<par>.
72If I<evaluate> return with *I<userbreak> set to a negative value, B<lmmin()> will interrupt the fitting and terminate.
73
74=item I<control>
75
76Parameter collection for tuning the fit procedure.
77In most cases, the default &I<lm_control_double> is adequate.
78If I<f> is only computed with single-precision accuracy,
79I<&lm_control_float> should be used.
80See also below, NOTES on initializing parameter records.
81
82I<control> has the following members (for more details, see the source file I<lmstruct.h>):
83
84=over
85
86=item B<double> I<control.ftol>
87
88Relative error desired in the sum of squares.
89Recommended setting: somewhat above machine precision; less if I<fvec> is computed with reduced accuracy.
90
91=item B<double> I<control.xtol>
92
93Relative error between last two approximations.
94Recommended setting: as I<ftol>.
95
96=item B<double> I<control.gtol>
97
98A measure for degeneracy.
99Recommended setting: as I<ftol>.
100
101=item B<double> I<control.epsilon>
102
103Step used to calculate the Jacobian.
104Recommended setting: as I<ftol>, but definitely less than the accuracy of I<fvec>.
105
106=item B<double> I<control.stepbound>
107
108Initial bound to steps in the outer loop, generally between 0.01 and 100; recommended value is 100.
109
110=item B<int> I<control.patience>
111
112Used to set the maximum number of function evaluations to patience*n_par.
113
114=item B<int> I<control.scale_diag>
115
116Logical switch (0 or 1).
117If 1, then scale parameters to their initial value.
118This is the recommended setting.
119
120=item B<FILE*> I<control.msgfile>
121
122Progress messages will be written to this file.
123Typically I<stdout> or I<stderr>.
124The value I<NULL> will be interpreted as I<stdout>.
125
126=item B<int> I<control.verbosity>
127
128If nonzero, some progress information from within the LM algorithm
129is written to control.stream.
130
131=item B<int> I<control.n_maxpri>
132
133-1, or maximum number of parameters to print.
134
135=item B<int> I<control.m_maxpri>
136
137-1, or maximum number of residuals to print.
138
139=back
140
141=item I<status>
142
143A record used to return information about the minimization process:
144
145=over
146
147=item B<double> I<status.fnorm>
148
149Norm of the vector I<fvec>;
150
151=item B<int> I<status.nfev>
152
153Actual number of iterations;
154
155=item B<int> I<status.outcome>
156
157Status of minimization;
158for the corresponding text message, print I<lm_infmsg>B<[>I<status.outcome>B<]>;
159for a short code, print I<lm_shortmsg>B<[>I<status.outcome>B<]>.
160
161=item B<int> I<status.userbreak>
162
163Set when termination has been forced by the user-supplied routine I<evaluate>.
164
165=back
166
167=back
168
169
170=head1 NOTES
171
172=head2 Initializing parameter records.
173
174The parameter record I<control> should always be initialized
175from supplied default records:
176
177    lm_control_struct control = lm_control_double; /* or _float */
178
179After this, parameters may be overwritten:
180
181    control.patience = 500; /* allow more iterations */
182    control.verbosity = 15; /* for verbose monitoring */
183
184An application written this way is guaranteed to work even if new parameters
185are added to I<lm_control_struct>.
186
187Conversely, addition of parameters is not considered an API change; it may happen without increment of the major version number.
188
189=head1 EXAMPLES
190
191=head2 Fitting a surface
192
193Fit a data set y(t) by a function f(t;p) where t is a two-dimensional vector:
194
195    #include "lmmin.h"
196    #include <stdio.h>
197
198    /* fit model: a plane p0 + p1*tx + p2*tz */
199    double f( double tx, double tz, const double *p )
200    {
201        return p[0] + p[1]*tx + p[2]*tz;
202    }
203
204    /* data structure to transmit data arays and fit model */
205    typedef struct {
206        double *tx, *tz;
207        double *y;
208        double (*f)( double tx, double tz, const double *p );
209    } data_struct;
210
211    /* function evaluation, determination of residues */
212    void evaluate_surface( const double *par, int m_dat,
213        const void *data, double *fvec, int *userbreak )
214    {
215        /* for readability, explicit type conversion */
216        data_struct *D;
217        D = (data_struct*)data;
218
219        int i;
220        for ( i = 0; i < m_dat; i++ )
221    	fvec[i] = D->y[i] - D->f( D->tx[i], D->tz[i], par );
222    }
223
224    int main()
225    {
226        /* parameter vector */
227        int n_par = 3; /* number of parameters in model function f */
228        double par[3] = { -1, 0, 1 }; /* arbitrary starting value */
229
230        /* data points */
231        int m_dat = 4;
232        double tx[4] = { -1, -1,  1,  1 };
233        double tz[4] = { -1,  1, -1,  1 };
234        double y[4]  = {  0,  1,  1,  2 };
235
236        data_struct data = { tx, tz, y, f };
237
238        /* auxiliary parameters */
239        lm_status_struct status;
240        lm_control_struct control = lm_control_double;
241        control.verbosity = 3;
242
243        /* perform the fit */
244        printf( "Fitting:\n" );
245        lmmin( n_par, par, m_dat, (const void*) &data, evaluate_surface,
246               &control, &status );
247
248        /* print results */
249        printf( "\nResults:\n" );
250        printf( "status after %d function evaluations:\n  %s\n",
251                status.nfev, lm_infmsg[status.outcome] );
252
253        printf("obtained parameters:\n");
254        int i;
255        for ( i=0; i<n_par; ++i )
256    	printf("  par[%i] = %12g\n", i, par[i]);
257        printf("obtained norm:\n  %12g\n", status.fnorm );
258
259        printf("fitting data as follows:\n");
260        double ff;
261        for ( i=0; i<m_dat; ++i ){
262            ff = f(tx[i], tz[i], par);
263            printf( "  t[%2d]=%12g,%12g y=%12g fit=%12g residue=%12g\n",
264                    i, tx[i], tz[i], y[i], ff, y[i] - ff );
265        }
266
267        return 0;
268    }
269
270=head2 More examples
271
272For more examples, see the homepage and directories demo/ and test/ in the source distribution.
273
274=head1 COPYING
275
276Copyright (C):
277   1980-1999 University of Chicago
278   2004-2015 Joachim Wuttke, Forschungszentrum Juelich GmbH
279
280Software: FreeBSD License
281
282Documentation: Creative Commons Attribution Share Alike
283
284
285=head1 SEE ALSO
286
287=begin html
288
289<a href="http://apps.jcns.fz-juelich.de/man/lmcurve.html"><b>lmcurve</b>(3)</a>
290
291=end html
292
293=begin man
294
295\fBlmcurve\fR(3)
296.PP
297
298=end man
299
300Homepage: http://apps.jcns.fz-juelich.de/lmfit
301
302=head1 BUGS
303
304Please send bug reports and suggestions to the author <j.wuttke@fz-juelich.de>.
305