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