Re: [eigen] Instability in LLT and LDLT methods.

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

OK so now I have a patch (attached) making LU use only partial pivoting.

Results:

* for an invertible matrix, partial pivoting seems to be OK. So
partial pivoting can be considered for implementing inverse(), where
we assume the matrix to be invertible.

* for a non-invertible matrix, when one wants to compute the rank or
things that depend on it, such as a basis of the kernel, partial
pivoting is not enough.

E.g. apply my patch, and try the test_lu unit test with 100
repetitions, it fails.

So yes it would be potentially  interesting to add back the template
parameter "int Pivoting" to class LU and let inverse() use a partial
pivoting LU. All rank-based methods (such as kernel() and image())
should only be available in the full pivoting version.

Cheers,
Benoit

Index: Eigen/src/LU/LU.h
===================================================================
--- Eigen/src/LU/LU.h	(révision 917831)
+++ Eigen/src/LU/LU.h	(copie de travail)
@@ -340,34 +340,22 @@
   IntRowVectorType cols_transpositions(matrix.cols());
   int number_of_transpositions = 0;
 
-  RealScalar biggest = RealScalar(0);
+  RealScalar biggest = m_lu.cwise().abs().maxCoeff();
   m_rank = size;
   for(int k = 0; k < size; ++k)
   {
     int row_of_biggest_in_corner, col_of_biggest_in_corner;
     RealScalar biggest_in_corner;
 
-    biggest_in_corner = m_lu.corner(Eigen::BottomRight, rows-k, cols-k)
+    biggest_in_corner = m_lu.block(k,k, rows-k, 1)
                         .cwise().abs()
-                        .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
+                        .maxCoeff(&row_of_biggest_in_corner);
     row_of_biggest_in_corner += k;
-    col_of_biggest_in_corner += k;
-    if(k==0) biggest = biggest_in_corner;
+    col_of_biggest_in_corner = k;
 
-    // if the corner is negligible, then we have less than full rank, and we can finish early
-    if(ei_isMuchSmallerThan(biggest_in_corner, biggest))
-    {
-      m_rank = k;
-      for(int i = k; i < size; i++)
-      {
-        rows_transpositions.coeffRef(i) = i;
-        cols_transpositions.coeffRef(i) = i;
-      }
-      break;
-    }
-
     rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
     cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
+
     if(k != row_of_biggest_in_corner) {
       m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
       ++number_of_transpositions;
@@ -376,6 +364,13 @@
       m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
       ++number_of_transpositions;
     }
+
+    // if the corner is negligible, then we have less than full rank, and we can finish early
+    if(ei_isMuchSmallerThan(biggest_in_corner, biggest))
+    {
+      --m_rank; continue;
+    }
+
     if(k<rows-1)
       m_lu.col(k).end(rows-k-1) /= m_lu.coeff(k,k);
     if(k<size-1)