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