Smoothing a whole line at once

A direction the smoother cannot see leaves the method at 0.9565 a cycle and locates the failure precisely: the mode that oscillates across the weak direction and is smooth along the strong one is damped by 1 − 2ωε/(1 + ε), which is 0.99867 at ε = 10⁻³, and the coarse grid cannot take it because it is not smooth.

There are two ways out and this essay is the first. It changes what the smoother inverts.

A point relaxation asks, for each unknown in turn, what value would satisfy that unknown’s own equation given the current values of its neighbours. A line relaxation asks the same question of a whole grid line at once: what values along this column would satisfy every equation on it exactly, given the current values of the columns either side. That is a tridiagonal system, it costs a constant number of operations per unknown to solve, and it inverts the direction the anisotropy made strong instead of iterating along it.

The smoothing factor goes from 0.9993 to 0.3340, which is the one-dimensional ⅓ arriving on a problem that is not one-dimensional.

What a line smoother is, written as a splitting #

Every stationary method on this site is a splitting: choose a matrix M that is easy to invert, write A = M − (M − A), and iterate x ← x + ωM⁻¹(b − Ax). Jacobi takes M to be the diagonal; Gauss–Seidel takes it to be the lower triangle; the stationary essay measures the rates both of them achieve on the one-dimensional model problem against closed forms.

A y-line smoother takes M to be the entries of A that couple points within the same grid column. For the five-point operator that is the diagonal and the two vertical off-diagonals, so M is block diagonal with k tridiagonal blocks of size k, and inverting it is k tridiagonal solves. The horizontal couplings are outside M and are therefore lagged — treated as part of the right-hand side, exactly as Jacobi lags everything off the diagonal.

Written that way the cost is easy to state. A point sweep is one multiply-add per stored entry plus a division; a line sweep is the same defect computation plus a Thomas solve, which is about five operations a point. So a line sweep costs roughly three times a point sweep, and buys a smoothing factor of 0.334 where the point sweep gets 0.999. Three times the work for a component that works is not a difficult trade.

The implementation is written against the matrix rather than against the stencil — for each line it reads the within-line couplings straight out of the rows of A — and that choice does two things. It makes the same code smooth the nine-point rotated operator, which the last essay in this sequence needs. And it makes the failure on that operator a consequence of the definition rather than a separate claim: everything the rotation adds lies outside the line, so a line smoother lags it.

The symbol, and where the ε goes #

The analysis is the same three lines as before with one substitution. The symbol of A at the frequency pair (θx, θy) is 4ε sin²(θx/2) + 4 sin²(θy/2). The symbol of M, for a y-line smoother, is what is left when the horizontal couplings are dropped:

2ε + 4 sin²(θy/2)

and the damping factor is 1 − ω·A/M. Take the mode that broke the point smoother — θx = π, θy → 0. There A → 4ε and M → 2ε, so the ratio is 2 and the factor is 1 − 2ω, which at ω = 2/3 is −⅓.

The ε cancelled. That is the whole repair in one line: the weak coupling still sets the size of what the smoother has to remove, and it now also sets the size of what the smoother inverts, so the two divide out and the mode is damped like any other.

The maximum over the frequencies the coarse grid cannot see comes out

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

against the point smoother’s max(|1 − ωε/(1+ε)|, |1 − 2ω|). One coefficient differs, and it is the ε. A scan over 128×128 frequency pairs agrees with the expression to 10⁻⁹ at six values of ε and five weights.

Two numbers fall out that the one-dimensional field never had reason to write down. The optimum is where the two branches cross:

ω* = 2(1 + ε)/(3 + 2ε),    μ* = (1 + 2ε)/(3 + 2ε)

At ε = 1 that is ω = 4/5 and μ = 3/5, which is the two-dimensional isotropic answer from the previous essay. As ε → 0 it is exactly 2/3 and ⅓ — the one-dimensional numbers, recovered on a two-dimensional problem, because in the limit the operator decouples into k independent one-dimensional problems along the columns and the line smoother solves each of them exactly.

The direction has to be right, and nothing checks it #

Relaxing along x — the weak direction — repairs nothing. The smoothing factor is 0.9993, the same number the point smoother gets, because a line smoother in the weak direction inverts the couplings that were never the problem and lags the ones that were.

This is worth dwelling on because it is where the repair’s limitation lives. There is nothing in the code that reads the operator and decides which direction to relax along. Somebody wrote smoother: "liney", and if the anisotropy had run the other way that would have been the wrong choice, at no cost in warnings: the method converges either way, at 0.037 or at 0.967, and the residual history is the only thing that says which.

On a problem whose coefficients vary in space — which is every problem anybody builds this machinery for — no single choice is right everywhere. Alternating line relaxation, which sweeps the columns and then the rows, is the standard answer and costs twice as much; it is not measured here and it is a real omission rather than a hidden one.

What it achieves as a method #

The smoothing factor is a statement about one component. The V-cycle built on it, at ε = 10⁻³ on a 31×31 grid, converges at 0.0370 a cycle — better than the isotropic method’s 0.2016, which is worth explaining rather than celebrating.

A line solve is a stronger smoother than a point sweep in every direction, not only the strong one: it removes the whole of the error’s variation along each column exactly, so what it hands to the coarse grid is smoother than what a point smoother hands over on the isotropic problem. The cost is the factor of about three per sweep, so per unit of work the two are comparable and the anisotropic solve is not actually cheaper than the isotropic one. Reporting the cycle count without the cost per cycle would say it was.

