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(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; } 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 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 205 template<typename FunctorType, typename Scalar> 206 LevenbergMarquardtSpace::Status 207 LevenbergMarquardt<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 354 template<typename FunctorType, typename Scalar> 355 LevenbergMarquardtSpace::Status 356 LevenbergMarquardt<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 376 template<typename FunctorType, typename Scalar> 377 LevenbergMarquardtSpace::Status 378 LevenbergMarquardt<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 425 template<typename FunctorType, typename Scalar> 426 LevenbergMarquardtSpace::Status 427 LevenbergMarquardt<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 604 template<typename FunctorType, typename Scalar> 605 LevenbergMarquardtSpace::Status 606 LevenbergMarquardt<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 617 template<typename FunctorType, typename Scalar> 618 LevenbergMarquardtSpace::Status 619 LevenbergMarquardt<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 651