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?
$\begingroup$
$\endgroup$
4
2 Answers
$\begingroup$
$\endgroup$
2
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.
answered May 11, 2019 at 15:14
-
$\begingroup$ Could you please elaborate or give references for the bullet with Hessenberg reduction plus specialized solver for systems $A+\sigma I$? $\endgroup$– VorKirCommented May 11, 2019 at 16:05
-
$\begingroup$ @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. $\endgroup$ Commented May 11, 2019 at 16:13
$\begingroup$
$\endgroup$
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.
See also: https://scicomp.stackexchange.com/questions/1001/how-does-the-matlab-backslash-operator-solve-ax-b-for-square-matrices for
A \ borscipy.linalg.solveorDGESVor whatever your programming language offers. $\endgroup$