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

[Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>]

Hi Gael and everyone else,

On Wed, Jun 16, 2010 at 1:05 PM, Gael Guennebaud
<gael.guennebaud@xxxxxxxxx> wrote:
> On Sun, Jun 6, 2010 at 5:34 AM, Ben Goodrich <bgokgm@xxxxxxxxxxxxxx> wrote:
>> 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.
>
> When I redesigned LDLT to support blocking and to work fully inplace,
> I also thought it might be cool to have such a pivoting option,
> however I was unsure whether this should be a runtime and or compile
> time option. And if it is a runtime option, shall it be set once for
> all by the ctor, or modifiable when calling compute, or via a specific
> function... Finally, we can consider this option to be in pair with
> the compute* options for the other decompositions, so let's do the
> same here, i.e.:
>
> add a "int option" parameter to all LDLT ctor and to the compute
> method defaulting to "Pivoting". To disable pivoting will use the
> NoPivoting constant. That's better than meaningless booleans.

Here is a patch that uses the options=Pivoting or options=NoPivoting
syntax instead of pivoting=true or pivoting=false and avoids code
duplication.

>> 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.
>
> maybe this is fixed now ?

I am not sure I understand you. As far as I know, it was not broken in
the Pivoting case, but the decomposition is done in-place now. So, if
we finish early in the NoPivoting case, then we have to put zeros on
the diagonal to get vectorD() to extract the correct thing, right?

>> 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?
>
> why not.

I have not added the singular tests yet, but I can do that soon.

> You should also make the solve function skips the transpositions when
> no pivoting has been computed.

I have not done this yet either. What did you decide about making a
new class versus putting a flag in the class definition?

Ben

# HG changeset patch
# User Ben Goodrich <bgokGM@xxxxxxxxx>
# Date 1277444329 14400
# Node ID f02066dadf9a6fbecf13f5d0874600e37130a1a1
# Parent  e2c98a46ac5c2ec769a33efee9480b02d3f46af5
LDLT: Introduce and test the option to pivot

diff -r e2c98a46ac5c -r f02066dadf9a Eigen/src/Cholesky/LDLT.h
--- a/Eigen/src/Cholesky/LDLT.h	Thu Jun 24 23:21:58 2010 +0200
+++ b/Eigen/src/Cholesky/LDLT.h	Fri Jun 25 01:38:49 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 pivoting (by default)
   *
   * \param MatrixType the type of the matrix of which to compute the LDL^T Cholesky decomposition
   *
@@ -41,12 +41,15 @@
   * 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.
   *
+  * The option to avoid pivoting should only be used if \f$ A \f$ is positive semi-definite and
+  * may result in numerical inaccuracy when \f$ A \f$ is rank-deficient.
+  *
   * 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
   */
@@ -96,13 +99,20 @@
         m_isInitialized(false)
     {}
 
-    LDLT(const MatrixType& matrix)
+    /** \brief Constructor with computation
+      *
+      * Performs the decomposition of the \a matrix, with the \param options
+      * of Pivoting (default) or NoPivoting.
+      * \sa LDLT()
+      */
+
+    LDLT(const MatrixType& matrix, int options=Pivoting)
       : m_matrix(matrix.rows(), matrix.cols()),
         m_transpositions(matrix.rows()),
         m_temporary(matrix.rows()),
         m_isInitialized(false)
     {
-      compute(matrix);
+      compute(matrix, options);
     }
 
     /** \returns a view of the upper triangular matrix U */
@@ -167,7 +177,7 @@
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
 
