A direction the smoother cannot see

Everything the multigrid field has measured is on the Laplacian, in which the coupling to a point’s neighbour on the left is the same as the coupling to its neighbour above. That is not a property of diffusion problems generally; it is a property of isotropic diffusion on a square grid, which is the case every component of the method was designed for.

Change one number. Take

−ε ∂²u/∂x² − ∂²u/∂y²

with ε a small positive constant, so that the material conducts a thousand times better vertically than horizontally. The operator is still symmetric. It is still positive definite. It still has five nonzeros a row, the same eigenvectors as before, and eigenvalues given by a formula in the two frequencies. Nothing about it looks like a hard problem.

At ε = 10⁻³ the V-cycle’s convergence factor goes from 0.2016 to 0.9565 — each run at the best weight its own closed form gives, so the comparison is between two methods rather than between a method and a badly tuned one — with the same smoother, the same transfer operators, the same number of levels and the same code.

The problem is not harder, and this is the part that matters #

The first thing to check is whether the anisotropic problem is simply a worse-conditioned one, in which case there would be nothing to explain — every method on this site slows down as κ grows and the essays on conditioning and Krylov rates say so at length.

It is not. The eigenvalues are

λ(m, l) = 4ε sin²(θx/2) + 4 sin²(θy/2)

largest at (k, k) and smallest at (1, 1), and at both of those extremes the two sine terms are equal, so both are (1 + ε) times the isotropic value, the factor cancels, and the ratio is

κ = cot²(πh/2)

with no ε in it anywhere. Measured across four decades of ε at k = 15, the condition number is 103.08686891981742 every time; the relative spread over the whole sweep is 2.8·10⁻¹⁶, which is rounding.

So the anisotropic problem is, as far as any bound stated in κ can see, the same problem. A reader handed the condition number and asked to predict the trouble would predict none. Everything that goes wrong goes wrong in the method, and specifically in one component of it, and none of it is visible in the one number anybody thinks to look at.

That is a stronger statement than “multigrid is sensitive to anisotropy”. It says the sensitivity is not a property of the equations at all.

Every sufficient condition still holds #

The second thing to check is whether the smoother has stopped satisfying the hypotheses it is usually justified under, because if it had, the failure would be a hypothesis violation rather than anything deeper.

Written out, the stencil is

⎡  0   −1    0 ⎤
⎢ −ε  2ε+2  −ε ⎥
⎣  0   −1    0 ⎦

and it passes every classical test. The diagonal is 2ε + 2 and the off-diagonal magnitudes sum to 2ε + 2, so the matrix is weakly diagonally dominant, exactly as the isotropic Laplacian is. Every off-diagonal entry is negative and every diagonal entry positive, so it is an M-matrix. It is symmetric to the bit and positive definite, with its smallest eigenvalue 4(1 + ε)sin²(πh/2) > 0. Weighted Jacobi on it converges for every ω in (0, 1], as it does on any symmetric positive definite matrix with a positive diagonal.

So the smoother is a convergent iteration on a well-conditioned symmetric positive definite M-matrix, and it is useless. Convergence is not the property the method needs from it; what the method needs is that the modes it fails to damp are a set the coarse grid can carry, and that is a joint property of two components which neither of them can be checked for alone.

That is worth stating in the site’s usual form. An assertion about a smoother is not an assertion about a method. The three-line justification “the iteration matrix has spectral radius below one, so the smoother converges” is true here, is checked here, and predicts nothing about the number this essay is about.

Which half broke #

The method is two claims joined by an assumption.

Smoothing works on the rough half. Weighted Jacobi damps the oscillatory modes and leaves the smooth ones nearly untouched — the measurement the field opened with.

The coarse grid represents the smooth half. A function that is smooth on a fine grid is a function a grid with half the points can carry, so the error the smoother leaves can be handled somewhere cheaper.

The assumption joining them is that the modes the smoother cannot damp are the modes the coarse grid can represent. Both claims can be true and the assumption false, and that is precisely what happens here.

