$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org> 00005 // Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> 00006 // 00007 // The algorithm of this class initially comes from MINPACK whose original authors are: 00008 // Copyright Jorge More - Argonne National Laboratory 00009 // Copyright Burt Garbow - Argonne National Laboratory 00010 // Copyright Ken Hillstrom - Argonne National Laboratory 00011 // 00012 // This Source Code Form is subject to the terms of the Minpack license 00013 // (a BSD-like license) described in the campaigned CopyrightMINPACK.txt file. 00014 // 00015 // This Source Code Form is subject to the terms of the Mozilla 00016 // Public License v. 2.0. If a copy of the MPL was not distributed 00017 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00018 00019 #ifndef EIGEN_LEVENBERGMARQUARDT_H 00020 #define EIGEN_LEVENBERGMARQUARDT_H 00021 00022 00023 namespace Eigen { 00024 namespace LevenbergMarquardtSpace { 00025 enum Status { 00026 NotStarted = -2, 00027 Running = -1, 00028 ImproperInputParameters = 0, 00029 RelativeReductionTooSmall = 1, 00030 RelativeErrorTooSmall = 2, 00031 RelativeErrorAndReductionTooSmall = 3, 00032 CosinusTooSmall = 4, 00033 TooManyFunctionEvaluation = 5, 00034 FtolTooSmall = 6, 00035 XtolTooSmall = 7, 00036 GtolTooSmall = 8, 00037 UserAsked = 9 00038 }; 00039 } 00040 00041 template <typename _Scalar, int NX=Dynamic, int NY=Dynamic> 00042 struct DenseFunctor 00043 { 00044 typedef _Scalar Scalar; 00045 enum { 00046 InputsAtCompileTime = NX, 00047 ValuesAtCompileTime = NY 00048 }; 00049 typedef Matrix<Scalar,InputsAtCompileTime,1> InputType; 00050 typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType; 00051 typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType; 00052 typedef ColPivHouseholderQR<JacobianType> QRSolver; 00053 const int m_inputs, m_values; 00054 00055 DenseFunctor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {} 00056 DenseFunctor(int inputs, int values) : m_inputs(inputs), m_values(values) {} 00057 00058 int inputs() const { return m_inputs; } 00059 int values() const { return m_values; } 00060 00061 //int operator()(const InputType &x, ValueType& fvec) { } 00062 // should be defined in derived classes 00063 00064 //int df(const InputType &x, JacobianType& fjac) { } 00065 // should be defined in derived classes 00066 }; 00067 00068 template <typename _Scalar, typename _Index> 00069 struct SparseFunctor 00070 { 00071 typedef _Scalar Scalar; 00072 typedef _Index Index; 00073 typedef Matrix<Scalar,Dynamic,1> InputType; 00074 typedef Matrix<Scalar,Dynamic,1> ValueType; 00075 typedef SparseMatrix<Scalar, ColMajor, Index> JacobianType; 00076 typedef SparseQR<JacobianType, COLAMDOrdering<int> > QRSolver; 00077 enum { 00078 InputsAtCompileTime = Dynamic, 00079 ValuesAtCompileTime = Dynamic 00080 }; 00081 00082 SparseFunctor(int inputs, int values) : m_inputs(inputs), m_values(values) {} 00083 00084 int inputs() const { return m_inputs; } 00085 int values() const { return m_values; } 00086 00087 const int m_inputs, m_values; 00088 //int operator()(const InputType &x, ValueType& fvec) { } 00089 // to be defined in the functor 00090 00091 //int df(const InputType &x, JacobianType& fjac) { } 00092 // to be defined in the functor if no automatic differentiation 00093 00094 }; 00095 namespace internal { 00096 template <typename QRSolver, typename VectorType> 00097 void lmpar2(const QRSolver &qr, const VectorType &diag, const VectorType &qtb, 00098 typename VectorType::Scalar m_delta, typename VectorType::Scalar &par, 00099 VectorType &x); 00100 } 00109 template<typename _FunctorType> 00110 class LevenbergMarquardt : internal::no_assignment_operator 00111 { 00112 public: 00113 typedef _FunctorType FunctorType; 00114 typedef typename FunctorType::QRSolver QRSolver; 00115 typedef typename FunctorType::JacobianType JacobianType; 00116 typedef typename JacobianType::Scalar Scalar; 00117 typedef typename JacobianType::RealScalar RealScalar; 00118 typedef typename JacobianType::Index Index; 00119 typedef typename QRSolver::Index PermIndex; 00120 typedef Matrix<Scalar,Dynamic,1> FVectorType; 00121 typedef PermutationMatrix<Dynamic,Dynamic> PermutationType; 00122 public: 00123 LevenbergMarquardt(FunctorType& functor) 00124 : m_functor(functor),m_nfev(0),m_njev(0),m_fnorm(0.0),m_gnorm(0), 00125 m_isInitialized(false),m_info(InvalidInput) 00126 { 00127 resetParameters(); 00128 m_useExternalScaling=false; 00129 } 00130 00131 LevenbergMarquardtSpace::Status minimize(FVectorType &x); 00132 LevenbergMarquardtSpace::Status minimizeInit(FVectorType &x); 00133 LevenbergMarquardtSpace::Status minimizeOneStep(FVectorType &x); 00134 LevenbergMarquardtSpace::Status lmder1( 00135 FVectorType &x, 00136 const Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon()) 00137 ); 00138 static LevenbergMarquardtSpace::Status lmdif1( 00139 FunctorType &functor, 00140 FVectorType &x, 00141 Index *nfev, 00142 const Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon()) 00143 ); 00144 00146 void resetParameters() 00147 { 00148 m_factor = 100.; 00149 m_maxfev = 400; 00150 m_ftol = std::sqrt(NumTraits<RealScalar>::epsilon()); 00151 m_xtol = std::sqrt(NumTraits<RealScalar>::epsilon()); 00152 m_gtol = 0. ; 00153 m_epsfcn = 0. ; 00154 } 00155 00157 void setXtol(RealScalar xtol) { m_xtol = xtol; } 00158 00160 void setFtol(RealScalar ftol) { m_ftol = ftol; } 00161 00163 void setGtol(RealScalar gtol) { m_gtol = gtol; } 00164 00166 void setFactor(RealScalar factor) { m_factor = factor; } 00167 00169 void setEpsilon (RealScalar epsfcn) { m_epsfcn = epsfcn; } 00170 00172 void setMaxfev(Index maxfev) {m_maxfev = maxfev; } 00173 00175 void setExternalScaling(bool value) {m_useExternalScaling = value; } 00176 00178 FVectorType& diag() {return m_diag; } 00179 00181 Index iterations() { return m_iter; } 00182 00184 Index nfev() { return m_nfev; } 00185 00187 Index njev() { return m_njev; } 00188 00190 RealScalar fnorm() {return m_fnorm; } 00191 00193 RealScalar gnorm() {return m_gnorm; } 00194 00196 RealScalar lm_param(void) { return m_par; } 00197 00200 FVectorType& fvec() {return m_fvec; } 00201 00204 JacobianType& jacobian() {return m_fjac; } 00205 00209 JacobianType& matrixR() {return m_rfactor; } 00210 00213 PermutationType permutation() {return m_permutation; } 00214 00224 ComputationInfo info() const 00225 { 00226 00227 return m_info; 00228 } 00229 private: 00230 JacobianType m_fjac; 00231 JacobianType m_rfactor; // The triangular matrix R from the QR of the jacobian matrix m_fjac 00232 FunctorType &m_functor; 00233 FVectorType m_fvec, m_qtf, m_diag; 00234 Index n; 00235 Index m; 00236 Index m_nfev; 00237 Index m_njev; 00238 RealScalar m_fnorm; // Norm of the current vector function 00239 RealScalar m_gnorm; //Norm of the gradient of the error 00240 RealScalar m_factor; // 00241 Index m_maxfev; // Maximum number of function evaluation 00242 RealScalar m_ftol; //Tolerance in the norm of the vector function 00243 RealScalar m_xtol; // 00244 RealScalar m_gtol; //tolerance of the norm of the error gradient 00245 RealScalar m_epsfcn; // 00246 Index m_iter; // Number of iterations performed 00247 RealScalar m_delta; 00248 bool m_useExternalScaling; 00249 PermutationType m_permutation; 00250 FVectorType m_wa1, m_wa2, m_wa3, m_wa4; //Temporary vectors 00251 RealScalar m_par; 00252 bool m_isInitialized; // Check whether the minimization step has been called 00253 ComputationInfo m_info; 00254 }; 00255 00256 template<typename FunctorType> 00257 LevenbergMarquardtSpace::Status 00258 LevenbergMarquardt<FunctorType>::minimize(FVectorType &x) 00259 { 00260 LevenbergMarquardtSpace::Status status = minimizeInit(x); 00261 if (status==LevenbergMarquardtSpace::ImproperInputParameters) { 00262 m_isInitialized = true; 00263 return status; 00264 } 00265 do { 00266 // std::cout << " uv " << x.transpose() << "\n"; 00267 status = minimizeOneStep(x); 00268 } while (status==LevenbergMarquardtSpace::Running); 00269 m_isInitialized = true; 00270 return status; 00271 } 00272 00273 template<typename FunctorType> 00274 LevenbergMarquardtSpace::Status 00275 LevenbergMarquardt<FunctorType>::minimizeInit(FVectorType &x) 00276 { 00277 n = x.size(); 00278 m = m_functor.values(); 00279 00280 m_wa1.resize(n); m_wa2.resize(n); m_wa3.resize(n); 00281 m_wa4.resize(m); 00282 m_fvec.resize(m); 00283 //FIXME Sparse Case : Allocate space for the jacobian 00284 m_fjac.resize(m, n); 00285 // m_fjac.reserve(VectorXi::Constant(n,5)); // FIXME Find a better alternative 00286 if (!m_useExternalScaling) 00287 m_diag.resize(n); 00288 eigen_assert( (!m_useExternalScaling || m_diag.size()==n) || "When m_useExternalScaling is set, the caller must provide a valid 'm_diag'"); 00289 m_qtf.resize(n); 00290 00291 /* Function Body */ 00292 m_nfev = 0; 00293 m_njev = 0; 00294 00295 /* check the input parameters for errors. */ 00296 if (n <= 0 || m < n || m_ftol < 0. || m_xtol < 0. || m_gtol < 0. || m_maxfev <= 0 || m_factor <= 0.){ 00297 m_info = InvalidInput; 00298 return LevenbergMarquardtSpace::ImproperInputParameters; 00299 } 00300 00301 if (m_useExternalScaling) 00302 for (Index j = 0; j < n; ++j) 00303 if (m_diag[j] <= 0.) 00304 { 00305 m_info = InvalidInput; 00306 return LevenbergMarquardtSpace::ImproperInputParameters; 00307 } 00308 00309 /* evaluate the function at the starting point */ 00310 /* and calculate its norm. */ 00311 m_nfev = 1; 00312 if ( m_functor(x, m_fvec) < 0) 00313 return LevenbergMarquardtSpace::UserAsked; 00314 m_fnorm = m_fvec.stableNorm(); 00315 00316 /* initialize levenberg-marquardt parameter and iteration counter. */ 00317 m_par = 0.; 00318 m_iter = 1; 00319 00320 return LevenbergMarquardtSpace::NotStarted; 00321 } 00322 00323 template<typename FunctorType> 00324 LevenbergMarquardtSpace::Status 00325 LevenbergMarquardt<FunctorType>::lmder1( 00326 FVectorType &x, 00327 const Scalar tol 00328 ) 00329 { 00330 n = x.size(); 00331 m = m_functor.values(); 00332 00333 /* check the input parameters for errors. */ 00334 if (n <= 0 || m < n || tol < 0.) 00335 return LevenbergMarquardtSpace::ImproperInputParameters; 00336 00337 resetParameters(); 00338 m_ftol = tol; 00339 m_xtol = tol; 00340 m_maxfev = 100*(n+1); 00341 00342 return minimize(x); 00343 } 00344 00345 00346 template<typename FunctorType> 00347 LevenbergMarquardtSpace::Status 00348 LevenbergMarquardt<FunctorType>::lmdif1( 00349 FunctorType &functor, 00350 FVectorType &x, 00351 Index *nfev, 00352 const Scalar tol 00353 ) 00354 { 00355 Index n = x.size(); 00356 Index m = functor.values(); 00357 00358 /* check the input parameters for errors. */ 00359 if (n <= 0 || m < n || tol < 0.) 00360 return LevenbergMarquardtSpace::ImproperInputParameters; 00361 00362 NumericalDiff<FunctorType> numDiff(functor); 00363 // embedded LevenbergMarquardt 00364 LevenbergMarquardt<NumericalDiff<FunctorType> > lm(numDiff); 00365 lm.setFtol(tol); 00366 lm.setXtol(tol); 00367 lm.setMaxfev(200*(n+1)); 00368 00369 LevenbergMarquardtSpace::Status info = LevenbergMarquardtSpace::Status(lm.minimize(x)); 00370 if (nfev) 00371 * nfev = lm.nfev(); 00372 return info; 00373 } 00374 00375 } // end namespace Eigen 00376 00377 #endif // EIGEN_LEVENBERGMARQUARDT_H