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