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

New feature: add rank one update in Cholesky decomposition

Browse Source
This commit is contained in:
Gael Guennebaud 2011-06-20 15:05:50 +02:00
parent a55c27a15f
commit 2f32e48517
2 changed files with 74 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

52
Eigen/src/Cholesky/LLT.h
View File

@ -178,6 +178,9 @@ template<typename _MatrixType, int _UpLo> class LLT
inline Index rows() const { return m_matrix.rows(); } inline Index rows() const { return m_matrix.rows(); }
inline Index cols() const { return m_matrix.cols(); } inline Index cols() const { return m_matrix.cols(); }
template<typename VectorType>
void rankUpdate(const VectorType& vec);
protected: protected:
/** \internal /** \internal
* Used to compute and store L * Used to compute and store L
@ -254,6 +257,34 @@ template<> struct llt_inplace<Lower>
} }
return -1; return -1;
} }
template<typename MatrixType, typename VectorType>
static void rankUpdate(MatrixType& mat, const VectorType& vec)
{
typedef typename MatrixType::ColXpr ColXpr;
typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
typedef typename VectorType::SegmentReturnType VecSegment;
int n = mat.cols();
eigen_assert(mat.rows()==n && vec.size()==n);
typedef typename MatrixType::Scalar Scalar;
Matrix<Scalar,Dynamic,1> temp(vec);
for(int i=0; i<n; ++i)
{
JacobiRotation<Scalar> g;
g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
int rs = n-i-1;
if(rs>0)
{
ColXprSegment x(mat.col(i).tail(rs));
VecSegment y(temp.tail(rs));
apply_rotation_in_the_plane(x, y, g);
}
}
}
}; };
template<> struct llt_inplace<Upper> template<> struct llt_inplace<Upper>
@ -270,6 +301,12 @@ template<> struct llt_inplace<Upper>
Transpose<MatrixType> matt(mat); Transpose<MatrixType> matt(mat);
return llt_inplace<Lower>::blocked(matt); return llt_inplace<Lower>::blocked(matt);
} }
template<typename MatrixType, typename VectorType>
static void rankUpdate(MatrixType& mat, const VectorType& vec)
{
Transpose<MatrixType> matt(mat);
return llt_inplace<Lower>::rankUpdate(matt, vec);
}
}; };
template<typename MatrixType> struct LLT_Traits<MatrixType,Lower> template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
@ -314,6 +351,20 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
return *this; return *this;
} }
/** Performs a rank one update of the current decomposition.
* If A = LL^* before the rank one update,
* then after it we have LL^* = A + vv^* where \a v must be a vector
* of same dimension.
*
*/
template<typename MatrixType, int _UpLo>
template<typename VectorType>
void LLT<MatrixType,_UpLo>::rankUpdate(const VectorType& v)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
internal::llt_inplace<UpLo>::rankUpdate(m_matrix,v);
}
namespace internal { namespace internal {
template<typename _MatrixType, int UpLo, typename Rhs> template<typename _MatrixType, int UpLo, typename Rhs>
struct solve_retval<LLT<_MatrixType, UpLo>, Rhs> struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
@ -384,3 +435,4 @@ SelfAdjointView<MatrixType, UpLo>::llt() const
} }
#endif // EIGEN_LLT_H #endif // EIGEN_LLT_H

23
test/cholesky.cpp
View File

@ -166,6 +166,10 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0); VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0);
VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval()); VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval());
} }
// restore
if(sign == -1)
symm = -symm;
} }
// test some special use cases of SelfCwiseBinaryOp: // test some special use cases of SelfCwiseBinaryOp:
@ -182,6 +186,24 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
m2 = m1; m2 = m1;
m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB); m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB)); VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
// Cholesky update/downdate
{
MatrixType symmLo = symm.template triangularView<Lower>();
MatrixType symmUp = symm.template triangularView<Upper>();
VectorType vec = VectorType::Random(rows);
MatrixType symmCpy = symm + vec * vec.adjoint();
LLT<MatrixType,Lower> chollo(symmLo);
chollo.rankUpdate(vec);
VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix());
LLT<MatrixType,Upper> cholup(symmUp);
cholup.rankUpdate(vec);
VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix());
}
} }
@ -242,7 +264,6 @@ template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
// matX = ldltlo.solve(matB); // matX = ldltlo.solve(matB);
// VERIFY_IS_APPROX(symm * matX, matB); // VERIFY_IS_APPROX(symm * matX, matB);
} }
} }
template<typename MatrixType> void cholesky_verify_assert() template<typename MatrixType> void cholesky_verify_assert()