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

[Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>]

On Sun, Feb 28, 2010 at 9:48 AM, Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx> wrote:
> 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.

When it comes right down to it, I pretty much just want mutually
exclusive things. At the optimum, A(theta) will be positive definite
but singular. But I would like maximum precision, so the square root
free Cholesky is nice. Full pivoting would be nice too but that would
create discontinuities in f(theta) unless I unpivoted. But unpivoting
the decomposition would probably require more operations than the
decomposition itself, and I would like speed too. So, no free lunch
for me (again). With 500x500 random matrices that have rank 250, I was
getting zeros that were on the order of 1e-12 in the worst case, which
was better than I would have expected without pivoting.

And it would probably be good on principle to have a nonpivoting
option in the eigen3 LDLT just so it could be somewhat backward
compatible with eigen2. So, I'll send a cleaned up patch RSN.

Thanks,
Ben