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 11 lmmin - Levenberg-Marquardt least-squares minimization 12 13 14 =head1 SYNOPSIS 15 16 B<#include <lmmin.h>> 17 18 B<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 26 B<extern const lm_control_struct lm_control_double;> 27 28 B<extern const lm_control_struct lm_control_float;> 29 30 B<extern const char *lm_infmsg[];> 31 32 B<extern const char *lm_shortmsg[];> 33 34 =head1 DESCRIPTION 35 36 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>(). 37 On success, I<par> represents a local minimum, not necessarily a global one; it may depend on its starting value. 38 39 For applications in curve fitting, the wrapper function B<lmcurve(3)> offers a simplified API. 40 41 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. 42 43 Function arguments: 44 45 =over 46 47 =item I<n_par> 48 49 Number of free variables. 50 Length of parameter vector I<par>. 51 52 =item I<par> 53 54 Parameter vector. 55 On input, it must contain a reasonable guess. 56 On output, it contains the solution found to minimize ||I<fvec>||. 57 58 =item I<m_dat> 59 60 Length of vector I<fvec>. 61 Must statisfy I<n_par> <= I<m_dat>. 62 63 =item I<data> 64 65 This pointer is ignored by the fit algorithm, 66 except for appearing as an argument in all calls to the user-supplied 67 routine I<evaluate>. 68 69 =item I<evaluate> 70 71 Pointer to a user-supplied function that computes I<m_dat> elements of vector I<fvec> for a given parameter vector I<par>. 72 If 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 76 Parameter collection for tuning the fit procedure. 77 In most cases, the default &I<lm_control_double> is adequate. 78 If I<f> is only computed with single-precision accuracy, 79 I<&lm_control_float> should be used. 80 See also below, NOTES on initializing parameter records. 81 82 I<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 88 Relative error desired in the sum of squares. 89 Recommended setting: somewhat above machine precision; less if I<fvec> is computed with reduced accuracy. 90 91 =item B<double> I<control.xtol> 92 93 Relative error between last two approximations. 94 Recommended setting: as I<ftol>. 95 96 =item B<double> I<control.gtol> 97 98 A measure for degeneracy. 99 Recommended setting: as I<ftol>. 100 101 =item B<double> I<control.epsilon> 102 103 Step used to calculate the Jacobian. 104 Recommended setting: as I<ftol>, but definitely less than the accuracy of I<fvec>. 105 106 =item B<double> I<control.stepbound> 107 108 Initial 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 112 Used to set the maximum number of function evaluations to patience*n_par. 113 114 =item B<int> I<control.scale_diag> 115 116 Logical switch (0 or 1). 117 If 1, then scale parameters to their initial value. 118 This is the recommended setting. 119 120 =item B<FILE*> I<control.msgfile> 121 122 Progress messages will be written to this file. 123 Typically I<stdout> or I<stderr>. 124 The value I<NULL> will be interpreted as I<stdout>. 125 126 =item B<int> I<control.verbosity> 127 128 If nonzero, some progress information from within the LM algorithm 129 is 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 143 A record used to return information about the minimization process: 144 145 =over 146 147 =item B<double> I<status.fnorm> 148 149 Norm of the vector I<fvec>; 150 151 =item B<int> I<status.nfev> 152 153 Actual number of iterations; 154 155 =item B<int> I<status.outcome> 156 157 Status of minimization; 158 for the corresponding text message, print I<lm_infmsg>B<[>I<status.outcome>B<]>; 159 for a short code, print I<lm_shortmsg>B<[>I<status.outcome>B<]>. 160 161 =item B<int> I<status.userbreak> 162 163 Set 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 174 The parameter record I<control> should always be initialized 175 from supplied default records: 176 177 lm_control_struct control = lm_control_double; /* or _float */ 178 179 After this, parameters may be overwritten: 180 181 control.patience = 500; /* allow more iterations */ 182 control.verbosity = 15; /* for verbose monitoring */ 183 184 An application written this way is guaranteed to work even if new parameters 185 are added to I<lm_control_struct>. 186 187 Conversely, 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 193 Fit 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 272 For more examples, see the homepage and directories demo/ and test/ in the source distribution. 273 274 =head1 COPYING 275 276 Copyright (C): 277 1980-1999 University of Chicago 278 2004-2015 Joachim Wuttke, Forschungszentrum Juelich GmbH 279 280 Software: FreeBSD License 281 282 Documentation: 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 300 Homepage: http://apps.jcns.fz-juelich.de/lmfit 301 302 =head1 BUGS 303 304 Please send bug reports and suggestions to the author <j.wuttke (a] fz-juelich.de>. 305
