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

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

2009/1/28 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>:
> On Tue, 27 Jan 2009, Keir Mierle wrote:
>
>> Probably it's better to do full pivoting. Apparently cholesky is
>> stable for semidefinite matrices when full pivoting is used:
>>
>> http://eprints.ma.man.ac.uk/1101/01/covered/MIMS_ep2008_56.pdf
>
> Cholesky (LL^T) is stable without any pivoting for positive definite
> matrices. So I think the possibility of Cholesky without pivoting should
> remain.

Here is the thing that's got me thinking since yesterday and why I
didn't reply so far.
- on the one hand, since pivoting is essential for the stability of
LU, and LLt is a close variant of LU, I would assume that pivoting is
also important for the stability of LLt.
- on the other hand, something very special happens with LLt: the
matrix L is _unique_. Since it is well-determined mathematically, one
would expect it to be computable in a stable way. Now, if we do
pivoting, it's no longer a LLt decomposition but rather a (PL)(PL)t
decomposition with P a permutation matrix. So implicitly if the L of
LLt is computable in a stable way (and I believe it is), it has to be
doable without pivoting. This argument doesn't apply to LU as, since
the L and U aren't unique, one can't talk of "the L of LU" being
computable.

Just to say that I'm puzzled by LLt and I don't know whether or not
pivoting is needed here (by contrast at least for LU and SVD pivoting
is clearly needed).

> Why do you need an epsilon? Perhaps a better question is, why does the
> routine test whether the matrix is positive definite? As you say, it's not
> easy to write such a test, and it could be said that it's the user's
> responsibility that the matrix is positive definite - though of course a
> test is nice when the code is run in debug mode.

I agree. If an algorithm works only for positive definite matrices, it
has to be the responsibility of the user.
Moreover we should not check for that. Just like the compiler doesn't
check for x!=0 when we do 1/x.
Also, here by "positive definite" we have to mean "the ratio of the
biggest to the lowest eigenvalue is not too big" or if you prefer "the
lowest eigenvalue is not negligible when compared to the biggest
eigenvalue" i.e. the notion of positive definite has to be "fuzzy".

> A related issue I've been wondering about is: why does LU decomposition use
> complete pivoting? I believe that's quite a bit more expensive than partial
> pivoting

It's about 2x more expensive currently, as the pivot lookup part is
not yet vectorized. If it were, it would still be 1.5x more expensive.

> When solving a linear system with
> an invertible matrix, partial pivoting performs as well as complete
> pivoting, doesn't it? LAPACK uses partial pivoting, as does GSL.

No, partial pivoting fails even for invertible matrices of size
200x200, and fails consistently at sizes like 1000x1000.

I think that LAPACK and GSL do a hybrid of partial and full pivoting
that I have yet to understand. In any case, anything less than full
pivoting will fail in certain cases. I believe that LAPACK's
implementation will only fail on specially tailored matrices, but
still, if I understand how their hybrid pivoting works, I can
construct an invertible matrix making it fail, even without specially
weird coefficients. So I'm not comfortable with this, all this
insecurity for, what, at 50% speedup. Plus, a user on the forum
reports that our current implementation is already as fast as LAPACK
for size 1500x1500.

Cheers,
Benoit