One step past the zero

A zero no twist can step around found where the circulant route to a banded system stops keeping its promise. A band that is nearly circulant can be solved by the circulant that wraps it — diagonalised by a Fourier transform — plus a correction of low rank for the corners the wrap gets wrong. On the Laplacian a twist of half a step keeps the wrap’s frequency samples off the symbol’s zero and the corrected route lands within a few times elimination. On the squared Laplacian, the pentadiagonal band $[1, -4, 6 + \varepsilon, -4, 1]$ whose symbol $(2 + s - 2\cos\theta)^2$ vanishes to fourth order, no twist can: the nearest sample is at best $\pi/n$ from the zero, where the symbol is $(\pi/n)^4$, and the wrap is as ill conditioned as the matrix. The corrected solve’s error followed κ(wrap)·u, at a tenth to three tenths of it, and at n = 64 it was 34 times behind elimination.

The essay proposed a remedy it did not try. The corrected solve’s loss is a solver’s loss, not the problem’s: the band itself is no worse conditioned than elimination finds it. A solver that is fast and not accurate enough is what iterative refinement exists for. Form the residual $r = b - Tz$ with the band matrix, which is cheap and exact up to rounding; solve $Td = r$ by the same fast route; add. If the fast route commits a relative error of about κ(wrap)·u, each step should multiply the error by about that, and since κ(wrap)·u is below one here the first step should take the error down to whatever floor the residual’s own rounding leaves. The measurement is whether it does.

The loss before refinement #

The starting point is the earlier picture. At a shift of $10^{-8}$ and a half-step twist, the squared Laplacian’s corrected solve is wrong by 2.6·10⁻¹³ at n = 16, 5.7·10⁻¹² at 32, 8.5·10⁻¹¹ at 64 and 1.4·10⁻⁹ at 128 — growing by a factor of fifteen to twenty a doubling, as $n^4$ does. Elimination with partial pivoting on the same band is wrong by 4.1·10⁻¹⁴, 3.2·10⁻¹³, 2.5·10⁻¹² and 4.2·10⁻¹¹. The corrected route is behind by 6, 18, 33 and 34 times. The band’s condition number is 3.3·10³, 4.3·10⁴, 6.1·10⁵ and 9.1·10⁶, and the wrap’s three to five times larger; elimination sits nine to twenty-six times below κ(T)·u, the familiar gap between the worst case a condition number bounds and the typical right-hand side, which the condition number is an amplifier measured directly.

One step #

At n = 64 one step takes the error from 8.5·10⁻¹¹ to 1.8·10⁻¹², below elimination’s 2.5·10⁻¹². The second step gives 1.9·10⁻¹² and the third 1.3·10⁻¹²: nothing further, the error wandering about the floor. Turn the dial and the same happens at every size. At n = 16 one step gives 2.4·10⁻¹⁴ against elimination’s 4.1·10⁻¹⁴; at 32, 2.2·10⁻¹³ against 3.2·10⁻¹³; at 128, 5.8·10⁻¹¹ against 4.2·10⁻¹¹. After one step the corrected route is 0.59, 0.68, 0.73 and 1.40 times elimination’s error — better than elimination at three of the four sizes, and within half again of it at the fourth.

Drawn against the size, the unrefined route rides κ(wrap)·u up at its $n^4$ slope, and the refined route lies on elimination’s line. The zero had put a factor of up to thirty-four between them, and one residual removes it. Nothing about the wrap was changed to do it: the same twist, the same four correction columns, the same four-by-four system.

Why one step is enough #

Refinement multiplies the error by the relative error of the correction’s solve, and that is where the zero’s cost went. The fast route commits an error of about a quarter of κ(wrap)·u relative to what it solves — at n = 64 about 7·10⁻¹¹ — so one step multiplies the error by that, taking 8.5·10⁻¹¹ to below 10⁻²⁰ in exact terms. It cannot get there, because the residual $b - Tz$ is formed in floating point: its rounding is about u times $\lVert T\rVert\,\lVert z\rVert$, and the correction solved from it carries that rounding multiplied by the band’s own conditioning. That is the floor — κ(T)·u times whatever fraction the right-hand side excites — and it is the same floor elimination sits at, because elimination’s backward error is also u times $\lVert T\rVert$ and its forward error is that times the same conditioning.

So after one step the corrected route and elimination are limited by the same thing: the problem, rounded once. The wrap’s larger condition number has dropped out entirely, because it multiplied an error that the first step made negligible. That is also why the second and third steps do nothing: they refine against a residual whose rounding is already the whole of the error, and each solves for a correction that is itself mostly rounding.

