Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices

[Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>]

Hi,

On Sun, Feb 28, 2010 at 1:44 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
> If I understand correctly your other email, you'd be happy with
> instead a non-pivoting LDLt. That sounds more reasonable than
> "unpivoting" a pivoting LDLt.
>
> So: OK to add a non-pivoting option to LDLt.
>
> Benoit

Here is a unified patch to reimplement the non-pivoting option. Sorry
for the more than 3 month delay; had to write words in dissertation
rather than code :( .

Some thoughts:

0) If you changed your mind and don't want this option anymore, that's
fine. I can apply the patch locally. That said, having a pivot option
makes it possible to be compatible with eigen2 and is probably useful
to someone besides just me. And it is easy to do.

1) Gael made some big changes to LDLT a couple of days ago that I
think I understand but might have misunderstood something when I
reimplemented the non-pivoting option on top of them.

2) For example, in the positive semi-definite but (numerically)
rank-deficient case, I had to put exact zeros on the diagonal of D to
get the reconstruction right. This contravenes Benoit's statement
earlier in this thread that Eigen does not "finish" decompositions of
rank-deficient matrices. But I didn't immediately see any other good
way to do it, in light of Gael's recent changes to LDLT.

3) I implemented the option like

if(pivot)
{
.....
}
else
{
.... // mostly copy-and-paste from above
}

which resulted in some code duplication but avoids the need to
repeatedly evaluate if(pivot) inside the loop. I guess that is
slightly faster and clearer, but if you would rather not have the code
duplication, I can do it the other way.

4) I made some changes to the documentation explaining the pivot
option, but just noticed (and haven't really read yet) today's long
documentation thread.

5) I copied-and-pasted a block inside /test/cholesky.cpp and exercised
the pivot=false option. It seems to work when you do ./check.sh
cholesky. I did some other tests locally with singular matrices, but
/test/cholesky.cpp does not seem to have any tests with singular
matrices, so maybe some should be added?

Thanks,
Ben

# HG changeset patch
# User Ben Goodrich <bgokGM@xxxxxxxxx>
# Date 1275793116 14400
# Branch pivot
# Node ID b65ed5328b2d5df3da0c6cd0dcc54dbe30ec5bba
# Parent  343f266cff99038dd710acaf0cacb88595930671
LDLT: Introduce and test pivot option, which defaults to true, but when false decomposes
the matrix without pivoting, similar to eigen 2 behavior.

diff -r 343f266cff99 -r b65ed5328b2d Eigen/src/Cholesky/LDLT.h
--- a/Eigen/src/Cholesky/LDLT.h	Fri Jun 04 23:17:57 2010 +0200
+++ b/Eigen/src/Cholesky/LDLT.h	Sat Jun 05 22:58:36 2010 -0400
@@ -33,7 +33,7 @@
   *
   * \class LDLT
   *
-  * \brief Robust Cholesky decomposition of a matrix with pivoting
+  * \brief Robust Cholesky decomposition of a matrix with (by default) pivoting
   *
   * \param MatrixType the type of the matrix of which to compute the LDL^T Cholesky decomposition
   *
@@ -41,12 +41,14 @@
   * matrix \f$ A \f$ such that \f$ A =  P^TLDL^*P \f$, where P is a permutation matrix, L
   * is lower triangular with a unit diagonal and D is a diagonal matrix.
   *
-  * The decomposition uses pivoting to ensure stability, so that L will have
+  * The decomposition uses pivoting (by default) to ensure stability, so that L will have
   * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root
-  * on D also stabilizes the computation.
+  * on D also stabilizes the computation. If the pivoting option is not utilized, \f$ A \f$
+  * must be positive semi-definite by construction, which is NOT checked. In that positive semi-definite
+  * case, numerical problems may arise if \f$ A \f$ is rank-deficient and pivoting is not utilized.
   *
   * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky
-	* decomposition to determine whether a system of equations has a solution.
+  * decomposition to determine whether a system of equations has a solution.
   *
   * \sa MatrixBase::ldlt(), class LLT
   */
