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

[Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>]

Hi,

Thank you for creating eigen. It is great.

I mentioned a bug about ldlt() with rank-deficient matrices and gave a
test case here:

http://forum.kde.org/viewtopic.php?f=74&t=85449&sid=43724496b3e1705a95b41fa17c0f95eb&start=10

I believe the problem is that when ldlt() decides to abort due to a
diagonal entry of D being too close to zero, it does not change the
remaining diagonal entries to be 0.0 . As a result, when the original
matrix has a nullity of 2 or more, multiplying the factors together
does not result in the original matrix.

The attached patch fixes this issue for me in the test case by zeroing
out the relevant diagonal entries. It applies to the development
branch. I will reply to this message with a patch and a question about
the stable branch.

Thank you,
Ben

# HG changeset patch
# User Ben Goodrich <bgokGM@xxxxxxxxx>
# Date 1266888211 18000
# Node ID 1c47a9153fef51ef667f90910de4cb86474d8fdc
# Parent  a4b00df290033caece6c90dfdd9c214d5bc20e51
set some diagonal elements of D as 0.0 if matrix is rank-deficient

diff -r a4b00df29003 -r 1c47a9153fef Eigen/src/Cholesky/LDLT.h
--- a/Eigen/src/Cholesky/LDLT.h	Mon Feb 22 15:39:17 2010 +0100
+++ b/Eigen/src/Cholesky/LDLT.h	Mon Feb 22 20:23:31 2010 -0500
@@ -214,7 +214,10 @@
     // Finish early if the matrix is not full rank.
     if(biggest_in_corner < cutoff)
     {
-      for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i;
+      for(int i = j; i < size; i++) {
+        m_transpositions.coeffRef(i) = i;
+        m_matrix.coeffRef(i,i) = 0.0;
+      }
       break;
     }
 
@@ -238,7 +241,10 @@
     // Finish early if the matrix is not full rank.
     if(ei_abs(Djj) < cutoff)
     {
-      for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i;
+      for(int i = j; i < size; i++) {
+        m_transpositions.coeffRef(i) = i;
+        m_matrix.coeffRef(i,i) = 0.0;
+      }
       break;
     }