It is the result the correction lost to its own two-by-two solve found for a different loss — Cramer’s rule inside the capacitance system, removed by refinement against the band — and the reason is the same. A solver’s error, however it arises, is removed by refinement against the true matrix as long as the solver is accurate enough to contract, which here means κ(wrap)·u well below one. The fourth-order zero makes κ(wrap) large; it does not make it $1/u$ at any size here. At n = 128 κ(wrap) is 4.4·10⁷ and κ(wrap)·u about 5·10⁻⁹; the zero would have to push it past $10^{16}$ — a size of about 16,000 at this growth — before one step stopped being enough.

Backward stable after one step #

The forward errors say the two routes end in the same place. The backward errors say why, and they are the cleaner measurement.

Before refinement the corrected route’s backward error is 1.3·10⁻¹⁴ at n = 16 and 1.3·10⁻¹¹ at 128 — a hundred to a hundred thousand times the unit roundoff, growing with the wrap’s condition number. The route is not backward stable for the band: its answer is the exact solution of a problem measurably different from the one posed, and that is why its forward error exceeded elimination’s by more than the band’s conditioning could explain. After one step the backward error is 5.5·10⁻¹⁷ to 8.1·10⁻¹⁷ at every size, and elimination’s is 6.1·10⁻¹⁷ to 6.8·10⁻¹⁷. Both are below the unit roundoff, and they are the same number to within the scatter of five answers.

That is the most exact statement of what refinement did. It made the corrected route backward stable, at elimination’s own level, and a backward-stable solve of a given matrix has the forward error that matrix’s conditioning allows, whatever algorithm produced it. The exact answer to a nearby problem set out that division of blame between the algorithm and the problem; before refinement the corrected route took some of the blame for itself, and after one step it takes none.

Against the shift #

At n = 64 the shift barely matters past $10^{-4}$: the unrefined error is 4.6·10⁻¹¹, 7.6·10⁻¹¹ and 8.5·10⁻¹¹ at $10^{-4}$, $10^{-6}$ and $10^{-8}$, because the smallest symbol value the wrap sees is set by the grid’s $(\pi/n)^4$ once the shift is below it. At $10^{-2}$ the shift dominates, the symbol is bounded away from zero, and the unrefined error is 2.3·10⁻¹², only twice elimination’s 1.1·10⁻¹². After one step the refined route is 0.83, 2.9, 0.51 and 0.73 times elimination at the four shifts. The one at $10^{-4}$, three times elimination, is the exception in the whole measurement: the refined error there is 2.8·10⁻¹² against elimination’s unusually good 9.8·10⁻¹³, and a second step takes it to 1.6·10⁻¹², within two. The floor both are measured against is a median over five answers, and a factor of two at the floor is the size of its scatter.

What the step reads, and what it never needs #

It is worth being precise about what refinement asks of the band, because the answer is: almost nothing. The circulant route exists because a circulant is diagonalised by the Fourier transform — the matrix that is one row drew why — so every solve with it is two transforms and a division. The band itself is never factorised by this route, and refinement does not change that: it needs the band only to multiply a vector, five multiplications and four additions a row for the squared Laplacian, to form the residual. No pivot is chosen, no fill appears, nothing about the band’s own conditioning is exercised except through the rounding of that product.

That is the property that makes the combination work. The fast solve is inaccurate in a way that depends on the wrap’s conditioning; the residual is accurate in a way that depends only on the band’s entries and the arithmetic of a product. Refinement uses each for the thing it does well — the wrap to find a correction of roughly the right size, the band to say exactly how far off the current answer is — and the accuracy of the combination is the residual’s, not the solve’s.

A second consequence is that the refinement step does not care where the solve’s error came from. It removed the fourth-order zero’s amplification here; it removed Cramer’s rule’s cancellation in the earlier measurement; it would remove a wrap computed in lower precision, or a transform with its own rounding, by the same mechanism, as long as the solve contracts — an expectation the mechanism supports and none of which is measured here. The earlier essays spent their effort identifying which part of the corrected route was losing accuracy, and each identification was right. One residual makes the question moot for every one of them at once.

The same step, seen as an iteration #

The earlier essay observed that a wrap is used more often as a preconditioner than as a direct solver, and that an iteration measuring its residual against the band at every step never inherits the wrap’s error. Four orders of conditioning and four steps found a circulant-preconditioned iteration’s step count indifferent to four orders of the matrix’s conditioning, for the same reason: the preconditioner only has to be close, because the iteration corrects it.

