What guidelines should I use when searching for good preconditioning methods for a specific problem?

Question

For the solution of large linear systems $Ax=b$ using iterative methods, it is often of interest to introduce preconditioning, e.g. solve instead $M^{-1}(Ax=b)$, where $M$ is here used for left-preconditioning of the system. Typically, we should have that $M^{-1}\approx A^{-1}$ and provide the basis for (much more) efficient solution or reduction in computational resources (e.g. memory storage) in comparison with solution of the original system (i.e. when $M=A$). However, what guidelines should we use to choose the preconditioner? How do practioneers do this for their specific problem?

Answer 1

I originally didn't want to give an answer because this deserves a very long treatment, and hopefully someone else will still give it. However, I can certainly give a very brief overview of the recommended approach:

1. Perform a thorough literature search.
2. If that fails, try every preconditioner that makes sense that you can get your hands on. MATLAB, PETSc, and Trilinos are nice environments for this.
3. If that fails, you should think carefully about the physics of your problem and see if it is possible to come up with a cheap approximate solution, perhaps even to a slightly changed version of your problem.

Examples of 3 are shifted Laplacian versions of Helmholtz, and Jinchao Xu's recent work on preconditioning the biharmonic operator with Laplacians.

Answer 2

Others have already commented on the issue of preconditioning what I will call "monolithic" matrices, i.e. for example the discretized form of a scalar equation such as the Laplace equation, the Helmholtz equation or, if you want to generalize it, the vector-valued elasticity equation. For these things, it is clear that multigrid (either algebraic or geometric) is the winner if the equation is elliptic, and for other equations it isn't quite so clear -- but something like SSOR often works reasonably well (for some meaning of "reasonable").

To me, the big revelation has been what to do about problems that are not monolithic, for example for the Stokes operator $$ \begin{pmatrix} A & B \\ B^T & 0 \end{pmatrix}. $$ When I started with numerical analysis some 15 years ago, I think people had the hope that the same techniques could be applied to such matrices as above, and the direction of research was to either try multigrid directly or to use generalizations of SSOR (using "point smoothers" like Vanka) and similar methods. But this has faded since it doesn't work very well.

What has come to replace this was what was initially called "physics-based preconditioners" and later simply (and maybe more accurately) "block preconditioners" like the one by Silvester and Wathen. These are often based on block eliminations or Schur complements and the idea is to build a preconditioner in such a way that one can reuse preconditioners for individual blocks that are known to work well. In the case of the Stokes equation, for example, the Silvester/Wathen preconditioner uses that the matrix $$ \begin{pmatrix} A & B \\ 0 & B^T A^{-1} B \end{pmatrix}^{-1} $$ when used as a preconditioner with GMRES would result in convergence in exactly two iterations. Since it is triangular, the inversion is also much simpler, but we still have the problem of what to do with the diagonal blocks, and there one uses approximations: $$ \begin{pmatrix} \widetilde{A^{-1}} & B \\ 0 & \widetilde{(B^T A^{-1} B)^{-1}} \end{pmatrix} $$ where the tilde means to replace the exact inverse by an approximation. This is often much simpler: because the $A$ block is an elliptic operator, $\widetilde{A^{-1}}$ is well approximated by a multigrid V-cycle, for example, and it turns out that here, $\widetilde{(B^T A^{-1} B)^{-1}}$ is well approximated by an ILU of a mass matrix.

This idea of working with the individual blocks that comprise the matrix and re-using preconditioners on individual ones has proven to be enormously powerful and has completely changed how we think of preconditioning systems of equations today. Of course, this is relevant because most actual problems are, in fact, systems of equations.

Answer 3

Jack has given a good procedure for finding a preconditioner. I will try an address the question, "What makes a good preconditioner?". The operational definition is:

A Preconditioner M accelerates the iterative solution of $A x = b$, and $M^{-1}$ can be applied cheaply compared to $A^{-1}$.

however this does not give us any insight into designing a preconditioner. Most preconditioners are based upon manipulation of the operator spectrum. Generically, Krylov methods converge faster when eigenvalues are clustered, see Matrix Iterations or Meromorphic Functions and Linear Algebra. Sometimes we can prove a preconditioner results is just a few unique eigenvalues, e.g. A Note on Preconditioning for Indefinite Linear Systems.

A common strategy is exemplified by Multigrid. Relaxational preconditioners (here smoothers) like SOR remove high frequency components in the error. When the residual is projected onto a coarse grid, lower frequency error components become higher frequency and can again be attacked by SOR. This basic strategy underlies more complicated versions of MG, such as AMG. Note that at bottom, the solver must resolve the lowest frequencies in the error accurately.

Another strategy involves solving the equation in small subspaces, which is exactly what Krylov solvers are doing. In the simplest form, this is the Kaczmarz method or the Additive Schwarz Method. The advanced strain of theory here, Domain Decomposition, concentrates on spectral approximation of the error on the interface, since the domains are assumed to be solved fairly accurately.

Then there are a bunch of things produce operators close to $A$ is some sense, which are easy to invert, but give no guarantees of actually improving convergence of the Krylov iteration, as far as I know. These preconditioners include incomplete factorization (ICC, ILU), sparse approximate inverse (SPAI), and some low rank approximations.