All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Groups Pages
NonLinearOptimization/LevenbergMarquardt.h
1 // -*- coding: utf-8
2 // vim: set fileencoding=utf-8
3 
4 // This file is part of Eigen, a lightweight C++ template library
5 // for linear algebra.
6 //
7 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
8 //
9 // This Source Code Form is subject to the terms of the Mozilla
10 // Public License v. 2.0. If a copy of the MPL was not distributed
11 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
12 
13 #ifndef EIGEN_LEVENBERGMARQUARDT__H
14 #define EIGEN_LEVENBERGMARQUARDT__H
15 
16 namespace Eigen {
17 
18 namespace LevenbergMarquardtSpace {
19  enum Status {
20  NotStarted = -2,
21  Running = -1,
22  ImproperInputParameters = 0,
23  RelativeReductionTooSmall = 1,
24  RelativeErrorTooSmall = 2,
25  RelativeErrorAndReductionTooSmall = 3,
26  CosinusTooSmall = 4,
27  TooManyFunctionEvaluation = 5,
28  FtolTooSmall = 6,
29  XtolTooSmall = 7,
30  GtolTooSmall = 8,
31  UserAsked = 9
32  };
33 }
34 
35 
36 
45 template<typename FunctorType, typename Scalar=double>
46 class LevenbergMarquardt
47 {
48  static Scalar sqrt_epsilon()
49  {
50  using std::sqrt;
51  return sqrt(NumTraits<Scalar>::epsilon());
52  }
53 
54 public:
55  LevenbergMarquardt(FunctorType &_functor)
56  : functor(_functor) { nfev = njev = iter = 0; fnorm = gnorm = 0.; useExternalScaling=false; }
57 
58  typedef DenseIndex Index;
59 
60  struct Parameters {
61  Parameters()
62  : factor(Scalar(100.))
63  , maxfev(400)
64  , ftol(sqrt_epsilon())
65  , xtol(sqrt_epsilon())
66  , gtol(Scalar(0.))
67  , epsfcn(Scalar(0.)) {}
68  Scalar factor;
69  Index maxfev; // maximum number of function evaluation
70  Scalar ftol;
71  Scalar xtol;
72  Scalar gtol;
73  Scalar epsfcn;
74  };
75 
76  typedef Matrix< Scalar, Dynamic, 1 > FVectorType;
77  typedef Matrix< Scalar, Dynamic, Dynamic > JacobianType;
78 
79  LevenbergMarquardtSpace::Status lmder1(
80  FVectorType &x,
81  const Scalar tol = sqrt_epsilon()
82  );
83 
84  LevenbergMarquardtSpace::Status minimize(FVectorType &x);
85  LevenbergMarquardtSpace::Status minimizeInit(FVectorType &x);
86  LevenbergMarquardtSpace::Status minimizeOneStep(FVectorType &x);
87 
88  static LevenbergMarquardtSpace::Status lmdif1(
89  FunctorType &functor,
90  FVectorType &x,
91  Index *nfev,
92  const Scalar tol = sqrt_epsilon()
93  );
94 
95  LevenbergMarquardtSpace::Status lmstr1(
96  FVectorType &x,
97  const Scalar tol = sqrt_epsilon()
98  );
99 
100  LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType &x);
101  LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(FVectorType &x);
102  LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(FVectorType &x);
103 
104  void resetParameters(void) { parameters = Parameters(); }
105 
106  Parameters parameters;
107  FVectorType fvec, qtf, diag;
108  JacobianType fjac;
109  PermutationMatrix<Dynamic,Dynamic> permutation;
110  Index nfev;
111  Index njev;
112  Index iter;
113  Scalar fnorm, gnorm;
114  bool useExternalScaling;
115 
116  Scalar lm_param(void) { return par; }
117 private:
118 
119  FunctorType &functor;
120  Index n;
121  Index m;
122  FVectorType wa1, wa2, wa3, wa4;
123 
124  Scalar par, sum;
125  Scalar temp, temp1, temp2;
126  Scalar delta;
127  Scalar ratio;
128  Scalar pnorm, xnorm, fnorm1, actred, dirder, prered;
129 
130  LevenbergMarquardt& operator=(const LevenbergMarquardt&);
131 };
132 
133 template<typename FunctorType, typename Scalar>
134 LevenbergMarquardtSpace::Status
135 LevenbergMarquardt<FunctorType,Scalar>::lmder1(
136  FVectorType &x,
137  const Scalar tol
138  )
139 {
140  n = x.size();
141  m = functor.values();
142 
143  /* check the input parameters for errors. */
144  if (n <= 0 || m < n || tol < 0.)
145  return LevenbergMarquardtSpace::ImproperInputParameters;
146 
147  resetParameters();
148  parameters.ftol = tol;
149  parameters.xtol = tol;
150  parameters.maxfev = 100*(n+1);
151 
152  return minimize(x);
153 }
154 
155 
156 template<typename FunctorType, typename Scalar>
157 LevenbergMarquardtSpace::Status
158 LevenbergMarquardt<FunctorType,Scalar>::minimize(FVectorType &x)
159 {
160  LevenbergMarquardtSpace::Status status = minimizeInit(x);
161  if (status==LevenbergMarquardtSpace::ImproperInputParameters)
162  return status;
163  do {
164  status = minimizeOneStep(x);
165  } while (status==LevenbergMarquardtSpace::Running);
166  return status;
167 }
168 
169 template<typename FunctorType, typename Scalar>
170 LevenbergMarquardtSpace::Status
171 LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(FVectorType &x)
172 {
173  n = x.size();
174  m = functor.values();
175 
176  wa1.resize(n); wa2.resize(n); wa3.resize(n);
177  wa4.resize(m);
178  fvec.resize(m);
179  fjac.resize(m, n);
180  if (!useExternalScaling)
181  diag.resize(n);
182  eigen_assert( (!useExternalScaling || diag.size()==n) && "When useExternalScaling is set, the caller must provide a valid 'diag'");
183  qtf.resize(n);
184 
185  /* Function Body */
186  nfev = 0;
187  njev = 0;
188 
189  /* check the input parameters for errors. */
190  if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
191  return LevenbergMarquardtSpace::ImproperInputParameters;
192 
193  if (useExternalScaling)
194  for (Index j = 0; j < n; ++j)
195  if (diag[j] <= 0.)
196  return LevenbergMarquardtSpace::ImproperInputParameters;
197 
198  /* evaluate the function at the starting point */
199  /* and calculate its norm. */
200  nfev = 1;
201  if ( functor(x, fvec) < 0)
202  return LevenbergMarquardtSpace::UserAsked;
203  fnorm = fvec.stableNorm();
204 
205  /* initialize levenberg-marquardt parameter and iteration counter. */
206  par = 0.;
207  iter = 1;
208 
209  return LevenbergMarquardtSpace::NotStarted;
210 }
211 
212 template<typename FunctorType, typename Scalar>
213 LevenbergMarquardtSpace::Status
214 LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType &x)
215 {
216  using std::abs;
217  using std::sqrt;
218 
219  eigen_assert(x.size()==n); // check the caller is not cheating us
220 
221  /* calculate the jacobian matrix. */
222  Index df_ret = functor.df(x, fjac);
223  if (df_ret<0)
224  return LevenbergMarquardtSpace::UserAsked;
225  if (df_ret>0)
226  // numerical diff, we evaluated the function df_ret times
227  nfev += df_ret;
228  else njev++;
229 
230  /* compute the qr factorization of the jacobian. */
231  wa2 = fjac.colwise().blueNorm();
232  ColPivHouseholderQR<JacobianType> qrfac(fjac);
233  fjac = qrfac.matrixQR();
234  permutation = qrfac.colsPermutation();
235 
236  /* on the first iteration and if external scaling is not used, scale according */
237  /* to the norms of the columns of the initial jacobian. */
238  if (iter == 1) {
239  if (!useExternalScaling)
240  for (Index j = 0; j < n; ++j)
241  diag[j] = (wa2[j]==0.)? 1. : wa2[j];
242 
243  /* on the first iteration, calculate the norm of the scaled x */
244  /* and initialize the step bound delta. */
245  xnorm = diag.cwiseProduct(x).stableNorm();
246  delta = parameters.factor * xnorm;
247  if (delta == 0.)
248  delta = parameters.factor;
249  }
250 
251  /* form (q transpose)*fvec and store the first n components in */
252  /* qtf. */
253  wa4 = fvec;
254  wa4.applyOnTheLeft(qrfac.householderQ().adjoint());
255  qtf = wa4.head(n);
256 
257  /* compute the norm of the scaled gradient. */
258  gnorm = 0.;
259  if (fnorm != 0.)
260  for (Index j = 0; j < n; ++j)
261  if (wa2[permutation.indices()[j]] != 0.)
262  gnorm = (std::max)(gnorm, abs( fjac.col(j).head(j+1).dot(qtf.head(j+1)/fnorm) / wa2[permutation.indices()[j]]));
263 
264  /* test for convergence of the gradient norm. */
265  if (gnorm <= parameters.gtol)
266  return LevenbergMarquardtSpace::CosinusTooSmall;
267 
268  /* rescale if necessary. */
269  if (!useExternalScaling)
270  diag = diag.cwiseMax(wa2);
271 
272  do {
273 
274  /* determine the levenberg-marquardt parameter. */
275  internal::lmpar2<Scalar>(qrfac, diag, qtf, delta, par, wa1);
276 
277  /* store the direction p and x + p. calculate the norm of p. */
278  wa1 = -wa1;
279  wa2 = x + wa1;
280  pnorm = diag.cwiseProduct(wa1).stableNorm();
281 
282  /* on the first iteration, adjust the initial step bound. */
283  if (iter == 1)
284  delta = (std::min)(delta,pnorm);
285 
286  /* evaluate the function at x + p and calculate its norm. */
287  if ( functor(wa2, wa4) < 0)
288  return LevenbergMarquardtSpace::UserAsked;
289  ++nfev;
290  fnorm1 = wa4.stableNorm();
291 
292  /* compute the scaled actual reduction. */
293  actred = -1.;
294  if (Scalar(.1) * fnorm1 < fnorm)
295  actred = 1. - numext::abs2(fnorm1 / fnorm);
296 
297  /* compute the scaled predicted reduction and */
298  /* the scaled directional derivative. */
299  wa3 = fjac.template triangularView<Upper>() * (qrfac.colsPermutation().inverse() *wa1);
300  temp1 = numext::abs2(wa3.stableNorm() / fnorm);
301  temp2 = numext::abs2(sqrt(par) * pnorm / fnorm);
302  prered = temp1 + temp2 / Scalar(.5);
303  dirder = -(temp1 + temp2);
304 
305  /* compute the ratio of the actual to the predicted */
306  /* reduction. */
307  ratio = 0.;
308  if (prered != 0.)
309  ratio = actred / prered;
310 
311  /* update the step bound. */
312  if (ratio <= Scalar(.25)) {
313  if (actred >= 0.)
314  temp = Scalar(.5);
315  if (actred < 0.)
316  temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
317  if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1))
318  temp = Scalar(.1);
319  /* Computing MIN */
320  delta = temp * (std::min)(delta, pnorm / Scalar(.1));
321  par /= temp;
322  } else if (!(par != 0. && ratio < Scalar(.75))) {
323  delta = pnorm / Scalar(.5);
324  par = Scalar(.5) * par;
325  }
326 
327  /* test for successful iteration. */
328  if (ratio >= Scalar(1e-4)) {
329  /* successful iteration. update x, fvec, and their norms. */
330  x = wa2;
331  wa2 = diag.cwiseProduct(x);
332  fvec = wa4;
333  xnorm = wa2.stableNorm();
334  fnorm = fnorm1;
335  ++iter;
336  }
337 
338  /* tests for convergence. */
339  if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm)
340  return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall;
341  if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.)
342  return LevenbergMarquardtSpace::RelativeReductionTooSmall;
343  if (delta <= parameters.xtol * xnorm)
344  return LevenbergMarquardtSpace::RelativeErrorTooSmall;
345 
346  /* tests for termination and stringent tolerances. */
347  if (nfev >= parameters.maxfev)
348  return LevenbergMarquardtSpace::TooManyFunctionEvaluation;
349  if (abs(actred) <= NumTraits<Scalar>::epsilon() && prered <= NumTraits<Scalar>::epsilon() && Scalar(.5) * ratio <= 1.)
350  return LevenbergMarquardtSpace::FtolTooSmall;
351  if (delta <= NumTraits<Scalar>::epsilon() * xnorm)
352  return LevenbergMarquardtSpace::XtolTooSmall;
353  if (gnorm <= NumTraits<Scalar>::epsilon())
354  return LevenbergMarquardtSpace::GtolTooSmall;
355 
356  } while (ratio < Scalar(1e-4));
357 
358  return LevenbergMarquardtSpace::Running;
359 }
360 
361 template<typename FunctorType, typename Scalar>
362 LevenbergMarquardtSpace::Status
363 LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
364  FVectorType &x,
365  const Scalar tol
366  )
367 {
368  n = x.size();
369  m = functor.values();
370 
371  /* check the input parameters for errors. */
372  if (n <= 0 || m < n || tol < 0.)
373  return LevenbergMarquardtSpace::ImproperInputParameters;
374 
375  resetParameters();
376  parameters.ftol = tol;
377  parameters.xtol = tol;
378  parameters.maxfev = 100*(n+1);
379 
380  return minimizeOptimumStorage(x);
381 }
382 
383 template<typename FunctorType, typename Scalar>
384 LevenbergMarquardtSpace::Status
385 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(FVectorType &x)
386 {
387  n = x.size();
388  m = functor.values();
389 
390  wa1.resize(n); wa2.resize(n); wa3.resize(n);
391  wa4.resize(m);
392  fvec.resize(m);
393  // Only R is stored in fjac. Q is only used to compute 'qtf', which is
394  // Q.transpose()*rhs. qtf will be updated using givens rotation,
395  // instead of storing them in Q.
396  // The purpose it to only use a nxn matrix, instead of mxn here, so
397  // that we can handle cases where m>>n :
398  fjac.resize(n, n);
399  if (!useExternalScaling)
400  diag.resize(n);
401  eigen_assert( (!useExternalScaling || diag.size()==n) && "When useExternalScaling is set, the caller must provide a valid 'diag'");
402  qtf.resize(n);
403 
404  /* Function Body */
405  nfev = 0;
406  njev = 0;
407 
408  /* check the input parameters for errors. */
409  if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
410  return LevenbergMarquardtSpace::ImproperInputParameters;
411 
412  if (useExternalScaling)
413  for (Index j = 0; j < n; ++j)
414  if (diag[j] <= 0.)
415  return LevenbergMarquardtSpace::ImproperInputParameters;
416 
417  /* evaluate the function at the starting point */
418  /* and calculate its norm. */
419  nfev = 1;
420  if ( functor(x, fvec) < 0)
421  return LevenbergMarquardtSpace::UserAsked;
422  fnorm = fvec.stableNorm();
423 
424  /* initialize levenberg-marquardt parameter and iteration counter. */
425  par = 0.;
426  iter = 1;
427 
428  return LevenbergMarquardtSpace::NotStarted;
429 }
430 
431 
432 template<typename FunctorType, typename Scalar>
433 LevenbergMarquardtSpace::Status
434 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(FVectorType &x)
435 {
436  using std::abs;
437  using std::sqrt;
438 
439  eigen_assert(x.size()==n); // check the caller is not cheating us
440 
441  Index i, j;
442  bool sing;
443 
444  /* compute the qr factorization of the jacobian matrix */
445  /* calculated one row at a time, while simultaneously */
446  /* forming (q transpose)*fvec and storing the first */
447  /* n components in qtf. */
448  qtf.fill(0.);
449  fjac.fill(0.);
450  Index rownb = 2;
451  for (i = 0; i < m; ++i) {
452  if (functor.df(x, wa3, rownb) < 0) return LevenbergMarquardtSpace::UserAsked;
453  internal::rwupdt<Scalar>(fjac, wa3, qtf, fvec[i]);
454  ++rownb;
455  }
456  ++njev;
457 
458  /* if the jacobian is rank deficient, call qrfac to */
459  /* reorder its columns and update the components of qtf. */
460  sing = false;
461  for (j = 0; j < n; ++j) {
462  if (fjac(j,j) == 0.)
463  sing = true;
464  wa2[j] = fjac.col(j).head(j).stableNorm();
465  }
466  permutation.setIdentity(n);
467  if (sing) {
468  wa2 = fjac.colwise().blueNorm();
469  // TODO We have no unit test covering this code path, do not modify
470  // until it is carefully tested
471  ColPivHouseholderQR<JacobianType> qrfac(fjac);
472  fjac = qrfac.matrixQR();
473  wa1 = fjac.diagonal();
474  fjac.diagonal() = qrfac.hCoeffs();
475  permutation = qrfac.colsPermutation();
476  // TODO : avoid this:
477  for(Index ii=0; ii< fjac.cols(); ii++) fjac.col(ii).segment(ii+1, fjac.rows()-ii-1) *= fjac(ii,ii); // rescale vectors
478 
479  for (j = 0; j < n; ++j) {
480  if (fjac(j,j) != 0.) {
481  sum = 0.;
482  for (i = j; i < n; ++i)
483  sum += fjac(i,j) * qtf[i];
484  temp = -sum / fjac(j,j);
485  for (i = j; i < n; ++i)
486  qtf[i] += fjac(i,j) * temp;
487  }
488  fjac(j,j) = wa1[j];
489  }
490  }
491 
492  /* on the first iteration and if external scaling is not used, scale according */
493  /* to the norms of the columns of the initial jacobian. */
494  if (iter == 1) {
495  if (!useExternalScaling)
496  for (j = 0; j < n; ++j)
497  diag[j] = (wa2[j]==0.)? 1. : wa2[j];
498 
499  /* on the first iteration, calculate the norm of the scaled x */
500  /* and initialize the step bound delta. */
501  xnorm = diag.cwiseProduct(x).stableNorm();
502  delta = parameters.factor * xnorm;
503  if (delta == 0.)
504  delta = parameters.factor;
505  }
506 
507  /* compute the norm of the scaled gradient. */
508  gnorm = 0.;
509  if (fnorm != 0.)
510  for (j = 0; j < n; ++j)
511  if (wa2[permutation.indices()[j]] != 0.)
512  gnorm = (std::max)(gnorm, abs( fjac.col(j).head(j+1).dot(qtf.head(j+1)/fnorm) / wa2[permutation.indices()[j]]));
513 
514  /* test for convergence of the gradient norm. */
515  if (gnorm <= parameters.gtol)
516  return LevenbergMarquardtSpace::CosinusTooSmall;
517 
518  /* rescale if necessary. */
519  if (!useExternalScaling)
520  diag = diag.cwiseMax(wa2);
521 
522  do {
523 
524  /* determine the levenberg-marquardt parameter. */
525  internal::lmpar<Scalar>(fjac, permutation.indices(), diag, qtf, delta, par, wa1);
526 
527  /* store the direction p and x + p. calculate the norm of p. */
528  wa1 = -wa1;
529  wa2 = x + wa1;
530  pnorm = diag.cwiseProduct(wa1).stableNorm();
531 
532  /* on the first iteration, adjust the initial step bound. */
533  if (iter == 1)
534  delta = (std::min)(delta,pnorm);
535 
536  /* evaluate the function at x + p and calculate its norm. */
537  if ( functor(wa2, wa4) < 0)
538  return LevenbergMarquardtSpace::UserAsked;
539  ++nfev;
540  fnorm1 = wa4.stableNorm();
541 
542  /* compute the scaled actual reduction. */
543  actred = -1.;
544  if (Scalar(.1) * fnorm1 < fnorm)
545  actred = 1. - numext::abs2(fnorm1 / fnorm);
546 
547  /* compute the scaled predicted reduction and */
548  /* the scaled directional derivative. */
549  wa3 = fjac.topLeftCorner(n,n).template triangularView<Upper>() * (permutation.inverse() * wa1);
550  temp1 = numext::abs2(wa3.stableNorm() / fnorm);
551  temp2 = numext::abs2(sqrt(par) * pnorm / fnorm);
552  prered = temp1 + temp2 / Scalar(.5);
553  dirder = -(temp1 + temp2);
554 
555  /* compute the ratio of the actual to the predicted */
556  /* reduction. */
557  ratio = 0.;
558  if (prered != 0.)
559  ratio = actred / prered;
560 
561  /* update the step bound. */
562  if (ratio <= Scalar(.25)) {
563  if (actred >= 0.)
564  temp = Scalar(.5);
565  if (actred < 0.)
566  temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
567  if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1))
568  temp = Scalar(.1);
569  /* Computing MIN */
570  delta = temp * (std::min)(delta, pnorm / Scalar(.1));
571  par /= temp;
572  } else if (!(par != 0. && ratio < Scalar(.75))) {
573  delta = pnorm / Scalar(.5);
574  par = Scalar(.5) * par;
575  }
576 
577  /* test for successful iteration. */
578  if (ratio >= Scalar(1e-4)) {
579  /* successful iteration. update x, fvec, and their norms. */
580  x = wa2;
581  wa2 = diag.cwiseProduct(x);
582  fvec = wa4;
583  xnorm = wa2.stableNorm();
584  fnorm = fnorm1;
585  ++iter;
586  }
587 
588  /* tests for convergence. */
589  if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm)
590  return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall;
591  if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.)
592  return LevenbergMarquardtSpace::RelativeReductionTooSmall;
593  if (delta <= parameters.xtol * xnorm)
594  return LevenbergMarquardtSpace::RelativeErrorTooSmall;
595 
596  /* tests for termination and stringent tolerances. */
597  if (nfev >= parameters.maxfev)
598  return LevenbergMarquardtSpace::TooManyFunctionEvaluation;
599  if (abs(actred) <= NumTraits<Scalar>::epsilon() && prered <= NumTraits<Scalar>::epsilon() && Scalar(.5) * ratio <= 1.)
600  return LevenbergMarquardtSpace::FtolTooSmall;
601  if (delta <= NumTraits<Scalar>::epsilon() * xnorm)
602  return LevenbergMarquardtSpace::XtolTooSmall;
603  if (gnorm <= NumTraits<Scalar>::epsilon())
604  return LevenbergMarquardtSpace::GtolTooSmall;
605 
606  } while (ratio < Scalar(1e-4));
607 
608  return LevenbergMarquardtSpace::Running;
609 }
610 
611 template<typename FunctorType, typename Scalar>
612 LevenbergMarquardtSpace::Status
613 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(FVectorType &x)
614 {
615  LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x);
616  if (status==LevenbergMarquardtSpace::ImproperInputParameters)
617  return status;
618  do {
619  status = minimizeOptimumStorageOneStep(x);
620  } while (status==LevenbergMarquardtSpace::Running);
621  return status;
622 }
623 
624 template<typename FunctorType, typename Scalar>
625 LevenbergMarquardtSpace::Status
626 LevenbergMarquardt<FunctorType,Scalar>::lmdif1(
627  FunctorType &functor,
628  FVectorType &x,
629  Index *nfev,
630  const Scalar tol
631  )
632 {
633  Index n = x.size();
634  Index m = functor.values();
635 
636  /* check the input parameters for errors. */
637  if (n <= 0 || m < n || tol < 0.)
638  return LevenbergMarquardtSpace::ImproperInputParameters;
639 
640  NumericalDiff<FunctorType> numDiff(functor);
641  // embedded LevenbergMarquardt
642  LevenbergMarquardt<NumericalDiff<FunctorType>, Scalar > lm(numDiff);
643  lm.parameters.ftol = tol;
644  lm.parameters.xtol = tol;
645  lm.parameters.maxfev = 200*(n+1);
646 
647  LevenbergMarquardtSpace::Status info = LevenbergMarquardtSpace::Status(lm.minimize(x));
648  if (nfev)
649  * nfev = lm.nfev;
650  return info;
651 }
652 
653 } // end namespace Eigen
654 
655 #endif // EIGEN_LEVENBERGMARQUARDT__H
656 
657 //vim: ai ts=4 sts=4 et sw=4
Index maxfev() const
Definition: LevenbergMarquardt/LevenbergMarquardt.h:194
RealScalar gtol() const
Definition: LevenbergMarquardt/LevenbergMarquardt.h:185
RealScalar fnorm()
Definition: LevenbergMarquardt/LevenbergMarquardt.h:209
RealScalar ftol() const
Definition: LevenbergMarquardt/LevenbergMarquardt.h:182
FVectorType & fvec()
Definition: LevenbergMarquardt/LevenbergMarquardt.h:219
PermutationType permutation()
Definition: LevenbergMarquardt/LevenbergMarquardt.h:232
RealScalar xtol() const
Definition: LevenbergMarquardt/LevenbergMarquardt.h:179
Index nfev()
Definition: LevenbergMarquardt/LevenbergMarquardt.h:203
Index njev()
Definition: LevenbergMarquardt/LevenbergMarquardt.h:206
RealScalar lm_param(void)
Definition: LevenbergMarquardt/LevenbergMarquardt.h:215
void resetParameters()
Definition: LevenbergMarquardt/LevenbergMarquardt.h:145
RealScalar factor() const
Definition: LevenbergMarquardt/LevenbergMarquardt.h:188
FVectorType & diag()
Definition: LevenbergMarquardt/LevenbergMarquardt.h:197
RealScalar gnorm()
Definition: LevenbergMarquardt/LevenbergMarquardt.h:212