Consider the mode that oscillates as fast as the grid allows across x and is nearly constant along y — θx = π, θy ≈ 0. Apply one weighted Jacobi sweep. The diagonal of the operator is 2ε + 2, and the part of the operator that acts on variation across x is 4ε sin²(θx/2), so the mode is multiplied by

1 − 2ωε/(1 + ε)

which at ε = 10⁻³ and ω = 2/3 is 0.99867. Three hundred sweeps would halve it. The smoother is not damping this mode in any useful sense, because the only coupling that can see it is the weak one.

And the coarse grid cannot help, because the mode is not smooth. It oscillates at the grid frequency in x, which is exactly what a coarser grid cannot carry.

The smoothing factor, in closed form #

The measurement is the smoothing factor: the largest damping over the frequencies the coarse grid cannot see. Scanning 128×128 frequency pairs gives it; so does an expression, and this phase found the expression by reading it off the scan rather than the other way round.

The symbol of the iteration is monotone in both sine terms, so its extremes over the region outside the coarse grid’s quarter are attained at that region’s corners, and only two of them matter — the least-damped mode and the mode at (π, π) where the damping overshoots to 1 − 2ω. So

μ = max( |1 − ωε/(1+ε)| , |1 − 2ω| )

and the scan agrees with it to 10⁻⁹ at every ε and ω tested. At ω = 2/3 that gives 2/3 when ε = 1, 0.9394 at ε = 0.1, 0.9934 at ε = 0.01 and 0.9993 at ε = 10⁻³.

The ε in the first branch is the whole failure. As ε → 0 the least-damped mode’s factor goes to one, and there is no weight that rescues it: raising ω to attack that branch lowers 1 − 2ω past −1 and the other branch takes over. The best available at ε = 10⁻³ is ω = 0.9995 and μ = 0.9990, which is not a repair, it is a rounding.

Both directions, and the symmetry that is not there #

There is a temptation to read all this as “the anisotropic problem is hard”, and one measurement refuses it flatly.

Turn the operator through a right angle — put the ε on ∂²u/∂y² instead, so the strong direction runs along x rather than y. On a square grid that is the same problem reflected in the grid diagonal, and the point smoother returns the same convergence factor to two decimal places, because it is aimed at neither direction and cannot tell them apart.

The two repairs the next essays build do not. Both go from 0.037 and 0.107 to 0.967 and 0.966, on a problem that differs from the one they solved only in which axis somebody chose to call y. That is the subject of the third essay in this sequence and it is worth flagging here, because it is the reason this failure is interesting rather than merely inconvenient: the obvious fixes encode a direction, and the direction is written into the code rather than read off the problem.

What the failure looks like from inside a solve #

Convergence factors are summary numbers, and it is worth saying what the residual history behind them looks like, because the shape is unusual.

The V-cycle at ε = 10⁻³ does not diverge, stall or oscillate. It converges, at 0.951 a cycle, which is to say it removes about 5% of the residual each time round. Ninety-two cycles would take it to 1e-2 of where it started. Every individual component behaves exactly as its own analysis says it should: the smoother damps what it can damp, the restriction restricts, the coarse solve is exact at the bottom, the prolongation prolongs.

There is no error message available. Nothing in the computation knows that the smooth error it is handing to the coarse grid is not the error that was left; the coarse grid solves the problem it is given, returns a correction that is right for the part of the error it can see, and the part it cannot see comes back untouched. Run enough cycles and the answer is correct — this is a convergent method, and the residual falls monotonically to 10⁻¹⁰.

What is lost is the property the entire field exists for. A method at 0.951 a cycle is a method whose cost is not proportional to the number of unknowns, and at that point one might as well be running conjugate gradients, which on this problem is simpler, needs no hierarchy and has a bound.

What a reader can check without running anything #

One consequence worth extracting, because it turns the essay into a diagnostic rather than a case study.

The quantity that predicts the trouble is the ratio of the largest to the smallest off-diagonal coupling in a row, and it can be read off the assembled matrix without solving anything, without computing an eigenvalue, and without knowing what the operator discretises. On the isotropic Laplacian that ratio is 1 in every interior row; at ε = 10⁻³ it is 1,000.

