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