[Math] Methods of solving linear system of equations, how to select the appropriate method

gaussian eliminationlinear algebrana.numerical-analysisnumerical linear algebra

A linear system of equations Ax=b can be solved using various methods, namely, inverse method, Gauss/Gauss-Jordan elimination, LU factorization, EVD (Eigenvalue Decomposition), and SVD (Singular Value Decomposition).
I know that there are several disadvantages of using inverse method; for example, with ill conditioned matrix A the solution can not be computed with inverse method. Moreover, I know that with changing vector b LU factorization has advantage over Gauss/Gauss-Jordan elimination.
How to decide between LU, SVD, and EVD?
Is there any scenario where Gauss/Gauss-Jordan elimination has advantage over LU, SVD, and EVD?

Best Answer

Disclaimer 1: Treating these topics properly would require a quick course in numerical analysis.

Disclaimer 2: If you are using any sane computer system, it's already going to have a library function to solve linear systems implemented, which is going to be better of what you can code yourself if you don't have a solid grasp of numerical linear algebra and scientific computing. Unless you have special requirements, just use A \ b or scipy.linalg.solve or DGESV or whatever your programming language offers.

Disclaimer 3: What follows applies to solving dense, small-scale (let's say $n\leq 10000$) square linear systems in IEEE floating point arithmetic. Things may be different if you use integer / rational / symbolic arithmetic in a computer algebra system, or if you compute by hand.

That said:

  • LU factorization (with partial pivoting): industry standard. It has the drawback that on some examples (very unlikely, basically just happens in counterexamples), the entries of $U$ may grow as $\sim 2^n$ times the entries of $A$, and thus it may become unstable. Usually not a concern in practice.

  • QR factorization: costs 2x as LU factorization, but avoids the problem mentioned above. If you don't mind paying 2x the cost just for the added peace of mind, go for it.

  • Gaussian elimination: basically, it's the LU factorization method, but you avoid storing L in memory (so you can't reuse the factorization in case you have multiple right-hand sides $b$ coming at different times). People usually just go with LU, because the few memory writes saved aren't typically worth the trouble.

  • Gaussian elimination without partial pivoting: avoid. Can be ridiculously unstable.

  • inverse method: not super-clear what you mean here since there are different ways to compute a matrix inverse; in any case, avoid. It is more expensive ($2n^3$ rather than $4/3n^3$ flops) and has worse stability properties (unstable in cases in which $\|b\| \ll \|A\|\|x\|$).

  • Gauss-Jordan elimination: Avoid. No advantages over LU, and has worse stability properties (it's forward stable but typically not backward stable).

  • Cholesky and LDL factorization (with suitable pivoting schemes, e.g. Bunch-Kaufmann): specialized solvers for symmetric matrices. Use them over LU if you have symmetric matrices; they are cheaper by a factor 2.

  • Eigenvalue decomposition / Jordan decomposition: Avoid. Can be ridiculously unstable.

  • Schur factorization: Avoid. Way more expensive; if you just have to solve a linear system it's not worth the trouble.

  • Hessenberg reduction + specialized Hessenberg solver: might be useful in the special case in which you have to solve systems with $A+\sigma_i I$ for many values of $\sigma_i$, since in this case you can re-use the same factorization to solve all the systems. Otherwise, avoid; slower than LU and QR without other advantages.

  • Singular value decomposition: Avoid. Way more expensive; if you just have to solve a square linear system it's not worth the trouble. If you have a least-squares problem, well, it's a whole other story.

Note that a performance-aware implementation of any of these methods will require using specialized libraries for the linear algebra primitives (e.g., matrix-vector products or rank-1 updates) and implementing in-place factorizations, compressed representations of orthogonal matrices (for algorithms in which they appear), and (most importantly) blocked versions of these algorithms.

TL;DR: use LU factorization for the best speed/stability tradeoff. Don't roll your own code if performance matters.

Reference: Higham, Accuracy and stability of numerical algorithms; Golub, Van Loan, Matrix Computations, 4th ed.

Related Question