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

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?

• Even for one particular class of equations, this would require a very long and detailed answer... – Jack Poulson Feb 1 '12 at 21:42
• Should be possible to suggest heuristic strategies for how preconditioners can be chosen. For example, given a problem, what do practioners do in practice to try and find a good preconditioner? Just start with a basic diagonal preconditioner based on extracting the diagonal from $A$? or? – Allan P. Engsig-Karup Feb 1 '12 at 21:52
• I'm going to channel @MattKnepley and say that the appropriate action is to do a literature search. If that fails, try all easily available options on a reasonably large problem. If that fails, think deeply about the physics and how you can come up with an approximate solution to the problem cheaply, and use that as the preconditioner. – Jack Poulson Feb 1 '12 at 21:54
• @JackPoulson: Since this question is in a similar vein to Which sparse linear system solver to use? and How to choose a scalable linear solver, it seems on topic to me (although broad). Since your comment is basically an answer, could you please convert it to an answer? – Geoff Oxberry Feb 1 '12 at 23:08
• I've started a bounty on this question, but I'm also interested in seeing more questions in this vein that might be better-posed or restricted to a specific class of problems. – Aron Ahmadia Apr 3 '12 at 10:34

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.

• Thanks! The rest of this comment meets the minimum character limit. – Geoff Oxberry Feb 2 '12 at 6:24

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.

• Man, yeah, I so wanted the bounty! ;-) – Wolfgang Bangerth Apr 10 '12 at 20:00
• In your second paragraph: "But this has faded since it doesn't work very well." Can you give some intuition about why it doesn't work very well? Are there circumstances under which it may work? – Andrew T. Barker Apr 11 '12 at 18:30
• The reason why direct multigrid applied to whole systems hasn't proven to be so successful is that the smoother needs to conserve the structural properties of the equation, and that's nontrivial to achieve. For example, if you want to apply multigrid to the Stokes equations, you have to have a smoother that given a divergence free vector gives you a divergence free vector. There are such smoothers for Stokes, but it's nontrivial to construct and it usually takes away from the quality as a smoother/solver. It becomes much more difficult to preserve properties in more egneral cases. – Wolfgang Bangerth Apr 14 '12 at 17:16
• As for generalizing things like Jacobi/SSOR/etc to systems: most of these methods require that the diagonal entries of the matrix are nonzero. That's obviously not the case for Stokes. So the next simplest method is to not look at individual matrix rows but blocks of rows, e.g., all rows for DoFs associated with a single vertex. These are call "point-smoothers" (point as in vertex) and they work to a degree, but they suffer from the same degradation of performance as Jacobi/SSOR once problems become large.To avoid that, a preconditioner needs to globally exchange information like multigrid. – Wolfgang Bangerth Apr 14 '12 at 17:19
• Multigrid is famously ineffective at solving Helmholtz because, mostly because the low energy oscillatory modes are difficult to smooth or to represent in a coarse space. There has been some work on wave-ray multigrid, but the formulation is very technical and it's not a mature methodology at this point. Note that non-symmetric systems can also be solved using this sort of block decomposition. Depending on the choice of variables, (e.g. primitive vs. conservative), a change of basis may be required inside the preconditioner to expose the blocked structure. – Jed Brown Apr 20 '12 at 15:40

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.

• thank you for your reply. Any experiences regarding how far we should go to actually make a proof that the preconditioning works for large systems -- and possibly how this can or should be done in practice. It is my experience that for many systems we have to rely on intuition, heuristics, etc. – Allan P. Engsig-Karup Apr 4 '12 at 11:30
• I think intuition is going too far. What I see in practice is a proof for a simple system. Then an argument that some modification should be insensitive to a parameter, or a certain kind of variation. Then numerical experiments showing it works even outside this model of variation. – Matt Knepley Apr 4 '12 at 13:01