Nothing else in the standard toolkit reports it. The condition number is identical either way, as this essay’s first measurement establishes. The sparsity pattern is identical. The symmetry, the definiteness and the diagonal dominance are all unchanged. A matrix norm sees a scaling and not a direction.

So a practitioner whose multigrid has slowed down has a cheap first question — print the ratio of row couplings, at a few interior rows — and it distinguishes this failure from every other one in the field in a single pass over the matrix. That is also, not coincidentally, the exact quantity the algebraic method is built on: its strength-of-connection test is a threshold on that ratio, and the whole of what it does with the answer is decide which couplings to coarsen along.

Where the transition happens #

The failure is not a cliff, and the intermediate values are worth having because they say which quantity the trouble tracks.

At ω = 2/3 the smoothing factor runs 0.6667 at ε = 1, 0.9394 at ε = 0.1, 0.9934 at ε = 0.01 and 0.9993 at ε = 10⁻³ — which is 1 − 2ωε/3 to the accuracy of the arithmetic once ε is small, so the distance from one falls in proportion to ε. Each decade of anisotropy costs a decade of whatever margin the smoother had.

Allowing each ε its own best weight makes the picture no better, which is the point of computing the best weight at all. The optimum climbs from 0.8000 to 0.9565, 0.9951 and 0.9995 as ε falls, and the factor it buys climbs with it: 0.6000, 0.9130, 0.9901, 0.9990. The V-cycle’s convergence factor, run at each of those weights, comes out 0.2016, 0.7205, 0.9204 and 0.9565. Both quantities are monotone in ε, both saturate, and the second lags the first — which is the coarse-grid correction still doing its share of the work while the smoother stops doing its own.

The relationship between the two numbers is not a formula anybody has here. The smoothing factor bounds the two-grid factor under hypotheses the anisotropic problem does not satisfy, so reporting the pair rather than a ratio between them is the honest form, and the pair is more informative than either. What it establishes is that the smoothing factor is the leading indicator: it is at 0.91 where the V-cycle is still at 0.72, so a practitioner watching only the residual history sees a method that is merely slow, and would have seen, in a quantity that costs nothing to compute and requires no solve at all, a component already most of the way to useless.

What is asserted, and what refuses #

The condition number does not move. Four decades of ε, one number to twelve digits, checked against cot²(πh/2) by a route that computes no eigenvalues.

The measured damping is the symbol. One sweep on a pure mode multiplies it by a scalar, the scalar is a formula in the two frequencies, and they agree to 10⁻¹⁵ over five modes.

The scan is the closed form. Six values of ε and five weights, every one agreeing to 10⁻⁹, which is the site’s two-routes habit applied to a quantity that could easily have been left as a scan.

The point smoother fails, refused as a success. The assertion “point relaxation smooths the anisotropic problem” is fed to rejects() and must fail; it reports μ = 0.9993.

And coarsening the wrong direction is refused too. The assertion that semi-coarsening in x — the weak direction — repairs the smoother is fed to rejects(), and reports the same 0.9993. Half of the next two essays’ content is that the repairs have to be aimed, and an assertion that would have accepted either aim would not be checking that half.

What is left #

Two repairs, and they are the next two essays. Line relaxation solves a whole grid line at once, so the strong direction is inverted rather than swept over; semi-coarsening halves only the strong direction, so the modes the smoother cannot damp are handed to a coarse grid that still holds them. Both restore the smoothing factor to a third, from opposite ends, and the two of them turn out to have identical smoothing factors at every ε and every ω — which is a coincidence until the closed form is written down.

The rotated case, where the strong direction lies along neither grid axis. Neither repair applies there, and the phase’s honest stopping point is that nothing measured in it works — including the algebraic method built for exactly this reason.

Convection, which is the other way a second-order operator stops being self-adjoint and is not touched anywhere in this phase. An operator with a first-derivative term large enough to matter is not symmetric, so the whole apparatus of Fourier symbols on a symmetric stencil applies to it only approximately, and it deserves its own measurements rather than a paragraph in somebody else’s essay.