aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
-rw-r--r--Eigen/src/Cholesky/LDLT.h17
-rw-r--r--Eigen/src/Cholesky/LLT.h90
-rw-r--r--test/cholesky.cpp70
3 files changed, 118 insertions, 59 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 3f0e63b68..6a2d2994b 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -39,6 +39,8 @@ template<typename MatrixType, int UpLo> struct LDLT_Traits;
* \brief Robust Cholesky decomposition of a matrix with pivoting
*
* \param MatrixType the type of the matrix of which to compute the LDL^T Cholesky decomposition
+ * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
+ * The other triangular part won't be read.
*
* Perform a robust Cholesky decomposition of a positive semidefinite or negative semidefinite
* matrix \f$ A \f$ such that \f$ A = P^TLDL^*P \f$, where P is a permutation matrix, L
@@ -53,10 +55,6 @@ template<typename MatrixType, int UpLo> struct LDLT_Traits;
*
* \sa MatrixBase::ldlt(), class LLT
*/
- /* THIS PART OF THE DOX IS CURRENTLY DISABLED BECAUSE INACCURATE BECAUSE OF BUG IN THE DECOMPOSITION CODE
- * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
- * the strict lower part does not have to store correct values.
- */
template<typename _MatrixType, int _UpLo> class LDLT
{
public:
@@ -228,6 +226,17 @@ template<typename _MatrixType, int _UpLo> class LDLT
inline Index rows() const { return m_matrix.rows(); }
inline Index cols() const { return m_matrix.cols(); }
+ /** \brief Reports whether previous computation was successful.
+ *
+ * \returns \c Success if computation was succesful,
+ * \c NumericalIssue if the matrix.appears to be negative.
+ */
+ ComputationInfo info() const
+ {
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
+ return Success;
+ }
+
protected:
/** \internal
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 0c7093375..2140f3d5c 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -36,6 +36,8 @@ template<typename MatrixType, int UpLo> struct LLT_Traits;
* \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
*
* \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
+ * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
+ * The other triangular part won't be read.
*
* This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
* matrix A such that A = LL^* = U^*U, where L is lower triangular.
@@ -182,7 +184,7 @@ template<typename _MatrixType, int _UpLo> class LLT
inline Index cols() const { return m_matrix.cols(); }
template<typename VectorType>
- void rankUpdate(const VectorType& vec);
+ LLT& rankUpdate(const VectorType& vec, const RealScalar& sigma = 1);
protected:
/** \internal
@@ -200,11 +202,11 @@ template<typename Scalar, int UpLo> struct llt_inplace;
template<typename Scalar> struct llt_inplace<Scalar, Lower>
{
+ typedef typename NumTraits<Scalar>::Real RealScalar;
template<typename MatrixType>
static typename MatrixType::Index unblocked(MatrixType& mat)
{
typedef typename MatrixType::Index Index;
- typedef typename MatrixType::RealScalar RealScalar;
eigen_assert(mat.rows()==mat.cols());
const Index size = mat.rows();
@@ -261,8 +263,9 @@ template<typename Scalar> struct llt_inplace<Scalar, Lower>
}
template<typename MatrixType, typename VectorType>
- static void rankUpdate(MatrixType& mat, const VectorType& vec)
+ static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
{
+ typedef typename MatrixType::Index Index;
typedef typename MatrixType::ColXpr ColXpr;
typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
@@ -271,26 +274,67 @@ template<typename Scalar> struct llt_inplace<Scalar, Lower>
int n = mat.cols();
eigen_assert(mat.rows()==n && vec.size()==n);
- TempVectorType temp(vec);
- for(int i=0; i<n; ++i)
+ TempVectorType temp;
+
+ if(sigma>0)
{
- JacobiRotation<Scalar> g;
- g.makeGivens(mat(i,i), -temp(i), &mat(i,i));
+ // This version is based on Givens rotations.
+ // It is faster than the other one below, but only works for updates,
+ // i.e., for sigma > 0
+ temp = sqrt(sigma) * vec;
- int rs = n-i-1;
- if(rs>0)
+ 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));
+ TempVecSegment y(temp.tail(rs));
+ apply_rotation_in_the_plane(x, y, g);
+ }
+ }
+ }
+ else
+ {
+ temp = vec;
+ RealScalar beta = 1;
+ for(int j=0; j<n; ++j)
{
- ColXprSegment x(mat.col(i).tail(rs));
- TempVecSegment y(temp.tail(rs));
- apply_rotation_in_the_plane(x, y, g);
+ RealScalar Ljj = real(mat.coeff(j,j));
+ RealScalar dj = abs2(Ljj);
+ Scalar wj = temp.coeff(j);
+ RealScalar swj2 = sigma*abs2(wj);
+ RealScalar gamma = dj*beta + swj2;
+
+ RealScalar x = dj + swj2/beta;
+ if (x<=RealScalar(0))
+ return j;
+ RealScalar nLjj = sqrt(x);
+ mat.coeffRef(j,j) = nLjj;
+ beta += swj2/dj;
+
+ // Update the terms of L
+ Index rs = n-j-1;
+ if(rs)
+ {
+ temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
+ if(gamma != 0)
+ mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs);
+ }
}
}
+ return -1;
}
};
template<typename Scalar> struct llt_inplace<Scalar, Upper>
{
+ typedef typename NumTraits<Scalar>::Real RealScalar;
+
template<typename MatrixType>
static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
{
@@ -304,10 +348,10 @@ template<typename Scalar> struct llt_inplace<Scalar, Upper>
return llt_inplace<Scalar, Lower>::blocked(matt);
}
template<typename MatrixType, typename VectorType>
- static void rankUpdate(MatrixType& mat, const VectorType& vec)
+ static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma)
{
Transpose<MatrixType> matt(mat);
- return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate());
+ return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma);
}
};
@@ -343,7 +387,7 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
template<typename MatrixType, int _UpLo>
LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
{
- assert(a.rows()==a.cols());
+ eigen_assert(a.rows()==a.cols());
const Index size = a.rows();
m_matrix.resize(size, size);
m_matrix = a;
@@ -355,18 +399,24 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
return *this;
}
-/** Performs a rank one update of the current decomposition.
+/** Performs a rank one update (or dowdate) 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
+ * then after it we have LL^* = A + sigma * v v^* 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)
+LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)
{
EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType);
- internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v);
+ eigen_assert(v.size()==m_matrix.cols());
+ eigen_assert(m_isInitialized);
+ if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0)
+ m_info = NumericalIssue;
+ else
+ m_info = Success;
+
+ return *this;
}
namespace internal {
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index d9806e5c3..1a1b2eeb5 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -41,6 +41,38 @@ static int nb_temporaries;
VERIFY( (#XPR) && nb_temporaries==N ); \
}
+template<typename MatrixType,template <typename,int> class CholType> void test_chol_update(const MatrixType& symm)
+{
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
+
+ MatrixType symmLo = symm.template triangularView<Lower>();
+ MatrixType symmUp = symm.template triangularView<Upper>();
+ MatrixType symmCpy = symm;
+
+ CholType<MatrixType,Lower> chollo(symmLo);
+ CholType<MatrixType,Upper> cholup(symmUp);
+
+ for (int k=0; k<10; ++k)
+ {
+ VectorType vec = VectorType::Random(symm.rows());
+ RealScalar sigma = internal::random<RealScalar>();
+ symmCpy += sigma * vec * vec.adjoint();
+
+ // we are doing some downdates, so it might be the case that the matrix is not SPD anymore
+ CholType<MatrixType,Lower> chol(symmCpy);
+ if(chol.info()!=Success)
+ break;
+
+ chollo.rankUpdate(vec, sigma);
+ VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix());
+
+ cholup.rankUpdate(vec, sigma);
+ VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix());
+ }
+}
+
template<typename MatrixType> void cholesky(const MatrixType& m)
{
typedef typename MatrixType::Index Index;
@@ -155,41 +187,9 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
- // LLT 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());
- }
-
- // LDLT 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();
-
- LDLT<MatrixType,Lower> chollo(symmLo);
- chollo.rankUpdate(vec);
- VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix());
-
- LDLT<MatrixType,Upper> cholup(symmUp);
- cholup.rankUpdate(vec);
- VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix());
- }
+ // update/downdate
+ CALL_SUBTEST(( test_chol_update<SquareMatrixType,LLT>(symm) ));
+ CALL_SUBTEST(( test_chol_update<SquareMatrixType,LDLT>(symm) ));
}
template<typename MatrixType> void cholesky_cplx(const MatrixType& m)