LevenbergMarquardt.h
Go to the documentation of this file.
00001 // -*- coding: utf-8
00002 // vim: set fileencoding=utf-8
00003 
00004 // This file is part of Eigen, a lightweight C++ template library
00005 // for linear algebra.
00006 //
00007 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
00008 //
00009 // This Source Code Form is subject to the terms of the Mozilla
00010 // Public License v. 2.0. If a copy of the MPL was not distributed
00011 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00012 
00013 #ifndef EIGEN_LEVENBERGMARQUARDT__H
00014 #define EIGEN_LEVENBERGMARQUARDT__H
00015 
00016 namespace Eigen { 
00017 
00018 namespace LevenbergMarquardtSpace {
00019     enum Status {
00020         NotStarted = -2,
00021         Running = -1,
00022         ImproperInputParameters = 0,
00023         RelativeReductionTooSmall = 1,
00024         RelativeErrorTooSmall = 2,
00025         RelativeErrorAndReductionTooSmall = 3,
00026         CosinusTooSmall = 4,
00027         TooManyFunctionEvaluation = 5,
00028         FtolTooSmall = 6,
00029         XtolTooSmall = 7,
00030         GtolTooSmall = 8,
00031         UserAsked = 9
00032     };
00033 }
00034 
00035 
00036 
00045 template<typename FunctorType, typename Scalar=double>
00046 class LevenbergMarquardt
00047 {
00048 public:
00049     LevenbergMarquardt(FunctorType &_functor)
00050         : functor(_functor) { nfev = njev = iter = 0;  fnorm = gnorm = 0.; useExternalScaling=false; }
00051 
00052     typedef DenseIndex Index;
00053 
00054     struct Parameters {
00055         Parameters()
00056             : factor(Scalar(100.))
00057             , maxfev(400)
00058             , ftol(std::sqrt(NumTraits<Scalar>::epsilon()))
00059             , xtol(std::sqrt(NumTraits<Scalar>::epsilon()))
00060             , gtol(Scalar(0.))
00061             , epsfcn(Scalar(0.)) {}
00062         Scalar factor;
00063         Index maxfev;   // maximum number of function evaluation
00064         Scalar ftol;
00065         Scalar xtol;
00066         Scalar gtol;
00067         Scalar epsfcn;
00068     };
00069 
00070     typedef Matrix< Scalar, Dynamic, 1 > FVectorType;
00071     typedef Matrix< Scalar, Dynamic, Dynamic > JacobianType;
00072 
00073     LevenbergMarquardtSpace::Status lmder1(
00074             FVectorType &x,
00075             const Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon())
00076             );
00077 
00078     LevenbergMarquardtSpace::Status minimize(FVectorType &x);
00079     LevenbergMarquardtSpace::Status minimizeInit(FVectorType &x);
00080     LevenbergMarquardtSpace::Status minimizeOneStep(FVectorType &x);
00081 
00082     static LevenbergMarquardtSpace::Status lmdif1(
00083             FunctorType &functor,
00084             FVectorType &x,
00085             Index *nfev,
00086             const Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon())
00087             );
00088 
00089     LevenbergMarquardtSpace::Status lmstr1(
00090             FVectorType  &x,
00091             const Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon())
00092             );
00093 
00094     LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType  &x);
00095     LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(FVectorType  &x);
00096     LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(FVectorType  &x);
00097 
00098     void resetParameters(void) { parameters = Parameters(); }
00099 
00100     Parameters parameters;
00101     FVectorType  fvec, qtf, diag;
00102     JacobianType fjac;
00103     PermutationMatrix<Dynamic,Dynamic> permutation;
00104     Index nfev;
00105     Index njev;
00106     Index iter;
00107     Scalar fnorm, gnorm;
00108     bool useExternalScaling; 
00109 
00110     Scalar lm_param(void) { return par; }
00111 private:
00112     FunctorType &functor;
00113     Index n;
00114     Index m;
00115     FVectorType wa1, wa2, wa3, wa4;
00116 
00117     Scalar par, sum;
00118     Scalar temp, temp1, temp2;
00119     Scalar delta;
00120     Scalar ratio;
00121     Scalar pnorm, xnorm, fnorm1, actred, dirder, prered;
00122 
00123     LevenbergMarquardt& operator=(const LevenbergMarquardt&);
00124 };
00125 
00126 template<typename FunctorType, typename Scalar>
00127 LevenbergMarquardtSpace::Status
00128 LevenbergMarquardt<FunctorType,Scalar>::lmder1(
00129         FVectorType  &x,
00130         const Scalar tol
00131         )
00132 {
00133     n = x.size();
00134     m = functor.values();
00135 
00136     /* check the input parameters for errors. */
00137     if (n <= 0 || m < n || tol < 0.)
00138         return LevenbergMarquardtSpace::ImproperInputParameters;
00139 
00140     resetParameters();
00141     parameters.ftol = tol;
00142     parameters.xtol = tol;
00143     parameters.maxfev = 100*(n+1);
00144 
00145     return minimize(x);
00146 }
00147 
00148 
00149 template<typename FunctorType, typename Scalar>
00150 LevenbergMarquardtSpace::Status
00151 LevenbergMarquardt<FunctorType,Scalar>::minimize(FVectorType  &x)
00152 {
00153     LevenbergMarquardtSpace::Status status = minimizeInit(x);
00154     if (status==LevenbergMarquardtSpace::ImproperInputParameters)
00155         return status;
00156     do {
00157         status = minimizeOneStep(x);
00158     } while (status==LevenbergMarquardtSpace::Running);
00159     return status;
00160 }
00161 
00162 template<typename FunctorType, typename Scalar>
00163 LevenbergMarquardtSpace::Status
00164 LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(FVectorType  &x)
00165 {
00166     n = x.size();
00167     m = functor.values();
00168 
00169     wa1.resize(n); wa2.resize(n); wa3.resize(n);
00170     wa4.resize(m);
00171     fvec.resize(m);
00172     fjac.resize(m, n);
00173     if (!useExternalScaling)
00174         diag.resize(n);
00175     eigen_assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
00176     qtf.resize(n);
00177 
00178     /* Function Body */
00179     nfev = 0;
00180     njev = 0;
00181 
00182     /*     check the input parameters for errors. */
00183     if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
00184         return LevenbergMarquardtSpace::ImproperInputParameters;
00185 
00186     if (useExternalScaling)
00187         for (Index j = 0; j < n; ++j)
00188             if (diag[j] <= 0.)
00189                 return LevenbergMarquardtSpace::ImproperInputParameters;
00190 
00191     /*     evaluate the function at the starting point */
00192     /*     and calculate its norm. */
00193     nfev = 1;
00194     if ( functor(x, fvec) < 0)
00195         return LevenbergMarquardtSpace::UserAsked;
00196     fnorm = fvec.stableNorm();
00197 
00198     /*     initialize levenberg-marquardt parameter and iteration counter. */
00199     par = 0.;
00200     iter = 1;
00201 
00202     return LevenbergMarquardtSpace::NotStarted;
00203 }
00204 
00205 template<typename FunctorType, typename Scalar>
00206 LevenbergMarquardtSpace::Status
00207 LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType  &x)
00208 {
00209     using std::abs;
00210     using std::sqrt;
00211     
00212     eigen_assert(x.size()==n); // check the caller is not cheating us
00213 
00214     /* calculate the jacobian matrix. */
00215     Index df_ret = functor.df(x, fjac);
00216     if (df_ret<0)
00217         return LevenbergMarquardtSpace::UserAsked;
00218     if (df_ret>0)
00219         // numerical diff, we evaluated the function df_ret times
00220         nfev += df_ret;
00221     else njev++;
00222 
00223     /* compute the qr factorization of the jacobian. */
00224     wa2 = fjac.colwise().blueNorm();
00225     ColPivHouseholderQR<JacobianType> qrfac(fjac);
00226     fjac = qrfac.matrixQR();
00227     permutation = qrfac.colsPermutation();
00228 
00229     /* on the first iteration and if external scaling is not used, scale according */
00230     /* to the norms of the columns of the initial jacobian. */
00231     if (iter == 1) {
00232         if (!useExternalScaling)
00233             for (Index j = 0; j < n; ++j)
00234                 diag[j] = (wa2[j]==0.)? 1. : wa2[j];
00235 
00236         /* on the first iteration, calculate the norm of the scaled x */
00237         /* and initialize the step bound delta. */
00238         xnorm = diag.cwiseProduct(x).stableNorm();
00239         delta = parameters.factor * xnorm;
00240         if (delta == 0.)
00241             delta = parameters.factor;
00242     }
00243 
00244     /* form (q transpose)*fvec and store the first n components in */
00245     /* qtf. */
00246     wa4 = fvec;
00247     wa4.applyOnTheLeft(qrfac.householderQ().adjoint());
00248     qtf = wa4.head(n);
00249 
00250     /* compute the norm of the scaled gradient. */
00251     gnorm = 0.;
00252     if (fnorm != 0.)
00253         for (Index j = 0; j < n; ++j)
00254             if (wa2[permutation.indices()[j]] != 0.)
00255                 gnorm = (std::max)(gnorm, abs( fjac.col(j).head(j+1).dot(qtf.head(j+1)/fnorm) / wa2[permutation.indices()[j]]));
00256 
00257     /* test for convergence of the gradient norm. */
00258     if (gnorm <= parameters.gtol)
00259         return LevenbergMarquardtSpace::CosinusTooSmall;
00260 
00261     /* rescale if necessary. */
00262     if (!useExternalScaling)
00263         diag = diag.cwiseMax(wa2);
00264 
00265     do {
00266 
00267         /* determine the levenberg-marquardt parameter. */
00268         internal::lmpar2<Scalar>(qrfac, diag, qtf, delta, par, wa1);
00269 
00270         /* store the direction p and x + p. calculate the norm of p. */
00271         wa1 = -wa1;
00272         wa2 = x + wa1;
00273         pnorm = diag.cwiseProduct(wa1).stableNorm();
00274 
00275         /* on the first iteration, adjust the initial step bound. */
00276         if (iter == 1)
00277             delta = (std::min)(delta,pnorm);
00278 
00279         /* evaluate the function at x + p and calculate its norm. */
00280         if ( functor(wa2, wa4) < 0)
00281             return LevenbergMarquardtSpace::UserAsked;
00282         ++nfev;
00283         fnorm1 = wa4.stableNorm();
00284 
00285         /* compute the scaled actual reduction. */
00286         actred = -1.;
00287         if (Scalar(.1) * fnorm1 < fnorm)
00288             actred = 1. - numext::abs2(fnorm1 / fnorm);
00289 
00290         /* compute the scaled predicted reduction and */
00291         /* the scaled directional derivative. */
00292         wa3 = fjac.template triangularView<Upper>() * (qrfac.colsPermutation().inverse() *wa1);
00293         temp1 = numext::abs2(wa3.stableNorm() / fnorm);
00294         temp2 = numext::abs2(sqrt(par) * pnorm / fnorm);
00295         prered = temp1 + temp2 / Scalar(.5);
00296         dirder = -(temp1 + temp2);
00297 
00298         /* compute the ratio of the actual to the predicted */
00299         /* reduction. */
00300         ratio = 0.;
00301         if (prered != 0.)
00302             ratio = actred / prered;
00303 
00304         /* update the step bound. */
00305         if (ratio <= Scalar(.25)) {
00306             if (actred >= 0.)
00307                 temp = Scalar(.5);
00308             if (actred < 0.)
00309                 temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
00310             if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1))
00311                 temp = Scalar(.1);
00312             /* Computing MIN */
00313             delta = temp * (std::min)(delta, pnorm / Scalar(.1));
00314             par /= temp;
00315         } else if (!(par != 0. && ratio < Scalar(.75))) {
00316             delta = pnorm / Scalar(.5);
00317             par = Scalar(.5) * par;
00318         }
00319 
00320         /* test for successful iteration. */
00321         if (ratio >= Scalar(1e-4)) {
00322             /* successful iteration. update x, fvec, and their norms. */
00323             x = wa2;
00324             wa2 = diag.cwiseProduct(x);
00325             fvec = wa4;
00326             xnorm = wa2.stableNorm();
00327             fnorm = fnorm1;
00328             ++iter;
00329         }
00330 
00331         /* tests for convergence. */
00332         if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm)
00333             return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall;
00334         if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.)
00335             return LevenbergMarquardtSpace::RelativeReductionTooSmall;
00336         if (delta <= parameters.xtol * xnorm)
00337             return LevenbergMarquardtSpace::RelativeErrorTooSmall;
00338 
00339         /* tests for termination and stringent tolerances. */
00340         if (nfev >= parameters.maxfev)
00341             return LevenbergMarquardtSpace::TooManyFunctionEvaluation;
00342         if (abs(actred) <= NumTraits<Scalar>::epsilon() && prered <= NumTraits<Scalar>::epsilon() && Scalar(.5) * ratio <= 1.)
00343             return LevenbergMarquardtSpace::FtolTooSmall;
00344         if (delta <= NumTraits<Scalar>::epsilon() * xnorm)
00345             return LevenbergMarquardtSpace::XtolTooSmall;
00346         if (gnorm <= NumTraits<Scalar>::epsilon())
00347             return LevenbergMarquardtSpace::GtolTooSmall;
00348 
00349     } while (ratio < Scalar(1e-4));
00350 
00351     return LevenbergMarquardtSpace::Running;
00352 }
00353 
00354 template<typename FunctorType, typename Scalar>
00355 LevenbergMarquardtSpace::Status
00356 LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
00357         FVectorType  &x,
00358         const Scalar tol
00359         )
00360 {
00361     n = x.size();
00362     m = functor.values();
00363 
00364     /* check the input parameters for errors. */
00365     if (n <= 0 || m < n || tol < 0.)
00366         return LevenbergMarquardtSpace::ImproperInputParameters;
00367 
00368     resetParameters();
00369     parameters.ftol = tol;
00370     parameters.xtol = tol;
00371     parameters.maxfev = 100*(n+1);
00372 
00373     return minimizeOptimumStorage(x);
00374 }
00375 
00376 template<typename FunctorType, typename Scalar>
00377 LevenbergMarquardtSpace::Status
00378 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(FVectorType  &x)
00379 {
00380     n = x.size();
00381     m = functor.values();
00382 
00383     wa1.resize(n); wa2.resize(n); wa3.resize(n);
00384     wa4.resize(m);
00385     fvec.resize(m);
00386     // Only R is stored in fjac. Q is only used to compute 'qtf', which is
00387     // Q.transpose()*rhs. qtf will be updated using givens rotation,
00388     // instead of storing them in Q.
00389     // The purpose it to only use a nxn matrix, instead of mxn here, so
00390     // that we can handle cases where m>>n :
00391     fjac.resize(n, n);
00392     if (!useExternalScaling)
00393         diag.resize(n);
00394     eigen_assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
00395     qtf.resize(n);
00396 
00397     /* Function Body */
00398     nfev = 0;
00399     njev = 0;
00400 
00401     /*     check the input parameters for errors. */
00402     if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
00403         return LevenbergMarquardtSpace::ImproperInputParameters;
00404 
00405     if (useExternalScaling)
00406         for (Index j = 0; j < n; ++j)
00407             if (diag[j] <= 0.)
00408                 return LevenbergMarquardtSpace::ImproperInputParameters;
00409 
00410     /*     evaluate the function at the starting point */
00411     /*     and calculate its norm. */
00412     nfev = 1;
00413     if ( functor(x, fvec) < 0)
00414         return LevenbergMarquardtSpace::UserAsked;
00415     fnorm = fvec.stableNorm();
00416 
00417     /*     initialize levenberg-marquardt parameter and iteration counter. */
00418     par = 0.;
00419     iter = 1;
00420 
00421     return LevenbergMarquardtSpace::NotStarted;
00422 }
00423 
00424 
00425 template<typename FunctorType, typename Scalar>
00426 LevenbergMarquardtSpace::Status
00427 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(FVectorType  &x)
00428 {
00429     using std::abs;
00430     using std::sqrt;
00431     
00432     eigen_assert(x.size()==n); // check the caller is not cheating us
00433 
00434     Index i, j;
00435     bool sing;
00436 
00437     /* compute the qr factorization of the jacobian matrix */
00438     /* calculated one row at a time, while simultaneously */
00439     /* forming (q transpose)*fvec and storing the first */
00440     /* n components in qtf. */
00441     qtf.fill(0.);
00442     fjac.fill(0.);
00443     Index rownb = 2;
00444     for (i = 0; i < m; ++i) {
00445         if (functor.df(x, wa3, rownb) < 0) return LevenbergMarquardtSpace::UserAsked;
00446         internal::rwupdt<Scalar>(fjac, wa3, qtf, fvec[i]);
00447         ++rownb;
00448     }
00449     ++njev;
00450 
00451     /* if the jacobian is rank deficient, call qrfac to */
00452     /* reorder its columns and update the components of qtf. */
00453     sing = false;
00454     for (j = 0; j < n; ++j) {
00455         if (fjac(j,j) == 0.)
00456             sing = true;
00457         wa2[j] = fjac.col(j).head(j).stableNorm();
00458     }
00459     permutation.setIdentity(n);
00460     if (sing) {
00461         wa2 = fjac.colwise().blueNorm();
00462         // TODO We have no unit test covering this code path, do not modify
00463         // until it is carefully tested
00464         ColPivHouseholderQR<JacobianType> qrfac(fjac);
00465         fjac = qrfac.matrixQR();
00466         wa1 = fjac.diagonal();
00467         fjac.diagonal() = qrfac.hCoeffs();
00468         permutation = qrfac.colsPermutation();
00469         // TODO : avoid this:
00470         for(Index ii=0; ii< fjac.cols(); ii++) fjac.col(ii).segment(ii+1, fjac.rows()-ii-1) *= fjac(ii,ii); // rescale vectors
00471 
00472         for (j = 0; j < n; ++j) {
00473             if (fjac(j,j) != 0.) {
00474                 sum = 0.;
00475                 for (i = j; i < n; ++i)
00476                     sum += fjac(i,j) * qtf[i];
00477                 temp = -sum / fjac(j,j);
00478                 for (i = j; i < n; ++i)
00479                     qtf[i] += fjac(i,j) * temp;
00480             }
00481             fjac(j,j) = wa1[j];
00482         }
00483     }
00484 
00485     /* on the first iteration and if external scaling is not used, scale according */
00486     /* to the norms of the columns of the initial jacobian. */
00487     if (iter == 1) {
00488         if (!useExternalScaling)
00489             for (j = 0; j < n; ++j)
00490                 diag[j] = (wa2[j]==0.)? 1. : wa2[j];
00491 
00492         /* on the first iteration, calculate the norm of the scaled x */
00493         /* and initialize the step bound delta. */
00494         xnorm = diag.cwiseProduct(x).stableNorm();
00495         delta = parameters.factor * xnorm;
00496         if (delta == 0.)
00497             delta = parameters.factor;
00498     }
00499 
00500     /* compute the norm of the scaled gradient. */
00501     gnorm = 0.;
00502     if (fnorm != 0.)
00503         for (j = 0; j < n; ++j)
00504             if (wa2[permutation.indices()[j]] != 0.)
00505                 gnorm = (std::max)(gnorm, abs( fjac.col(j).head(j+1).dot(qtf.head(j+1)/fnorm) / wa2[permutation.indices()[j]]));
00506 
00507     /* test for convergence of the gradient norm. */
00508     if (gnorm <= parameters.gtol)
00509         return LevenbergMarquardtSpace::CosinusTooSmall;
00510 
00511     /* rescale if necessary. */
00512     if (!useExternalScaling)
00513         diag = diag.cwiseMax(wa2);
00514 
00515     do {
00516 
00517         /* determine the levenberg-marquardt parameter. */
00518         internal::lmpar<Scalar>(fjac, permutation.indices(), diag, qtf, delta, par, wa1);
00519 
00520         /* store the direction p and x + p. calculate the norm of p. */
00521         wa1 = -wa1;
00522         wa2 = x + wa1;
00523         pnorm = diag.cwiseProduct(wa1).stableNorm();
00524 
00525         /* on the first iteration, adjust the initial step bound. */
00526         if (iter == 1)
00527             delta = (std::min)(delta,pnorm);
00528 
00529         /* evaluate the function at x + p and calculate its norm. */
00530         if ( functor(wa2, wa4) < 0)
00531             return LevenbergMarquardtSpace::UserAsked;
00532         ++nfev;
00533         fnorm1 = wa4.stableNorm();
00534 
00535         /* compute the scaled actual reduction. */
00536         actred = -1.;
00537         if (Scalar(.1) * fnorm1 < fnorm)
00538             actred = 1. - numext::abs2(fnorm1 / fnorm);
00539 
00540         /* compute the scaled predicted reduction and */
00541         /* the scaled directional derivative. */
00542         wa3 = fjac.topLeftCorner(n,n).template triangularView<Upper>() * (permutation.inverse() * wa1);
00543         temp1 = numext::abs2(wa3.stableNorm() / fnorm);
00544         temp2 = numext::abs2(sqrt(par) * pnorm / fnorm);
00545         prered = temp1 + temp2 / Scalar(.5);
00546         dirder = -(temp1 + temp2);
00547 
00548         /* compute the ratio of the actual to the predicted */
00549         /* reduction. */
00550         ratio = 0.;
00551         if (prered != 0.)
00552             ratio = actred / prered;
00553 
00554         /* update the step bound. */
00555         if (ratio <= Scalar(.25)) {
00556             if (actred >= 0.)
00557                 temp = Scalar(.5);
00558             if (actred < 0.)
00559                 temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
00560             if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1))
00561                 temp = Scalar(.1);
00562             /* Computing MIN */
00563             delta = temp * (std::min)(delta, pnorm / Scalar(.1));
00564             par /= temp;
00565         } else if (!(par != 0. && ratio < Scalar(.75))) {
00566             delta = pnorm / Scalar(.5);
00567             par = Scalar(.5) * par;
00568         }
00569 
00570         /* test for successful iteration. */
00571         if (ratio >= Scalar(1e-4)) {
00572             /* successful iteration. update x, fvec, and their norms. */
00573             x = wa2;
00574             wa2 = diag.cwiseProduct(x);
00575             fvec = wa4;
00576             xnorm = wa2.stableNorm();
00577             fnorm = fnorm1;
00578             ++iter;
00579         }
00580 
00581         /* tests for convergence. */
00582         if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm)
00583             return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall;
00584         if (abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.)
00585             return LevenbergMarquardtSpace::RelativeReductionTooSmall;
00586         if (delta <= parameters.xtol * xnorm)
00587             return LevenbergMarquardtSpace::RelativeErrorTooSmall;
00588 
00589         /* tests for termination and stringent tolerances. */
00590         if (nfev >= parameters.maxfev)
00591             return LevenbergMarquardtSpace::TooManyFunctionEvaluation;
00592         if (abs(actred) <= NumTraits<Scalar>::epsilon() && prered <= NumTraits<Scalar>::epsilon() && Scalar(.5) * ratio <= 1.)
00593             return LevenbergMarquardtSpace::FtolTooSmall;
00594         if (delta <= NumTraits<Scalar>::epsilon() * xnorm)
00595             return LevenbergMarquardtSpace::XtolTooSmall;
00596         if (gnorm <= NumTraits<Scalar>::epsilon())
00597             return LevenbergMarquardtSpace::GtolTooSmall;
00598 
00599     } while (ratio < Scalar(1e-4));
00600 
00601     return LevenbergMarquardtSpace::Running;
00602 }
00603 
00604 template<typename FunctorType, typename Scalar>
00605 LevenbergMarquardtSpace::Status
00606 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(FVectorType  &x)
00607 {
00608     LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x);
00609     if (status==LevenbergMarquardtSpace::ImproperInputParameters)
00610         return status;
00611     do {
00612         status = minimizeOptimumStorageOneStep(x);
00613     } while (status==LevenbergMarquardtSpace::Running);
00614     return status;
00615 }
00616 
00617 template<typename FunctorType, typename Scalar>
00618 LevenbergMarquardtSpace::Status
00619 LevenbergMarquardt<FunctorType,Scalar>::lmdif1(
00620         FunctorType &functor,
00621         FVectorType  &x,
00622         Index *nfev,
00623         const Scalar tol
00624         )
00625 {
00626     Index n = x.size();
00627     Index m = functor.values();
00628 
00629     /* check the input parameters for errors. */
00630     if (n <= 0 || m < n || tol < 0.)
00631         return LevenbergMarquardtSpace::ImproperInputParameters;
00632 
00633     NumericalDiff<FunctorType> numDiff(functor);
00634     // embedded LevenbergMarquardt
00635     LevenbergMarquardt<NumericalDiff<FunctorType>, Scalar > lm(numDiff);
00636     lm.parameters.ftol = tol;
00637     lm.parameters.xtol = tol;
00638     lm.parameters.maxfev = 200*(n+1);
00639 
00640     LevenbergMarquardtSpace::Status info = LevenbergMarquardtSpace::Status(lm.minimize(x));
00641     if (nfev)
00642         * nfev = lm.nfev;
00643     return info;
00644 }
00645 
00646 } // end namespace Eigen
00647 
00648 #endif // EIGEN_LEVENBERGMARQUARDT__H
00649 
00650 //vim: ai ts=4 sts=4 et sw=4


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:37:51