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

[Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>]

On Sun, 28 Feb 2010, Ben Goodrich wrote:

> Right. And I am basically just trying to figure out the fastest way to
> get a consistent D that has decent numerical precision even when
> A(theta) is nearly singular (but not indefinite). If I am
> understanding correctly, I need to unpivot in basically the same way
> that reconstructMatrix() does. But unpivoting is going to be
> computationally expensive when it has to be done repeatedly in an
> optimization context, so I was hoping someone had a short-cut relative
> to the steps that first came to my mind :) .

I'm still not sure I understand what you want. The decomposition computed 
by Eigen is A = P^T L D L^T P, while you want to compute the decomposition 
A = L D L^T. My first thought would be to compute this from the Cholesky 
decomposition. Suppose A = X X^T is the Cholesky decomposition (with X 
lower triangular). Let D' to be the diagonal matrix with the same entries 
on the diagonal. Set D = (D')^2 and L = X (D')^(-1). Then L is triangular 
and D is diagonal, so A = L D L^T is the decomposition you're looking for.


Cheers,
Jitse