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 
37 /**
38   * \ingroup NonLinearOptimization_Module
39   * \brief Performs non linear optimization over a non-linear function,
40   * using a variant of the Levenberg Marquardt algorithm.
41   *
42   * Check wikipedia for more information.
43   * http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
44   */
45 template<typename FunctorType, typename Scalar=double>
46 class LevenbergMarquardt
47 {
48 public:
LevenbergMarquardt(FunctorType & _functor)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 {
ParametersParameters55         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 
resetParameters(void)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 
lm_param(void)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
lmder1(FVectorType & x,const Scalar tol)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
minimize(FVectorType & x)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
minimizeInit(FVectorType & x)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
minimizeOneStep(FVectorType & x)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
lmstr1(FVectorType & x,const Scalar tol)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
minimizeOptimumStorageInit(FVectorType & x)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
minimizeOptimumStorageOneStep(FVectorType & x)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
minimizeOptimumStorage(FVectorType & x)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
lmdif1(FunctorType & functor,FVectorType & x,Index * nfev,const Scalar tol)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
651