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

[Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>]

Hi,

On Sat, Jun 5, 2010 at 11:34 PM, 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 :( .

Has anyone tried my patch or has comments on the things below? I
rebased the patch in the attachment to make sure it would still apply
and to improve the documentation a little more. -- Ben

> 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 1276629596 14400
# Branch pivot
# Node ID c6bf0801d6e0f89fff5d7168ce97ce00e70bf58c
# Parent  ee38d0fe25eca4b019f2544ab05c1aa09c213d67
LDLT: Introduce and test pivot option, which defaults to true, but when false decomposes
the matrix without pivoting, similar to eigen2 behavior.

diff -r ee38d0fe25ec -r c6bf0801d6e0 Eigen/src/Cholesky/LDLT.h
--- a/Eigen/src/Cholesky/LDLT.h	Sat Jun 05 22:58:36 2010 -0400
+++ b/Eigen/src/Cholesky/LDLT.h	Tue Jun 15 15:19:56 2010 -0400
@@ -44,7 +44,7 @@
   * 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. 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
+  * must be positive semi-definite by construction, which is NOT checked. In the 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
@@ -98,6 +98,11 @@
         m_isInitialized(false)
     {}
 
+    /** \brief Default Constructor with decomposition
+      *
+      * Performs the decomposition of \a matrix with pivoting by default, which
+      * can be changed by setting \a pivot to false.
+      */
     LDLT(const MatrixType& matrix, const bool pivot=true)
       : m_matrix(matrix.rows(), matrix.cols()),
         m_transpositions(matrix.rows()),
@@ -302,7 +307,11 @@
         if(sign)
           *sign = 1; // this is by assumption rather than calculated
       }
-      
+
+      // partition the matrix:
+      //       A00 |  -  |  -
+      // lu  = A10 | A11 |  -
+      //       A20 | A21 | A22
       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);
@@ -364,6 +373,7 @@
 };
 
 /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
+  * with pivoting by default.
   */
 template<typename MatrixType, int _UpLo>
 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a, const bool pivot)