How to choose a method for solving linear equations

Question

To my knowledge, there are 4 ways to solving a system of linear equations (correct me if there are more):

1. If the system matrix is a full-rank square matrix, you can use Cramer’s Rule;
2. Compute the inverse or the pseudoinverse of the system matrix;
3. Use matrix decomposition methods (Gaussian or Gauss-Jordan elimination is considered as LU decomposition);
4. Use iterative methods, such as the conjugate gradient method.

In fact, you almost never want to solving the equations by using Cramer's rule or computing the inverse or pseudoinverse, especially for high dimensional matrices, so the first question is when to use decomposition methods and iterative methods, respectively. I guess it depends on the size and properties of the system matrix.

The second question is, to your knowledge, what kind of decomposition methods or iterative methods are most suitable for certain system matrix in terms of numerical stability and efficiency.

For example, the conjugate gradient method is used to solve equations where the matrix is symmetric and positive definite, although it can also be applied to any linear equations by converting $\mathbf{A}x=b$ to $\mathbf{A}^{\rm T}\mathbf{A}x=\mathbf{A}^{\rm T}b$. Also for positive definite matrix, you can use Cholesky decomposition method to seek the solution. But I don't know when to choose the CG method and when to choose Cholesky decomposition. My feeling is we'd better use CG method for large matrices.

For rectangular matrices, we can either use QR decomposition or SVD, but again I don't know how to choose one of them.

For other matrices, I don't now how to choose the appropriate solver, such Hermitian/symmetric matrices, sparse matrices, band matrices etc.

Answer 1

Your question is a bit like asking for which screwdriver to choose depending on the drive (slot, Phillips, Torx, ...): Besides there being too many, the choice also depends on whether you want to just tighten one screw or assemble a whole set of library shelves. Nevertheless, in partial answer to your question, here are some of the issues you should keep in mind when choosing a method for solving the linear system $Ax=b$. I will also restrict myself to invertible matrices; the cases of over- or underdetermined systems are a different matter and should really be separate questions.

As you rightly noted, option 1 and 2 are right out: Computing and applying the inverse matrix is a tremendously bad idea, since it is much more expensive and often numerically less stable than applying one of the other algorithms. That leaves you with the choice between direct and iterative methods. The first thing to consider is not the matrix $A$, but what you expect from the numerical solution $\tilde x$:

1. How accurate does it have to be? Does $\tilde x$ have to solve the system up to machine precision, or are you satisfied with $\tilde x$ satisfying (say) $\|\tilde x - x^*\| < 10^{-3}$, where $x^*$ is the exact solution?
2. How fast do you need it? The only relevant metric here is clock time on your machine - a method which scales perfectly on a huge cluster might not be the best choice if you don't have one of those, but you do have one of those shiny new Tesla cards.

As there's no such thing as a free lunch, you usually have to decide on a trade-off between the two. After that, you start looking at the matrix $A$ (and your hardware) to decide on a good method (or rather, the method for which you can find a good implementation). (Note how I avoided writing "best" here...) The most relevant properties here are

• The structure: Is $A$ symmetric? Is it dense or sparse? Banded?
• The eigenvalues: Are they all positive (i.e., is $A$ positive definite)? Are they clustered? Do some of them have very small or very large magnitude?

With this in mind, you then have to trawl the (huge) literature and evaluate the different methods you find for your specific problem. Here are some general remarks:

• If you really need (close to) machine precision for your solution, or if your matrix is small (say, up to $1000$ rows), it is hard to beat direct methods, especially for dense systems (since in this case, every matrix multiplication will be $\mathcal{O}(n^2)$, and if you need a lot of iterations, this might not be far from the $\mathcal{O}(n^3)$ a direct method needs). Also, LU decomposition (with pivoting) works for any invertible matrix, as opposed to most iterative methods. (Of course, if $A$ is symmetric and positive definite, you'd use Cholesky.)

  This is also true for (large) sparse matrices if you don't run into memory problems: Sparse matrices in general do not have a sparse LU decomposition, and if the factors do not fit into (fast) memory, these methods becomes unusable.

  In addition, direct methods have been around for a long time, and very high quality software exists (e.g., UMFPACK, MUMPS, SuperLU for sparse matrices) which can automatically exploit the band structure of $A$.

• If you need less accuracy, or cannot use direct methods, choose a Krylov method (e.g., CG if $A$ is symmetric positive definite, GMRES or BiCGStab if not) instead of a stationary method (such as Jacobi or Gauss-Seidel): These usually work much better, since their convergence is not determined by the spectral radius of $A$ but by (the square root) of the condition number and does not depend on the structure of the matrix. However, to get really good performance from a Krylov method, you need to choose a good preconditioner for your matrix - and that is more a craft than a science...

• If you repeatedly need to solve linear systems with the same matrix and different right hand sides, direct methods can still be faster than iterative methods since you only need to compute the decomposition once. (This assumes sequential solution; if you have all the right hand sides at the same time, you can use block Krylov methods.)

Of course, these are just very rough guidelines: For any of the above statements, there likely exists a matrix for which the converse is true...

Since you asked for references in the comments, here are some textbooks and review papers to get you started. (Neither of these - nor the set - is comprehensive; this question is much too broad, and depends too much on your particular problem.)

• Golub, van Loan: Matrix Computations (still the classical reference on matrix algorithms; the much expanded fourth edition now also discusses sparse matrices and has Matlab code in place of Fortran as well as an extensive bibliography)
• Davis: Direct Methods for Sparse Linear Systems (a good introduction on decomposition methods for sparse matrices)
• Duff: Direct Methods (review paper; more details on modern "multifrontal" direct methods for sparse matrices)
• Saad: Iterative methods for sparse linear systems (the theory and - to a lesser extent - practice of Krylov methods; also covers preconditioning)

Answer 2

The decision tree in Section 4 of the relevant chapter in the NAG Library Manual answers (in part) some of your questions.

Answer 3

The Eigen Library Documentation also has a nice overview page with a lot of information about different matrix decompositions:

http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html