diff options
author | Desire NUENTSA <desire.nuentsa_wakam@inria.fr> | 2013-02-20 14:10:14 +0100 |
---|---|---|
committer | Desire NUENTSA <desire.nuentsa_wakam@inria.fr> | 2013-02-20 14:10:14 +0100 |
commit | febf8e5d7b69130e67025c79bf5c45500df570f4 (patch) | |
tree | f753f4c166c54621bfa2fad554c7062c2068a60d /unsupported/Eigen/src/LevenbergMarquardt | |
parent | dca7190e15d5f787cfe84830bc4169e013154c5b (diff) |
Set built-in sparse QR as the default sparse solver and add ComputationInfo for Levenberg Marquardt,
Diffstat (limited to 'unsupported/Eigen/src/LevenbergMarquardt')
-rw-r--r-- | unsupported/Eigen/src/LevenbergMarquardt/LMonestep.h | 66 | ||||
-rw-r--r-- | unsupported/Eigen/src/LevenbergMarquardt/LMpar.h | 18 | ||||
-rw-r--r-- | unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h | 77 |
3 files changed, 109 insertions, 52 deletions
diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LMonestep.h b/unsupported/Eigen/src/LevenbergMarquardt/LMonestep.h index 51ade866c..351c28c2c 100644 --- a/unsupported/Eigen/src/LevenbergMarquardt/LMonestep.h +++ b/unsupported/Eigen/src/LevenbergMarquardt/LMonestep.h @@ -40,13 +40,14 @@ LevenbergMarquardt<FunctorType>::minimizeOneStep(FVectorType &x) /* compute the qr factorization of the jacobian. */ for (int j = 0; j < x.size(); ++j) - m_wa2(j) = m_fjac.col(j).norm(); - //FIXME Implement bluenorm for sparse vectors -// m_wa2 = m_fjac.colwise().blueNorm(); - QRSolver qrfac(m_fjac); //FIXME Check if the QR decomposition succeed + m_wa2(j) = m_fjac.col(j).blueNorm(); + QRSolver qrfac(m_fjac); + if(qrfac.info() != Success) { + m_info = NumericalIssue; + return LevenbergMarquardtSpace::ImproperInputParameters; + } // Make a copy of the first factor with the associated permutation - JacobianType rfactor; - rfactor = qrfac.matrixQR(); + m_rfactor = qrfac.matrixR(); m_permutation = (qrfac.colsPermutation()); /* on the first iteration and if external scaling is not used, scale according */ @@ -75,11 +76,13 @@ LevenbergMarquardt<FunctorType>::minimizeOneStep(FVectorType &x) if (m_fnorm != 0.) for (Index j = 0; j < n; ++j) if (m_wa2[m_permutation.indices()[j]] != 0.) - m_gnorm = (std::max)(m_gnorm, abs( rfactor.col(j).head(j+1).dot(m_qtf.head(j+1)/m_fnorm) / m_wa2[m_permutation.indices()[j]])); + m_gnorm = (std::max)(m_gnorm, abs( m_rfactor.col(j).head(j+1).dot(m_qtf.head(j+1)/m_fnorm) / m_wa2[m_permutation.indices()[j]])); /* test for convergence of the gradient norm. */ - if (m_gnorm <= m_gtol) - return LevenbergMarquardtSpace::CosinusTooSmall; + if (m_gnorm <= m_gtol) { + m_info = Success; + return LevenbergMarquardtSpace::CosinusTooSmall; + } /* rescale if necessary. */ if (!m_useExternalScaling) @@ -111,7 +114,7 @@ LevenbergMarquardt<FunctorType>::minimizeOneStep(FVectorType &x) /* compute the scaled predicted reduction and */ /* the scaled directional derivative. */ - m_wa3 = rfactor.template triangularView<Upper>() * (m_permutation.inverse() *m_wa1); + m_wa3 = m_rfactor.template triangularView<Upper>() * (m_permutation.inverse() *m_wa1); temp1 = internal::abs2(m_wa3.stableNorm() / m_fnorm); temp2 = internal::abs2(sqrt(m_par) * pnorm / m_fnorm); prered = temp1 + temp2 / Scalar(.5); @@ -152,21 +155,42 @@ LevenbergMarquardt<FunctorType>::minimizeOneStep(FVectorType &x) /* tests for convergence. */ if (abs(actred) <= m_ftol && prered <= m_ftol && Scalar(.5) * ratio <= 1. && m_delta <= m_xtol * xnorm) - return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall; - if (abs(actred) <= m_ftol && prered <= m_ftol && Scalar(.5) * ratio <= 1.) - return LevenbergMarquardtSpace::RelativeReductionTooSmall; + { + m_info = Success; + return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall; + } + if (abs(actred) <= m_ftol && prered <= m_ftol && Scalar(.5) * ratio <= 1.) + { + m_info = Success; + return LevenbergMarquardtSpace::RelativeReductionTooSmall; + } if (m_delta <= m_xtol * xnorm) - return LevenbergMarquardtSpace::RelativeErrorTooSmall; + { + m_info = Success; + return LevenbergMarquardtSpace::RelativeErrorTooSmall; + } /* tests for termination and stringent tolerances. */ - if (m_nfev >= m_maxfev) - return LevenbergMarquardtSpace::TooManyFunctionEvaluation; + if (m_nfev >= m_maxfev) + { + m_info = NoConvergence; + return LevenbergMarquardtSpace::TooManyFunctionEvaluation; + } if (abs(actred) <= NumTraits<Scalar>::epsilon() && prered <= NumTraits<Scalar>::epsilon() && Scalar(.5) * ratio <= 1.) - return LevenbergMarquardtSpace::FtolTooSmall; - if (m_delta <= NumTraits<Scalar>::epsilon() * xnorm) - return LevenbergMarquardtSpace::XtolTooSmall; + { + m_info = Success; + return LevenbergMarquardtSpace::FtolTooSmall; + } + if (m_delta <= NumTraits<Scalar>::epsilon() * xnorm) + { + m_info = Success; + return LevenbergMarquardtSpace::XtolTooSmall; + } if (m_gnorm <= NumTraits<Scalar>::epsilon()) - return LevenbergMarquardtSpace::GtolTooSmall; + { + m_info = Success; + return LevenbergMarquardtSpace::GtolTooSmall; + } } while (ratio < Scalar(1e-4)); @@ -176,4 +200,4 @@ LevenbergMarquardt<FunctorType>::minimizeOneStep(FVectorType &x) } // end namespace Eigen -#endif // EIGEN_LMONESTEP_H +#endif // EIGEN_LMONESTEP_H
\ No newline at end of file diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h b/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h index af76a9799..9532042d9 100644 --- a/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h +++ b/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h @@ -40,11 +40,15 @@ namespace internal { Scalar temp, paru; Scalar gnorm; Scalar dxnorm; - + + // Make a copy of the triangular factor. + // This copy is modified during call the qrsolv + MatrixType s; + s = qr.matrixR(); /* Function Body */ const Scalar dwarf = (std::numeric_limits<Scalar>::min)(); - const Index n = qr.matrixQR().cols(); + const Index n = qr.matrixR().cols(); eigen_assert(n==diag.size()); eigen_assert(n==qtb.size()); @@ -58,8 +62,7 @@ namespace internal { wa1 = qtb; wa1.tail(n-rank).setZero(); //FIXME There is no solve in place for sparse triangularView - //qr.matrixQR().topLeftCorner(rank, rank).template triangularView<Upper>().solveInPlace(wa1.head(rank)); - wa1.head(rank) = qr.matrixQR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(qtb.head(rank)); + wa1.head(rank) = s.topLeftCorner(rank,rank).template triangularView<Upper>().solve(qtb.head(rank)); x = qr.colsPermutation()*wa1; @@ -81,14 +84,14 @@ namespace internal { parl = 0.; if (rank==n) { wa1 = qr.colsPermutation().inverse() * diag.cwiseProduct(wa2)/dxnorm; - qr.matrixQR().topLeftCorner(n, n).transpose().template triangularView<Lower>().solveInPlace(wa1); + s.topLeftCorner(n,n).transpose().template triangularView<Lower>().solveInPlace(wa1); temp = wa1.blueNorm(); parl = fp / m_delta / temp / temp; } /* calculate an upper bound, paru, for the zero of the function. */ for (j = 0; j < n; ++j) - wa1[j] = qr.matrixQR().col(j).head(j+1).dot(qtb.head(j+1)) / diag[qr.colsPermutation().indices()(j)]; + wa1[j] = s.col(j).head(j+1).dot(qtb.head(j+1)) / diag[qr.colsPermutation().indices()(j)]; gnorm = wa1.stableNorm(); paru = gnorm / m_delta; @@ -103,8 +106,6 @@ namespace internal { par = gnorm / dxnorm; /* beginning of an iteration. */ - MatrixType s; - s = qr.matrixQR(); while (true) { ++iter; @@ -130,7 +131,6 @@ namespace internal { /* compute the newton correction. */ wa1 = qr.colsPermutation().inverse() * diag.cwiseProduct(wa2/dxnorm); // we could almost use this here, but the diagonal is outside qr, in sdiag[] - // qr.matrixQR().topLeftCorner(n, n).transpose().template triangularView<Lower>().solveInPlace(wa1); for (j = 0; j < n; ++j) { wa1[j] /= sdiag[j]; temp = wa1[j]; diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h index bf2423ef0..8d5538d69 100644 --- a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h +++ b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h @@ -65,7 +65,6 @@ struct DenseFunctor // should be defined in derived classes }; -#ifdef EIGEN_SPQR_SUPPORT template <typename _Scalar, typename _Index> struct SparseFunctor { @@ -74,7 +73,11 @@ struct SparseFunctor typedef Matrix<Scalar,Dynamic,1> InputType; typedef Matrix<Scalar,Dynamic,1> ValueType; typedef SparseMatrix<Scalar, ColMajor, Index> JacobianType; - typedef SPQR<JacobianType> QRSolver; + typedef SparseQR<JacobianType, COLAMDOrdering<int> > QRSolver; + enum { + InputsAtCompileTime = Dynamic, + ValuesAtCompileTime = Dynamic + }; SparseFunctor(int inputs, int values) : m_inputs(inputs), m_values(values) {} @@ -89,7 +92,6 @@ struct SparseFunctor // to be defined in the functor if no automatic differentiation }; -#endif namespace internal { template <typename QRSolver, typename VectorType> void lmpar2(const QRSolver &qr, const VectorType &diag, const VectorType &qtb, @@ -119,7 +121,8 @@ class LevenbergMarquardt typedef PermutationMatrix<Dynamic,Dynamic> PermutationType; public: LevenbergMarquardt(FunctorType& functor) - : m_functor(functor),m_nfev(0),m_njev(0),m_fnorm(0.0),m_gnorm(0) + : m_functor(functor),m_nfev(0),m_njev(0),m_fnorm(0.0),m_gnorm(0), + m_isInitialized(false),m_info(InvalidInput) { resetParameters(); m_useExternalScaling=false; @@ -171,41 +174,61 @@ class LevenbergMarquardt /** Use an external Scaling. If set to true, pass a nonzero diagonal to diag() */ void setExternalScaling(bool value) {m_useExternalScaling = value; } - /** Get a reference to the diagonal of the jacobian */ + /** \returns a reference to the diagonal of the jacobian */ FVectorType& diag() {return m_diag; } - /** Number of iterations performed */ + /** \returns the number of iterations performed */ Index iterations() { return m_iter; } - /** Number of functions evaluation */ + /** \returns the number of functions evaluation */ Index nfev() { return m_nfev; } - /** Number of jacobian evaluation */ + /** \returns the number of jacobian evaluation */ Index njev() { return m_njev; } - /** Norm of current vector function */ + /** \returns the norm of current vector function */ RealScalar fnorm() {return m_fnorm; } - /** Norm of the gradient of the error */ + /** \returns the norm of the gradient of the error */ RealScalar gnorm() {return m_gnorm; } - /** the LevenbergMarquardt parameter */ + /** \returns the LevenbergMarquardt parameter */ RealScalar lm_param(void) { return m_par; } - /** reference to the current vector function + /** \returns a reference to the current vector function */ FVectorType& fvec() {return m_fvec; } - /** reference to the matrix where the current Jacobian matrix is stored + /** \returns a reference to the matrix where the current Jacobian matrix is stored */ - JacobianType& fjac() {return m_fjac; } + JacobianType& jacobian() {return m_fjac; } - /** the permutation used + /** \returns a reference to the triangular matrix R from the QR of the jacobian matrix. + * \sa jacobian() + */ + JacobianType& matrixR() {return m_rfactor; } + + /** the permutation used in the QR factorization */ PermutationType permutation() {return m_permutation; } + /** + * \brief Reports whether the minimization was successful + * \returns \c Success if the minimization was succesful, + * \c NumericalIssue if a numerical problem arises during the + * minimization process, for exemple during the QR factorization + * \c NoConvergence if the minimization did not converge after + * the maximum number of function evaluation allowed + * \c InvalidInput if the input matrix is invalid + */ + ComputationInfo info() const + { + + return m_info; + } private: JacobianType m_fjac; + JacobianType m_rfactor; // The triangular matrix R from the QR of the jacobian matrix m_fjac FunctorType &m_functor; FVectorType m_fvec, m_qtf, m_diag; Index n; @@ -226,6 +249,8 @@ class LevenbergMarquardt PermutationType m_permutation; FVectorType m_wa1, m_wa2, m_wa3, m_wa4; //Temporary vectors RealScalar m_par; + bool m_isInitialized; // Check whether the minimization step has been called + ComputationInfo m_info; }; template<typename FunctorType> @@ -233,13 +258,16 @@ LevenbergMarquardtSpace::Status LevenbergMarquardt<FunctorType>::minimize(FVectorType &x) { LevenbergMarquardtSpace::Status status = minimizeInit(x); - if (status==LevenbergMarquardtSpace::ImproperInputParameters) - return status; + if (status==LevenbergMarquardtSpace::ImproperInputParameters) { + m_isInitialized = true; + return status; + } do { // std::cout << " uv " << x.transpose() << "\n"; status = minimizeOneStep(x); } while (status==LevenbergMarquardtSpace::Running); - return status; + m_isInitialized = true; + return status; } template<typename FunctorType> @@ -265,13 +293,18 @@ LevenbergMarquardt<FunctorType>::minimizeInit(FVectorType &x) m_njev = 0; /* check the input parameters for errors. */ - if (n <= 0 || m < n || m_ftol < 0. || m_xtol < 0. || m_gtol < 0. || m_maxfev <= 0 || m_factor <= 0.) - return LevenbergMarquardtSpace::ImproperInputParameters; + if (n <= 0 || m < n || m_ftol < 0. || m_xtol < 0. || m_gtol < 0. || m_maxfev <= 0 || m_factor <= 0.){ + m_info = InvalidInput; + return LevenbergMarquardtSpace::ImproperInputParameters; + } if (m_useExternalScaling) for (Index j = 0; j < n; ++j) - if (m_diag[j] <= 0.) - return LevenbergMarquardtSpace::ImproperInputParameters; + if (m_diag[j] <= 0.) + { + return LevenbergMarquardtSpace::ImproperInputParameters; + m_info = InvalidInput; + } /* evaluate the function at the starting point */ /* and calculate its norm. */ |