| Re: [eigen] Instability in LLT and LDLT methods. |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
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? 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.
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
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.
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.
Cheers,
Jitse