One step of refinement is the simplest such iteration, stationary and unaccelerated, stopped after one correction. What the measurement adds is that for this operator one correction is all a direct solve needs — the second and third change nothing — so the direct route and the preconditioned iteration coincide after a single step, and the question of which to use becomes a question of how the solve is packaged rather than of how accurate it is.

What the route costs now #

The unrefined corrected solve applies the wrap’s inverse five times — once to the right-hand side and once to each of the four columns of the correction — and solves a four-by-four system. One refinement step adds a band multiplication, which is five multiplications a row, and one more application of the wrap to the residual; the four correction columns are the same and need not be recomputed. So the refined route costs about six applications of a Fourier transform pair where the unrefined one costs five.

That is the whole price of the thirty-fold accuracy, and it is small. What it does not change is the comparison with elimination on a one-dimensional band, which is linear in n and needs no transform at all: on the squared Laplacian in one dimension, elimination was always the cheaper route, and the circulant route’s case was never cost. Its case is the operator whose elimination is not linear — the two-dimensional analogue, whose band is as wide as the grid, where two dimensions and the cluster that thins measured the same wraps as preconditioners — and there the fourth-order zero arrives with the biharmonic operator, and the result here says its accuracy penalty is one refinement step, not a reason to abandon the route.

Where this leaves the zero #

A zero no twist can step around read its measurement as a limit on the route: every operator whose symbol has a multiple zero at the frequency of the smooth modes — the biharmonic operator, plate bending, higher-order smoothing penalties — would cost more than a constant through its wrap. Measured with one step of refinement, the limit is on the unrefined route only. The constant it costs becomes one extra wrap application and one band product.

The earlier essay also noted that a wrap is more often used as a preconditioner, where an iteration measures its residual against the band every step and so never inherits the wrap’s error. Refinement is the one-step version of that: a preconditioned iteration that stops after one correction because one is enough. Measured here, the direct route and the iterative one are the same route with different stopping tests, and the zero decided only how many corrections the direct one needed — one.

Two near-zeros cost less than one found that where the wrap’s samples fall relative to a symbol’s zeros decides the loss, and the twist was the lever for moving them. For a zero of order two it was enough; for order four it was not; with refinement after the solve it does not need to be. The twist still matters — it keeps the correction’s cancellation small, which keeps the first solve’s error small enough to contract — but the accuracy no longer rests on it.

What this does not settle #

One band, the squared Laplacian with a small shift, at four sizes up to 128 and one twist. The right-hand sides are generic, and a right-hand side concentrated on the symbol’s zero — the smooth mode the band barely damps — is the one on which the wrap’s conditioning would bite hardest; it is not singled out here. A zero of sixth or eighth order grows the wrap as $n^6$ or $n^8$ and would reach κ(wrap)·u of one at far smaller sizes, and there one step would no longer be enough and a second would be needed, or the iteration would not contract at all. That boundary is predicted by the argument above and not measured.

The cost accounting above is arithmetic, not measurement. The refinement step as measured recomputes the four correction columns rather than keeping them from the first solve, so it did six wrap applications’ work where a careful implementation would do one; the accuracy is identical either way, since the columns are the same numbers, but no timing here supports the claim that the refined route costs a sixth more than the unrefined one. The count is what the claim rests on, and it is the count a reader can check.

The residual is formed in working precision. In extended precision refinement would push past elimination to the forward error the rounded data allow, and the comparison with elimination would change sides; that is the classical mixed-precision refinement and it is not the question here.

Still open: the order at which one step stops sufficing, and the factorised route #

Zeros of higher order. The argument says one step suffices while κ(wrap)·u is well below one, and κ(wrap) grows as $n^p$ for a zero of order p. Cubing the Laplacian, or raising it to the fourth power, reaches that boundary at sizes that can be measured, and finding the size at which the refined route first needs a second step — and checking it against $n^p u \approx 1$ — would turn the argument into a rule for how many steps a given operator needs.

Two Laplacians instead of their square. The squared Laplacian is almost the product of two Laplacians, each of which the twisted route solves within a few times elimination; the band and the product differ only in their first and last rows. Solving the product by two twisted solves and correcting the boundary rows would compound two small errors instead of inheriting a fourth-order zero, and whether that is cheaper than one refinement step — it trades a second wrap application for a rank-two correction — is a comparison of the two ways to spend the same budget.