-    LDLT& compute(const MatrixType& matrix);
+    LDLT& compute(const MatrixType& matrix, int options=Pivoting);
 
     /** \returns the internal LDLT decomposition matrix
       *
@@ -204,7 +214,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, int options=Pivoting)
   {
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
@@ -220,51 +231,52 @@
       return true;
     }
 
-    RealScalar cutoff = 0, 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.
+    Index index_of_biggest_in_corner;
+    RealScalar biggest_in_corner = mat.diagonal().tail(size).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
+    RealScalar cutoff = ei_abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
+
+    if(sign)
+      *sign = ei_real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
 
     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;
+      
+      if(options==Pivoting)
+      {
+        // Find largest diagonal element
+        biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner);
+        index_of_biggest_in_corner += k;
 
-      if(k == 0)
-      {
-        // 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);
+        // Finish early if the matrix is not full rank.
+        if(biggest_in_corner < cutoff)
+        {
+          for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
+          break;
+        }
+        else transpositions.coeffRef(k) = index_of_biggest_in_corner;
 
-        if(sign)
-          *sign = ei_real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
+      
+        if(k != index_of_biggest_in_corner)
+        {
+          // apply the transposition while taking care to consider only
+          // the lower triangular part
+          Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
+          mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
+          mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
+          std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
+          for(int i=k+1;i<index_of_biggest_in_corner;++i)
+          {
+            Scalar tmp = mat.coeffRef(i,k);
+            mat.coeffRef(i,k) = ei_conj(mat.coeffRef(index_of_biggest_in_corner,i));
+            mat.coeffRef(index_of_biggest_in_corner,i) = ei_conj(tmp);
+          }
+          if(NumTraits<Scalar>::IsComplex)
+            mat.coeffRef(index_of_biggest_in_corner,k) = ei_conj(mat.coeff(index_of_biggest_in_corner,k));
+        }
       }
-
-      // Finish early if the matrix is not full rank.
-      if(biggest_in_corner < cutoff)
-      {
-        for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
-        break;
-      }
-
-      transpositions.coeffRef(k) = index_of_biggest_in_corner;
-      if(k != index_of_biggest_in_corner)
-      {
-        // apply the transposition while taking care to consider only
-        // the lower triangular part
-        Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
-        mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
-        mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
-        std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
-        for(int i=k+1;i<index_of_biggest_in_corner;++i)
-        {
-          Scalar tmp = mat.coeffRef(i,k);
-          mat.coeffRef(i,k) = ei_conj(mat.coeffRef(index_of_biggest_in_corner,i));
-          mat.coeffRef(index_of_biggest_in_corner,i) = ei_conj(tmp);
-        }
-        if(NumTraits<Scalar>::IsComplex)
-          mat.coeffRef(index_of_biggest_in_corner,k) = ei_conj(mat.coeff(index_of_biggest_in_corner,k));
-      }
+      else transpositions.coeffRef(k) = k;
 
       // partition the matrix:
       //       A00 |  -  |  -
@@ -282,10 +294,21 @@
         if(rs>0)
           A21.noalias() -= A20 * temp.head(k);
       }
-      if((rs>0) && (ei_abs(mat.coeffRef(k,k)) > cutoff))
-        A21 /= mat.coeffRef(k,k);
+      if(rs>0)
+      {
+        if(ei_abs(mat.coeffRef(k,k)) > cutoff)
+          A21 /= mat.coeffRef(k,k);
+        else if(options==NoPivoting) // finish early
+        {
+          for(Index i = k + 1; i < size; i++)
+          {
+            transpositions.coeffRef(i) = i;
+            mat.coeffRef(i,i) = 0; // NoPivoting assumes matrix is PSD
+          }
+          break;
+        }
+      }
     }
-
     return true;
   }
 };
@@ -293,10 +316,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, int options=Pivoting)
   {
     Transpose<MatrixType> matt(mat);
-    return ei_ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
+    return ei_ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign, options);
   }
 };
 
@@ -317,9 +341,10 @@
 };
 
 /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
+  * with the \param options of Pivoting (default) or NoPivoting.
   */
 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, int options)
 {
   ei_assert(a.rows()==a.cols());
   const Index size = a.rows();
@@ -330,7 +355,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, options);
 
   m_isInitialized = true;
   return *this;
@@ -413,23 +438,23 @@
 }
 
 /** \cholesky_module
-  * \returns the Cholesky decomposition with full pivoting without square root of \c *this
+  * \returns the Cholesky decomposition with full pivoting (by default) 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(int options) const
 {
-  return LDLT<PlainObject,UpLo>(m_matrix);
+  return LDLT<PlainObject,UpLo>(m_matrix, options);
 }
 
 /** \cholesky_module
-  * \returns the Cholesky decomposition with full pivoting without square root of \c *this
+  * \returns the Cholesky decomposition with full pivoting (by default) without square root of \c *this
   */
 template<typename Derived>
 inline const LDLT<typename MatrixBase<Derived>::PlainObject>
-MatrixBase<Derived>::ldlt() const
+MatrixBase<Derived>::ldlt(int options) const
 {
-  return LDLT<PlainObject>(derived());
+  return LDLT<PlainObject>(derived(), options);
 }
 
 #endif // EIGEN_LDLT_H
diff -r e2c98a46ac5c -r f02066dadf9a Eigen/src/Core/MatrixBase.h
--- a/Eigen/src/Core/MatrixBase.h	Thu Jun 24 23:21:58 2010 +0200
+++ b/Eigen/src/Core/MatrixBase.h	Fri Jun 25 01:38:49 2010 -0400
@@ -313,7 +313,7 @@
 /////////// Cholesky module ///////////
 
     const LLT<PlainObject>  llt() const;
-    const LDLT<PlainObject> ldlt() const;
+    const LDLT<PlainObject> ldlt(int options=Pivoting) const;
 
 /////////// QR module ///////////
 
diff -r e2c98a46ac5c -r f02066dadf9a Eigen/src/Core/SelfAdjointView.h
--- a/Eigen/src/Core/SelfAdjointView.h	Thu Jun 24 23:21:58 2010 +0200
+++ b/Eigen/src/Core/SelfAdjointView.h	Fri Jun 25 01:38:49 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(int options=Pivoting) const;
 
 /////////// Eigenvalue module ///////////
 
diff -r e2c98a46ac5c -r f02066dadf9a test/cholesky.cpp
--- a/test/cholesky.cpp	Thu Jun 24 23:21:58 2010 +0200
+++ b/test/cholesky.cpp	Fri Jun 25 01:38:49 2010 -0400
@@ -113,6 +113,13 @@
     matX = chollo.solve(matB);
     VERIFY_IS_APPROX(symm * matX, matB);
 
+    LDLT<SquareMatrixType,Lower> ldltlo(symmLo, NoPivoting);
+    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);
+    
     // test the upper mode
     LLT<SquareMatrixType,Upper> cholup(symmUp);
     VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix());
@@ -121,12 +128,19 @@
     matX = cholup.solve(matB);
     VERIFY_IS_APPROX(symm * matX, matB);
 
+    LDLT<SquareMatrixType,Upper> ldltup(symmUp, NoPivoting);
+    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);
+    
     MatrixType neg = -symmLo;
     chollo.compute(neg);
     VERIFY(chollo.info()==NumericalIssue);
   }
 
-  // LDLT
+  // more LDLT
   {
     int sign = ei_random<int>()%2 ? 1 : -1;