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