Re: [eigen] [RFC] Full pivoting LDLT

[Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>]

On Tue, 3 Feb 2009, Keir Mierle wrote:

> http://codereview.appspot.com/14042/diff/1/4
>
> There are a couple of issues.
>
> 1. Pivoting is not optional. Any objections?
> 2. If the matrix is only semidefinite, then matrixL() gives wrong results
> (1's in diagonal where they shouldn't be).
>
> Minus objections, I'll commit this.

I don't understand issue 2. I thought the L in the LDL^T decomposition has 
per definition 1s on the diagonal?

It looks like you assume that the matrix is positive semi-definite. 
However, I believe that the LDL^T decomposition is also useful for 
symmetric indefinite matrices. Solving Ax = b with A symmetric is twice as 
fast with LDL^T than with LU (if I remember correctly).

Why do you store P in two places: m_p and m_transpositions? Is using 
m_transpositions faster than m_p?

I was a bit surprised that you use both the triangles in m_matrix during 
the computation, since one half is just a rescaling of the other half. So 
I tried using only the lower triangle, but my timings (using g++ 4.1) are 
inconclusive on whether this is actually faster. I attach the diff for 
your entertainment, but I don't have any good reason for it.

Cheers,
Jitse

Index: Eigen/src/Cholesky/LDLT.h
===================================================================
--- Eigen/src/Cholesky/LDLT.h	(revision 921039)
+++ Eigen/src/Cholesky/LDLT.h	(working copy)
@@ -171,13 +171,14 @@
     }
 
     if (j == 0) {
-      m_matrix.row(0) = m_matrix.row(0).conjugate();
-      m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix.coeff(0,0);
+      m_matrix.col(0).end(size-1) /= m_matrix.coeff(0,0);
       continue;
     }
 
+    _temporary.start(j) = ( m_matrix.row(j).start(j).cwise() 
+			    * m_matrix.diagonal().start(j).transpose() ).conjugate();
     RealScalar Djj = ei_real(m_matrix.coeff(j,j) - (m_matrix.row(j).start(j)
-                                                  * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
+                                                  * _temporary.start(j)).coeff(0,0));
     m_matrix.coeffRef(j,j) = Djj;
 
     // Finish early if the matrix is not full rank or is indefinite.  This
@@ -194,12 +195,10 @@
     int endSize = size - j - 1;
     if (endSize > 0) {
       _temporary.end(endSize) = ( m_matrix.block(j+1,0, endSize, j)
-                                * m_matrix.col(j).start(j).conjugate() ).lazy();
+                                * _temporary.start(j) ).lazy();
 
-      m_matrix.row(j).end(endSize) = m_matrix.row(j).end(endSize).conjugate()
-                                   - _temporary.end(endSize).transpose();
-
-      m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / Djj;
+      m_matrix.col(j).end(endSize) = ( m_matrix.col(j).end(endSize) 
+				       - _temporary.end(endSize)    ) / Djj;
     }
   }