@@ -97,13 +99,13 @@
         m_isInitialized(false)
     {}
 
-    LDLT(const MatrixType& matrix)
+    LDLT(const MatrixType& matrix, const bool pivot=true)
       : m_matrix(matrix.rows(), matrix.cols()),
         m_transpositions(matrix.rows()),
         m_temporary(matrix.rows()),
         m_isInitialized(false)
     {
-      compute(matrix);
+      compute(matrix, pivot);
     }
 
     /** \returns a view of the upper triangular matrix U */
@@ -168,7 +170,7 @@
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
 
-    LDLT& compute(const MatrixType& matrix);
+    LDLT& compute(const MatrixType& matrix, const bool pivot=true);
 
     /** \returns the internal LDLT decomposition matrix
       *
@@ -205,7 +207,8 @@
 template<> struct ei_ldlt_inplace<Lower>
 {
   template<typename MatrixType, typename TranspositionType, typename Workspace>
-  static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
+  static bool unblocked(MatrixType& mat, TranspositionType& transpositions, 
+                        Workspace& temp, int* sign=0, bool pivot=true)
   {
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
@@ -222,11 +225,11 @@
     }
 
     RealScalar cutoff = 0, biggest_in_corner;
+    Index index_of_biggest_in_corner;
 
-    for (Index k = 0; k < size; ++k)
+    if(pivot) for (Index k = 0; k < size; ++k)
     {
       // Find largest diagonal element
-      Index index_of_biggest_in_corner;
       biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
       index_of_biggest_in_corner += k;
 
@@ -270,6 +273,49 @@
       if((rs>0) && (ei_abs(mat.coeffRef(k,k)) > cutoff))
         A21 /= mat.coeffRef(k,k);
     }
+    else for (Index k = 0; k < size; ++k) // no pivoting, which assumes matrix is positive semi-definite
+    {
+      if(k == 0)
+      {
+        biggest_in_corner = mat.diagonal().tail(size).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
+
+        // The biggest overall is the point of reference to which further diagonals
+        // are compared; if any diagonal is negligible compared
+        // to the largest overall, the algorithm bails.
+        cutoff = ei_abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
+
+        if(sign)
+          *sign = 1; // this is by assumption rather than calculated
+      }
+      
+      Index rs = size - k - 1;
+      Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
+      Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
+      Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
+
+      if(k>0)
+      {
+        temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
+        mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
+        if(rs>0)
+          A21.noalias() -= A20 * temp.head(k);
+      }
+      if(ei_abs(mat.coeffRef(k,k)) > cutoff)
+      {
+        if(rs>0)
+          A21 /= mat.coeffRef(k,k);
+      }
+      else // Finish early if the matrix is not full rank.
+      {
+        for(Index i = k; i < size; i++)
+        {
+          transpositions.coeffRef(i) = i;
+	  mat(i,i) = 0;
+	}
+        break;
+      }
+      transpositions.coeffRef(k) = k;
+    }
 
     return true;
   }
@@ -278,10 +324,11 @@
 template<> struct ei_ldlt_inplace<Upper>
 {
   template<typename MatrixType, typename TranspositionType, typename Workspace>
-  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
+  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, 
+                                            Workspace& temp, int* sign=0, bool pivot=true)
   {
     Transpose<MatrixType> matt(mat);
-    return ei_ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
+    return ei_ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign, pivot);
   }
 };
 
@@ -304,7 +351,7 @@
 /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
   */
 template<typename MatrixType, int _UpLo>
-LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
+LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a, const bool pivot)
 {
   ei_assert(a.rows()==a.cols());
   const Index size = a.rows();
@@ -315,7 +362,7 @@
   m_isInitialized = false;
   m_temporary.resize(size);
 
-  ei_ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
+  ei_ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign, pivot);
 
   m_isInitialized = true;
   return *this;
