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

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

2009/1/28 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>:
> On Wed, 28 Jan 2009, Benoit Jacob wrote:
>
>>> 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.
>
> How did you arrive at this conclusion?

When I implemented LU, I was initially hoping that partial pivoting
would be enough. What I say here is the results of my experiments. So
maybe I got something wrong, but aside from that possibility, I was
talking first-hand experience when saying that.

>I'm by no means an expert in
> numerical linear algebra, but the books I have (including Golub & Van Loan,
> Matrix Computations, 3e and Higham, Accuracy and Stability ..., 2e) say that
> partial pivoting should work fine.

At least Golub & van Loan do mention the numerical stability problem
of partial pivoting...
(plus I don't consider this book to be a solid reference on numerical
stability, e.g. with Keir we found a numerical stability issue in
their householder vector computation).

>> I think that LAPACK and GSL do a hybrid of partial and full pivoting
>> that I have yet to understand.
>
> That's not what their documentation says:
> http://netlib.org/lapack/lug/node38.html
> http://www.gnu.org/software/gsl/manual/html_node/LU-Decomposition.html

Ah OK. Here I was just repeating stuff I heared over IRC. So I don't
know. Perhaps they still do kind of a hybrid but don't mention it in
the docs. If they really do plain partial pivoting, then let's try
inverting, say, a 2000x2000 matrix, it will fail. For this reason I
doubt that they actually do plain partial pivoting.

> I think a 50% speedup is a worthwhile target. At least for vector solves,
> you can compute the residual and redo the computation with complete pivoting
> if necessary.

Here we're talking about a cubic complexity algorithm. So while a 50%
speedup sounds big, it just means that in the same time you handle a
matrix with 14% more rows and columns ( 1.5^(1/3) = 1.14). I want
Eigen to be easy to use and that implies prefering safeness over
performance.

Now let me retry my experiment with partial pivoting to see if it's
really unstable...

Cheers,
Benoit