The rate is flat in the size, which is the property the whole field is about: 0.0370, 0.0370 and 0.0371 at k = 15, 31 and 63. And it is worth noting what has not been repaired — nothing about the operator changed, the condition number is what it always was, and a Krylov method on this problem is exactly as fast or slow as it was before.

How much of that is the coarse grid, and the answer is: less and less #

A convergence factor of 0.037 is better than the isotropic method achieves, and a number that good invites the question of what is producing it. The way to ask is the one the one-dimensional field used: run the same cycle with the coarse-grid correction removed, leaving three line sweeps and nothing else, and see how much is lost.

At ε = 1 the answer is decisive. The V-cycle contracts at 0.271 and the smoother alone at 0.876, so the coarse grid is doing nearly all of the work, exactly as the field’s opening essays say.

As the anisotropy grows that stops being true, and the crossover is gradual and complete:

By ε = 10⁻³ the line smoother on its own converges at 0.083 a cycle, and the whole hierarchy below it improves that to 0.038. The coarse grid has stopped being the method and become a refinement.

The reason is not subtle and it is the same limit as before. As ε → 0 the operator decouples into k independent one-dimensional problems along the columns, a y-line solve solves each of them exactly, and a smoother that solves the problem is not a smoother — it is a direct method. So the impressive 0.037 is partly the anisotropy having made the problem easy for this particular smoother, and quoting it as multigrid’s achievement would be crediting the hierarchy for something the tridiagonal solve did.

The honest reading is the middle of that table. At ε between 0.03 and 0.3, both components are contributing and the V-cycle is several times better than either half, which is what a working multigrid method looks like.

The trade, in work rather than in cycles #

A convergence factor compares methods only if their cycles cost the same, and these do not. It is worth doing the arithmetic once.

A point sweep on the five-point operator touches five entries a row, so it is five multiply-adds and a division: about seven operations an unknown. A line sweep computes the same defect — five multiply-adds — and then runs a Thomas solve, which for a tridiagonal system is two operations per unknown going forward and two coming back, plus the divisions: about thirteen operations an unknown in total, so roughly twice a point sweep rather than the three times a rough count suggests. The tridiagonal factorisation itself is computed once per line per sweep here; a production code would factor the lines once and reuse them, which takes the ratio nearer 1.5.

Set against that, at ε = 10⁻³ the point method needs 51 cycles to reduce the residual by a decade and the line method needs 0.7. The trade is not close, and it is not close by a factor large enough that the exact cost ratio does not matter — which is the useful form of the statement, since the ratio depends on the implementation and the factor of seventy does not.

Where the ratio would matter is the isotropic problem, and there the answer goes the other way: a line smoother’s V-cycle contracts at 0.271 against a point smoother’s 0.202, so it is both slower per cycle and more expensive per cycle. A line smoother is not a better smoother; it is a smoother aimed at a particular difficulty, and on a problem without that difficulty it is a worse trade in both currencies at once.

The blocks are solvable, which is not free #

The Thomas algorithm used for each line solve does no pivoting, and a tridiagonal solve without pivoting is only safe on a matrix that does not need it. So that is asserted rather than assumed.

Every y-line block has diagonal 2ε + 2 and within-line off-diagonal magnitudes summing to at most 2, so the ratio is (2ε + 2)/2 = 1 + ε > 1 and every block is strictly diagonally dominant. The measurement runs at three parameter settings — isotropic, anisotropic and rotated — and reports the worst ratio found: 2 at ε = 1, 1.001 at ε = 10⁻³, and 2 for the rotated operator at 45°.

The 1.001 is the interesting one. The blocks are dominant by a margin that vanishes with ε, which is exactly the regime where a Thomas solve would be expected to lose accuracy, and it does not: the factorisation of a tridiagonal M-matrix has no growth, and the residual of the line solves is at rounding throughout. A margin of 1.001 is a margin.

What is asserted, and what refuses #

Every smoother’s measured damping is its symbol. A grid mode is an eigenvector of the line smoother’s iteration matrix as well as the point smoother’s — a fact that stops being true the moment the operator is rotated, and is therefore checked rather than quoted. Five modes, three smoothers, agreement to 10⁻¹⁵.

The scan is the closed form, at six values of ε and five weights, to 10⁻⁹.

Line relaxation in the weak direction does nothing, asserted as a positive statement rather than as an absence: μ > 0.9 for x-line relaxation at ε = 10⁻³, reported alongside the 0.334 for y.

The line blocks are strictly diagonally dominant, at three parameter settings, with the worst margin printed.

And the claim that a point smoother copes is refused. assertThat(μ < 0.5) for point relaxation at ε = 10⁻³ is fed to rejects() and reports 0.9993.

What is left #

Semi-coarsening, which is the other repair and is the next essay. It leaves the smoother alone and changes the coarse grid instead, and its smoothing factor is identical to this one’s at every ε and every ω — which the closed forms explain and which the convergence factors then contradict, by a factor of about three.

Alternating line relaxation, sweeping columns and then rows, which removes the need to know which direction is strong at the cost of doubling the smoother. It is the standard production answer to a spatially varying anisotropy and it is not measured here.

The rotated case. A line smoother cuts across a strong direction at 45°, which is a stronger statement than “it is not aligned”: the couplings the rotation adds are outside every line, so a line smoother lags all of them by construction. That is measured in the essay this sequence ends with, at 0.784 a cycle, and it is the measurement that makes the algebraic method necessary rather than merely elegant.