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