2
$\begingroup$

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?

$\endgroup$
  • 6
    $\begingroup$ 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. $\endgroup$ – Dirk May 11 at 8:06
  • 5
    $\begingroup$ 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... ? $\endgroup$ – Federico Poloni May 11 at 12:18
  • 3
    $\begingroup$ 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. $\endgroup$ – Federico Poloni May 11 at 12:21
  • $\begingroup$ @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. $\endgroup$ – MBK May 11 at 14:43
11
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Could you please elaborate or give references for the bullet with Hessenberg reduction plus specialized solver for systems $A+\sigma I$? $\endgroup$ – VorKir May 11 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$ – Federico Poloni May 11 at 16:13
2
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.