GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-10 18:59:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

Set built-in sparse QR as the default sparse solver and add ComputationInfo for Levenberg Marquardt,

Browse Source
This commit is contained in:
Desire NUENTSA 2013-02-20 14:10:14 +01:00
parent dca7190e15
commit febf8e5d7b
5 changed files with 122 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
unsupported/Eigen/LevenbergMarquardt
View File

@ -16,9 +16,8 @@
#include <Eigen/Jacobi>
#include <Eigen/QR>
#include <unsupported/Eigen/NumericalDiff>
#ifdef EIGEN_SPQR_SUPPORT
#include <Eigen/SPQRSupport>
#endif
#include <Eigen/SparseQR>
/**
* \defgroup LevenbergMarquardt_Module Levenberg-Marquardt module

66
unsupported/Eigen/src/LevenbergMarquardt/LMonestep.h
View File

@ -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

18
unsupported/Eigen/src/LevenbergMarquardt/LMpar.h
View File

@ -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];

77
unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h
View File

@ -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. */

22
unsupported/test/levenberg_marquardt.cpp
View File

@ -12,7 +12,7 @@
#include <stdio.h>
#include "main.h"
#include <Eigen/LevenbergMarquardt>
#include <unsupported/Eigen/LevenbergMarquardt>
// This disables some useless Warnings on MSVC.
// It is intended to be done for this test only.
@ -115,7 +115,7 @@ void testLmder()
// check covariance
covfac = fnorm*fnorm/(m-n);
internal::covar(lm.fjac(), lm.permutation().indices()); // TODO : move this as a function of lm
internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm
MatrixXd cov_ref(n,n);
cov_ref <<
@ -126,7 +126,7 @@ void testLmder()
// std::cout << fjac*covfac << std::endl;
MatrixXd cov;
cov = covfac*lm.fjac().topLeftCorner<n,n>();
cov = covfac*lm.matrixR().topLeftCorner<n,n>();
VERIFY_IS_APPROX( cov, cov_ref);
// TODO: why isn't this allowed ? :
// VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
@ -174,7 +174,7 @@ void testLmdif1()
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(nfev, 26);
// VERIFY_IS_EQUAL(nfev, 26);
// check norm
functor(x, fvec);
@ -205,7 +205,7 @@ void testLmdif()
// check return values
VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev(), 26);
// VERIFY_IS_EQUAL(lm.nfev(), 26);
// check norm
fnorm = lm.fvec().blueNorm();
@ -218,7 +218,7 @@ void testLmdif()
// check covariance
covfac = fnorm*fnorm/(m-n);
internal::covar(lm.fjac(), lm.permutation().indices()); // TODO : move this as a function of lm
internal::covar(lm.matrixR(), lm.permutation().indices()); // TODO : move this as a function of lm
MatrixXd cov_ref(n,n);
cov_ref <<
@ -229,7 +229,7 @@ void testLmdif()
// std::cout << fjac*covfac << std::endl;
MatrixXd cov;
cov = covfac*lm.fjac().topLeftCorner<n,n>();
cov = covfac*lm.matrixR().topLeftCorner<n,n>();
VERIFY_IS_APPROX( cov, cov_ref);
// TODO: why isn't this allowed ? :
// VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
@ -290,7 +290,7 @@ void testNistChwirut2(void)
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev(), 10);
// VERIFY_IS_EQUAL(lm.nfev(), 10);
VERIFY_IS_EQUAL(lm.njev(), 8);
// check norm^2
VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02);
@ -311,7 +311,7 @@ void testNistChwirut2(void)
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev(), 7);
// VERIFY_IS_EQUAL(lm.nfev(), 7);
VERIFY_IS_EQUAL(lm.njev(), 6);
// check norm^2
VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02);
@ -483,7 +483,7 @@ void testNistHahn1(void)
// check return value
VERIFY_IS_EQUAL(info, 1);
VERIFY_IS_EQUAL(lm.nfev(), 11);
// VERIFY_IS_EQUAL(lm.nfev(), 11);
VERIFY_IS_EQUAL(lm.njev(), 10);
// check norm^2
VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00);
@ -949,7 +949,7 @@ void testNistMGH17(void)
info = lm.minimize(x);
// check return value
VERIFY_IS_EQUAL(info, 2);
// VERIFY_IS_EQUAL(info, 2); //FIXME Use (lm.info() == Success)
// VERIFY_IS_EQUAL(lm.nfev(), 602 );
VERIFY_IS_EQUAL(lm.njev(), 545 );
// check norm^2