Re: [eigen] [RFC] Full pivoting LDLT 
[ Thread Index 
Date Index
 More lists.tuxfamily.org/eigen Archives
]
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 semidefinite.
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(size1) = m_matrix.row(0).end(size1) / m_matrix.coeff(0,0);
+ m_matrix.col(0).end(size1) /= 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;
}
}