# For which problems Krylov subspace methods are preferred over multigrid methods?

As multigrid methods are known to have grid independent convergence rates with $O(N)$ computational cost, then why would one be interested in using Krylov subspace methods at all, for which convergence rates deteriorate with grid refinement? Are there any problems where Krylov solvers perform better than MG methods?

• To date, multigrid methods have not been (and maybe can't even be) applied to matrix functions. Jun 15 '17 at 19:49
• @artificial_moonlet Could you elaborate on your comment? I was under the impression that multigrid methods can be used as preconditioners for Krylov based solvers (as noted in an answer by @JesseChan). Do you disagree? Jun 16 '17 at 3:38
• How does preconditioning work for $f(A)b$? In general, this is an open problem. You can certainly approximate $f(A)$ with quadrature or some other rational approximation and then apply multigrid on each resulting linear system, but such techniques only work for certain classes of matrix-function pairs. Jun 20 '17 at 12:17

The answer depends somewhat on the discretization; for example, some boundary integral discretizations result in very well-conditioned matrices, for which a Krylov solver works just fine without having to introduce the multi-level machinery of multigrid.

For ill-conditioned matrices arising from finite difference or finite element methods, Krylov solvers are almost always paired with a preconditioner, for which multigrid methods can work very well.

• An extra nice feature of Krylov methods is that you can get information about the spectrum of the system, e.g. max/min eigenvalues. You can then leverage this knowledge for multigrid. For example, if you use Chebyshev smoothers, you want to have some advance knowledge of where the spectrum of the operator lives in the complex plane. Krylov methods and multigrid methods can both benefit from each other. Jun 19 '17 at 13:49

This question is pretty well discussed in literature. However, there are lots of questions concerning multigrid on SciComp, so I decided to compose more or less detailed answer.

# I. When multigrid DOES work fine as a stand–alone solver

MG works fine for rather simple elliptic problems. In fact, both MG and PCG (with MG cycles as inner iterations) usually give similar results for such problems in terms of iterations / working time and stay robust as you refine the mesh.

I solved diffusion–reaction 2D problem $$-\epsilon \, \Delta u + u = f, \quad \mathbf x \in \Omega, \\ u|_{\Gamma_D} = g_D, \\ \left( \epsilon \, \nabla u \cdot \hat{\mathbf n} + u \right)|_{\Gamma_R} = g_R$$

with CATS’ PDEs on hierarchy of triangular meshes $\Omega_1 \subset \Omega_2 \subset \dots \subset \Omega_8$ of unit square with P1 FEs. I obtained the following results:

• $\text{MG(V)}$ and $\text{MG(W)}$ are V– / W–cycles as stand–alone solvers; $\mathbf P_{\text{MG(V)}}\text{CG}$ is preconditioned conjugate gradient with one V–cycle as an inner iteration; $\mathbf P_{\text{ILU(0)}}\text{CG}$ is preconditioned conjugate with incomplete $(\mathbf L + \mathbf I) \, \mathbf D \, (\mathbf I + \mathbf L^T)$ preconditioner (zero fill–in).
• I used 2 pre– and post–smoothing SSOR iterations inside MG–cycles.
• Numb of DOFs on the coarsest grid: 17; numb of DOFs on the finest grid: 656 897.
• In the table: “numb of iterations, $m$” / “solving time (in seconds), $t$.” I always started with a random initial guess $\mathbf x_0$ in order to introduce components of different frequencies in the initial error vector $(\mathbf x - \mathbf x_0)$; I stopped solving when $||\mathbf r_m||/||\mathbf r_0|| < 10^{-12}$.

Here’s comparison in terms of time $t$ as system size $n$ grows:

Indeed, you can see that $\text{MG}$ allows you to solve this problem in $O(n)$ time; so does $\mathbf P_{\text{MG}}\text{CG}$, but there is no difference between them in this case. I think this makes people ask questions like

As multigrid methods are known to have grid independent convergence rates with $O(n)$ computational cost, then why would one be interested in using Krylov subspace methods at all […]?

# II. When multigrid DOES NOT work fine as a stand–alone solver

Consider the following cases.

## II.1 More complex elliptic problems

It is usual to have jumping coefficients. For one, consider the following magnetostatic Poisson problem $$\nabla \cdot \left( \frac{1}{\hat{\mu}} \, \nabla A_z \right) = \mu_0 \, J_z.$$

It is natural to have different values for magnetic permeability $\hat{\mu}$ for a magnet (gray box on the picture above) and air; current density $J_z$ is only nonzero within wires (red and blue squares on the picture above).

MG as a stand–alone solver loses its robustness (or even diverges) with respect to jump sizes for such problems, although the configuration is very simple. However, it still works nice as a preconditioner for Krylov solvers. If I have some free time, I will add a similar table for this problem.

There are lots of simple elliptic problems for which MG does not work. It is useful to look through a paper I mentioned in this answer in order to get a deeper understanding.

## II.2 “Real–life” problems

Complex real–life problems involve different physics, several unknown fields and so forth. From the algebraic point of view it usually means block structure of the resulting (non-)linear systems. (Consider, for one, Stokes or Oseen problems, fluid–structure–interaction models typical for haemodynamics, thermoelasticity etc.)

It turns out that MG is in general a poor choice for such kind of problems. However, these problems usually “contain” some kind of “diffusion” (read: elliptic / Laplace blocks in the matrix), and one usually keeps this in mind while constructing preconditioners. Indeed, one often uses some kind of MG–cycles to approximate actions of inverses of elliptic blocks. One example I mentioned in this answer.

In order to get a better understanding, it is useful to check this lecture by @wolfgang-bangerth.

• Is MG generally unreliable for problems with large jumps in coefficients? Or does it depend on the discretization/coarsening strategies? (I am thinking of FV solvers with operator interpolation, like Black-box Multigrid, and transmissibility upscaling in in oil-reservoir simulation.) May 19 '18 at 18:24