# How to choose a method for solving linear equations

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.

• Hi @chaohuang, and welcome to SciComp! You may want to see this discussion: scicomp.stackexchange.com/questions/81/…
– Paul
Sep 11, 2012 at 16:37
• Your question has massive scope and may be a bit too broad for the Q&A format that we have here on the stackexchange... is there a particular class of matrix system that you are interested in?
– Paul
Sep 11, 2012 at 22:41
• @chaohuang There are numerous books on this subject. This question is a bit like asking a medical doctor how they choose treatments "in general". If you want to ask a question that is not specific to a certain class of problems, you should put in the work to become familiar enough with the field to ask something precise. Otherwise, explain the specific problem that you are concerned with. Sep 12, 2012 at 2:52
• From the FAQ: If you can imagine an entire book that answers your question, you’re asking too much. There are entire journals, and hundreds of books, associated with this question. Sep 12, 2012 at 5:09
• Can you guys recommend a book comprehensively covering the question? Sep 12, 2012 at 5:18

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.)

• I like your analogy of the screwdriver!
– Paul
Sep 16, 2012 at 17:25
• @chaohuang If this answered your question, you should accept it. (If it didn't, feel free to point out what is missing.) Nov 29, 2012 at 13:05
• @ChristianClason accepted it. I was waiting and hoping someone could shed some light on the issue of rectangular matrices. Since it has been a long time, I guess there will never be such an answer :( Nov 29, 2012 at 16:06
• @chaohuang Thank you. If you're still interested in rectangular matrices, you should pose a (linked) question on "How to choose a method for solving overdetermined systems". Nov 29, 2012 at 16:10
• Here a reference on use of iterative methods for solving large sparse systems of linear equations. Jan 11, 2013 at 2:27

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

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