LU decomposition

I found the following page that I think is a very good step-by-step example of how to do LU decomposition

LU Decomposition

Beware that LU factorisation on is a very poor method of solving general linear systems with floating-point data. By all means learn it, but don't use it (unmodified) for anything that matters, unless you want unpredictable and potentially catastrophic error accumulation.

Quote:

Original post by TheAdmiral Beware that LU factorisation on is a very poor method of solving general linear systems with floating-point data. By all means learn it, but don't use it (unmodified) for anything that matters, unless you want unpredictable and potentially catastrophic error accumulation.

The Cholesky decomposition makes use of the LU factorization. This is used in fitting curves to samples (for example, to use a B-spline curve to fit keyframe animation data as a means of data compression). I would certainly not call the use of LU factorization here "poor"...

Admiral,

What would you use instead?

Quote:

Original post by Dave Eberly The Cholesky decomposition makes use of the LU factorization. This is used in fitting curves to samples (for example, to use a B-spline curve to fit keyframe animation data as a means of data compression). I would certainly not call the use of LU factorization here "poor"...

You're right, sorry - I didn't qualify this claim fully. I meant that using LU factorisation and back-substitution to solve a general system of linear equations is a bad idea. In the symmetric positive-definite case (as in a Cholesky factorisation), numerical instability is rarely a problem.

I felt that this needed bringing up since LU factorisation is the first way everyone learns to solve a linear problem (as it's so simple), but not everybody is told that it is highly unstable even for many realistic problems.

Quote:

Original post by assaf What would you use instead?

This depends on your input data and resources. For 'small', dense systems (say up to 1000x1000) LU-factorisation with partial-pivoting, or QR-factorsation are excellent O(n³) choices. O(n²) alternatives exist, but they have much higher coefficients and tend not to be beneficial until the input volume gets very large. These also tend to take a lot more work to implement. GMRES is the favoured method for arbitrary large systems.

Do you have any other constraints? Tell us how large you expect the system to be, what guarantees you do or don't have on its non-singularity, and if the matrix takes on any special forms (block, banded or sparse in particular). Better yet, assuming it's not too deep, tell us the origins of the system needing to be solved.

Quote:

Original post by TheAdmiral You're right, sorry - I didn't qualify this claim fully. I meant that using LU factorisation and back-substitution to solve a general system of linear equations is a bad idea. In the symmetric positive-definite case (as in a Cholesky factorisation), numerical instability is rarely a problem. I felt that this needed bringing up since LU factorisation is the first way everyone learns to solve a linear problem (as it's so simple), but not everybody is told that it is highly unstable even for many realistic problems.

I agree with you that LU decomposition is (potentially very) unstable if you don't use any pivoting (you can not do that on any matrix though - it has to be a particularlly nice matrix for you to be able to avoid at least partial pivoting anyway :). Now, correct me if I'm wrong, but I was under the impression that LU decomposition is numerically stable with full pivoting (can't remember if partial pivoting is enough for stability here..).

Quote:

Original post by Devilogic Now, correct me if I'm wrong, but I was under the impression that LU decomposition is numerically stable with full pivoting (can't remember if partial pivoting is enough for stability here..).

It's true that situations exist where partial pivoting is insufficient to guarantee stable convergence, but in reality this is almost never a problem. Full pivoting does improve stability (over its partial counterpart) in general, but this improvement is very small considering the extra work required at each pass. For this reason, if stability is essential then shifted LU-factorisation with partial-pivoting is usually employed, as it's cheaper than full-pivoting without shifts (and I believe it to also be maximally stable, though I can't prove this).

In certain situations an 'intelligent' pivoting algorithm may work on both rows and columns so as to preserve sparse- or block-structure. As far as I know, this is the extent of its usefulness.

Thank you for the explanation TheAdmiral :)

Ummm...wow, big words. I just wanted to find the inverse of a 4x4 matrix. Its okay though, I found a solution. I decided to use the determinant method. I've learnt quite a bit though, thanks all.