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