# Methods of solving linear system of equations, how to select the appropriate method

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?

• On top of this there are iterative methods like Gauss-Seidel or CG (which is actually both iterative and exact in some sense). I would suspect that asking on scicomp.se would be better.
– Dirk
May 11, 2019 at 8:06
• How are you solving these linear systems? By hand, on a computer with IEEE double-precision arithmetic, on a computer algebra system with exact computations... ? May 11, 2019 at 12:18
• 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. May 11, 2019 at 12:21
• @FedericoPoloni lets suppose I am creating my own library in C/C++. The question is not about tool the question is about knowing the advantage. Like LU factorization is clearly better choice over Gauss/Gauss-Jorden elimination if our system has dynamic values in vector b. May 11, 2019 at 14:43

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.

• Could you please elaborate or give references for the bullet with Hessenberg reduction plus specialized solver for systems $A+\sigma I$? May 11, 2019 at 16:05
• @VorKir If $A=QHQ^*$, with $Q$ orthogonal and $H$ Hessenberg (and this factorization is produced in $O(n^3)$ as an initial step of the classical eigenvalue solvers), then $A+\sigma I =Q(H+\sigma I)Q^*$ for all $\sigma$, and you can use this factorization to solve a sequence linear systems of the form $(A+\sigma_i I)x=b_i$ in $O(n^2)$ each, since both orthogonal and Hessenberg linear systems can be solved in $O(n^2)$. It's a simple observation, but I don't know of a specific reference treating it. May 11, 2019 at 16:13

If you want to see what a group of experts did, Matlab has a flowchart for their mldivide (aka "backslash") : https://www.mathworks.com/help/matlab/ref/mldivide.html

The punchline appears to be that unless the matrix has structure, use LU for square and QR for non-square.