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

[Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>]

Hi,

On Fri, Jun 25, 2010 at 1:50 AM, Ben Goodrich <bgokgm@xxxxxxxxxxxxxx> wrote:
>>> 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.
>

This patch does so. It seems to work okay with both the Pivoting and
NoPivoting options from the previous patch.

>> 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?

Gael?

Thanks,
Ben

# HG changeset patch
# User Ben Goodrich <bgokGM@xxxxxxxxx>
# Date 1277574978 14400
# Node ID 0a9da1ff35cca179862ba9ec8d353cde4fab7f81
# Parent  f02066dadf9a6fbecf13f5d0874600e37130a1a1
LDLT: unit test decompositions of some singular matrices

diff -r f02066dadf9a -r 0a9da1ff35cc test/cholesky.cpp
--- a/test/cholesky.cpp	Fri Jun 25 01:38:49 2010 -0400
+++ b/test/cholesky.cpp	Sat Jun 26 13:56:18 2010 -0400
@@ -64,7 +64,7 @@
   MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols);
   SquareMatrixType symm =  a0 * a0.adjoint();
   // let's make sure the matrix is not singular or near singular
-  for (int k=0; k<3; ++k)
+  if(rows==cols) for (int k=0; k<3; ++k)
   {
     MatrixType a1 = MatrixType::Random(rows,cols);
     symm += a1 * a1.adjoint();
@@ -106,37 +106,48 @@
   #endif
 
   {
-    LLT<SquareMatrixType,Lower> chollo(symmLo);
-    VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
-    vecX = chollo.solve(vecB);
-    VERIFY_IS_APPROX(symm * vecX, vecB);
-    matX = chollo.solve(matB);
-    VERIFY_IS_APPROX(symm * matX, matB);
-
+    if(rows==cols)
+    {
+      LLT<SquareMatrixType,Lower> chollo(symmLo);
+      VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
+      vecX = chollo.solve(vecB);
+      VERIFY_IS_APPROX(symm * vecX, vecB);
+      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);
+    if(rows==cols)
+    {
+      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());
-    vecX = cholup.solve(vecB);
-    VERIFY_IS_APPROX(symm * vecX, vecB);
-    matX = cholup.solve(matB);
-    VERIFY_IS_APPROX(symm * matX, matB);
+    if(rows==cols)
+    {
+      LLT<SquareMatrixType,Upper> cholup(symmUp);
+      VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix());  
+      vecX = cholup.solve(vecB);
+      VERIFY_IS_APPROX(symm * vecX, vecB);
+      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);
+    if(rows==cols)
+    {
+      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);
+    LLT<SquareMatrixType,Lower> chollo(neg);
     VERIFY(chollo.info()==NumericalIssue);
   }
 
@@ -154,19 +165,25 @@
 
     LDLT<SquareMatrixType,Lower> ldltlo(symmLo);
     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);
+    if(rows==cols)
+    {
+      vecX = ldltlo.solve(vecB);
+      VERIFY_IS_APPROX(symm * vecX, vecB);
+      matX = ldltlo.solve(matB);
+      VERIFY_IS_APPROX(symm * matX, matB);
+    }
 
     LDLT<SquareMatrixType,Upper> ldltup(symmUp);
     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(rows==cols)
+    {
+      vecX = ldltup.solve(vecB);
+      VERIFY_IS_APPROX(symm * vecX, vecB);
+      matX = ldltup.solve(matB);
+      VERIFY_IS_APPROX(symm * matX, matB);
+    }
 
-    if(MatrixType::RowsAtCompileTime==Dynamic)
+    if(MatrixType::RowsAtCompileTime==Dynamic && rows==cols)
     {
       // note : each inplace permutation requires a small temporary vector (mask)
 
@@ -214,6 +231,9 @@
     CALL_SUBTEST_5( cholesky(Matrix4d()) );
     CALL_SUBTEST_2( cholesky(MatrixXd(200,200)) );
     CALL_SUBTEST_6( cholesky(MatrixXcd(100,100)) );
+    // test singular matrices
+     CALL_SUBTEST_10( cholesky(MatrixXd(200,200-i-1)) );
+     CALL_SUBTEST_11( cholesky(MatrixXcd(100,100-i-1)) );
   }
 
   CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() );