@@ -398,23 +445,23 @@
 }
 
 /** \cholesky_module
-  * \returns the Cholesky decomposition with full pivoting without square root of \c *this
+  * \returns the Cholesky decomposition with (by default) full pivoting without square root of \c *this
   */
 template<typename MatrixType, unsigned int UpLo>
 inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo>
-SelfAdjointView<MatrixType, UpLo>::ldlt() const
+SelfAdjointView<MatrixType, UpLo>::ldlt(const bool pivot) const
 {
-  return LDLT<PlainObject,UpLo>(m_matrix);
+  return LDLT<PlainObject,UpLo>(m_matrix, pivot);
 }
 
 /** \cholesky_module
-  * \returns the Cholesky decomposition with full pivoting without square root of \c *this
+  * \returns the Cholesky decomposition with (by default) full pivoting without square root of \c *this
   */
 template<typename Derived>
 inline const LDLT<typename MatrixBase<Derived>::PlainObject>
-MatrixBase<Derived>::ldlt() const
+MatrixBase<Derived>::ldlt(const bool pivot) const
 {
-  return LDLT<PlainObject>(derived());
+  return LDLT<PlainObject>(derived(), pivot);
 }
 
 #endif // EIGEN_LDLT_H
diff -r 343f266cff99 -r b65ed5328b2d Eigen/src/Core/MatrixBase.h
--- a/Eigen/src/Core/MatrixBase.h	Fri Jun 04 23:17:57 2010 +0200
+++ b/Eigen/src/Core/MatrixBase.h	Sat Jun 05 22:58:36 2010 -0400
@@ -299,7 +299,7 @@
 /////////// Cholesky module ///////////
 
     const LLT<PlainObject>  llt() const;
-    const LDLT<PlainObject> ldlt() const;
+    const LDLT<PlainObject> ldlt(const bool pivot=true) const;
 
 /////////// QR module ///////////
 
diff -r 343f266cff99 -r b65ed5328b2d Eigen/src/Core/SelfAdjointView.h
--- a/Eigen/src/Core/SelfAdjointView.h	Fri Jun 04 23:17:57 2010 +0200
+++ b/Eigen/src/Core/SelfAdjointView.h	Sat Jun 05 22:58:36 2010 -0400
@@ -153,7 +153,7 @@
 /////////// Cholesky module ///////////
 
     const LLT<PlainObject, UpLo> llt() const;
-    const LDLT<PlainObject, UpLo> ldlt() const;
+    const LDLT<PlainObject, UpLo> ldlt(const bool pivot=true) const;
 
 /////////// Eigenvalue module ///////////
 
diff -r 343f266cff99 -r b65ed5328b2d test/cholesky.cpp
--- a/test/cholesky.cpp	Fri Jun 04 23:17:57 2010 +0200
+++ b/test/cholesky.cpp	Sat Jun 05 22:58:36 2010 -0400
@@ -120,6 +120,38 @@
     matX = cholup.solve(matB);
     VERIFY_IS_APPROX(symm * matX, matB);
   }
+  
+  { // now check LDLT *without* pivoting
+    LDLT<SquareMatrixType,Lower> ldltlo(symm, false);
+    VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix());
+    vecX = ldltlo.solve(vecB);
+    VERIFY_IS_APPROX(symm * vecX, vecB);
+    matX = ldltlo.solve(matB);
+    VERIFY_IS_APPROX(symm * matX, matB);
+
+    LDLT<SquareMatrixType,Upper> ldltup(symm, false);
+    VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix());
+    vecX = ldltup.solve(vecB);
+    VERIFY_IS_APPROX(symm * vecX, vecB);
+    matX = ldltup.solve(matB);
+    VERIFY_IS_APPROX(symm * matX, matB);
+
+    if(MatrixType::RowsAtCompileTime==Dynamic)
+    {
+      // note : each inplace permutation requires a small temporary vector (mask)
+
+      // check inplace solve
+      matX = matB;
+      VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0);
+      VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval());
+
+
+      matX = matB;
+      VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0);
+      VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval());
+    }
+  }
+
 
   int sign = ei_random<int>()%2 ? 1 : -1;