source: pacpussensors/trunk/Vislab/lib3dv/eigen/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

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