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