A zero no twist can step around

The circulant the problem did not contain solved a banded matrix through the circulant that differs from it in its corners — a fast transform for the circulant, a small capacitance system for the corners — and found it losing accuracy against elimination. The correction lost to its own two-by-two solve found most of that loss was the small system’s, solved by Cramer’s rule, and that with a pivoted solve what remained was the cancellation: the intermediate solve through the wrap is larger than the answer, the correction’s last line subtracts two large vectors, and that ratio times u is the error. It drew the lesson plainly. “What is left is the answer, with the rounding of the cancellation and no more.”

Its last paragraph named the case that could break this. “A banded matrix whose symbol has a zero of high order — the square of the Laplacian, whose symbol vanishes to fourth order at the origin — brings several nearly null samples in at once from a single zero, and a correction of fixed rank then has fewer columns than near-null directions. Whether a stable small solve is still enough there, or whether the cancellation itself grows faster than the wrap’s condition number, is the measurement the rank-four band cannot make.”

Neither alternative is what happens. The cancellation stops growing. The error does not, and what it follows is the wrap.

Two symbols, and what a twist can reach #

The operators are the Laplacian, tridiag(−1, 2 + s, −1), whose symbol $2 + s - 2\cos\theta$ vanishes to second order at θ = 0 as s → 0, and its square, the pentadiagonal band with symbol $(2 + s - 2\cos\theta)^2$, which vanishes to fourth order there. Each is the circulant of its symbol with its corners removed, so each can be solved through the circulant and a correction — of rank two for the Laplacian, rank four for its square — and each circulant can be twisted by a phase on its corners, which moves the frequencies it samples from $2\pi k/n$ to $(2\pi k + \varphi)/n$.

The earlier essays twisted by half a step, φ = π, for the Laplacian, and it was the right repair: the untwisted wrap samples θ = 0 itself, where the symbol is s, and the twisted one samples π/n and nothing nearer. That is as far from the zero as a twist can move the nearest sample, since the samples are 2π/n apart. At π/n the Laplacian’s symbol is about $(\pi/n)^2$, 2.4·10⁻³ at n = 64, and the wrap’s condition number — the largest sample, 4, over the smallest — is 1,659, the same as the banded matrix’s own 1,710. The twisted wrap is as well conditioned as the problem.

For the square, the nearest sample sees $(\pi/n)^4$, 5.8·10⁻⁶ at n = 64, against a largest value of 16. The twisted wrap’s condition number is 2.8·10⁶, and the banded matrix’s own is 6.1·10⁵: the wrap is 4.6 times worse conditioned than the matrix it stands in for, and at every size measured it is between 3.2 and 4.8 times worse. No twist does better, because no twist puts a sample further than π/n from the zero, and the fourth power is what the symbol does there.

The corrected solve, twisted every way #

At n = 64 and s = 10⁻⁸, with every capacitance system solved by pivoted elimination and the median taken over five right-hand sides, the Laplacian’s corrected solve at a half-step twist has an error of 4.6·10⁻¹⁴ against elimination’s 1.1·10⁻¹⁴: four times. The squared Laplacian’s is 8.5·10⁻¹¹ against elimination’s 2.5·10⁻¹²: thirty-four times. Over the seven twists drawn, from a twentieth of a half-step to a whole one, the squared Laplacian’s best is 5.6·10⁻¹¹ at nine tenths of a half-step, twenty-two times elimination, and its worst, at the smallest twist, is 4.6·10⁻⁶.

Turn the dial. At n = 16 the squared Laplacian’s corrected solve is 6.3 times elimination’s; at 32, 18 times; at 64 and 128, 34 and 33 times. The Laplacian’s stays between 1.1 and 6 times at every order. The twist that repaired the second-order zero leaves the fourth-order one an order and a half of magnitude behind.

What the error follows #

That figure is where the cancellation law was established, and its regime is visible in it: the periodic wrap samples the zero, the cancellation grows as the shift falls, and the pivoted curve rides on the cancellation times u because the cancellation is by far the largest thing in the problem.

The preceding essay’s law would put the error at the cancellation times u. On the squared Laplacian at n = 64 the cancellation is 1.0·10⁴, so the law predicts 1.1·10⁻¹², and the error is 74 times that; at the other sizes it is 25 to 115 times. The open dots on the figure are that prediction, and they sit below every measured error.

What the error does follow is the wrap’s condition number. Across four sizes and seven twists, the squared Laplacian’s error is between 0.08 and 0.30 times κ(wrap)·u — every case, over seven orders of magnitude of κ(wrap). The Laplacian at a half-step twist is on the same line, at 0.11 to 0.25. And that explains the Laplacian’s own small excess over elimination too: its cancellation at n = 64 is 7.6, and its error is 55 times the cancellation times u, because the floor is the wrap’s condition number and not the cancellation. The cancellation law was measured in the regime where the cancellation was enormous — the periodic wrap nearly landing on the zero — and there it dominates. Keep it small, and the floor under it is the wrap’s.

That is a statement the earlier essays could not have made, because on the Laplacian the two floors nearly coincide. The twisted wrap’s condition number is the matrix’s, and a quarter of κ(T)·u is within a few times of what elimination achieves. On the fourth-order zero they separate: the wrap is worse conditioned than the matrix by a factor near five, and the corrected route inherits every bit of it.

The band’s conditioning, and the wrap’s #

The band matrices themselves grow ill conditioned with the order of their zero, and a limit the matrix never reaches measured how a banded Toeplitz family’s condition number approaches what its symbol predicts. Here the Laplacian’s grows as $n^2$: 116, 441, 1,710 and 6,740 at orders 16 to 128. Its square’s grows as $n^4$, near enough: 3.3·10³, 4.2·10⁴, 6.1·10⁵ and 9.1·10⁶. Both are what the symbol says, since the smallest eigenvalue of the band sits at the smallest frequency it can represent and the symbol there is a power of that frequency.

The wrap has the same growth and a different constant. Its smallest sample is at π/n, a little nearer the zero than the band’s smallest frequency, and for a second-order zero the difference is a factor between 0.89 and 0.98 in the condition number — nothing. For a fourth-order zero the same small difference in frequency is raised to the fourth power, and the wrap’s condition number is 3.2, 4.0, 4.6 and 4.8 times the band’s at the four orders. That factor, and the constant the corrected route pays against the constant elimination pays, are the whole of the thirty-fold gap; neither is large, and their product is.

The capacitance matrix, measured #

The small system is ill conditioned too, and for the reason the rest of this essay gives rather than the one the earlier essays analysed. At a half-step twist the Laplacian’s 2 × 2 capacitance matrix has condition numbers of 17, 33, 65 and 129 at the four orders — growing as n, since each is the wrap’s smallest sample against the correction’s scale. The squared Laplacian’s 4 × 4 capacitance matrix has 3.1·10³, 2.6·10⁴, 2.1·10⁵ and 1.7·10⁶, growing as $n^3$. The error is not a fixed multiple of that either: it is 0.8, 2.0, 3.6 and 7.3 times κ(S)·u at the four orders, so the capacitance matrix’s condition number under-predicts the loss by a growing factor, where the wrap’s predicts it within a constant.

That is consistent with where the loss comes from. The capacitance matrix is formed from the wrap’s inverse applied to the correction’s columns, so it carries the wrap’s ill-conditioning in a four-dimensional projection; the intermediate solve carries it in all n dimensions. The small system’s pivoted solve is backward stable, and whatever it loses is bounded by its own condition number; the intermediate solve’s rounding is amplified by the whole wrap, and nothing in the correction’s four columns undoes that.

Why a stable small solve cannot absorb it #

The preceding essay explained why an ill-conditioned capacitance matrix cost nothing once it was solved stably: as the wrap nears a zero, its inverse applied to anything has one enormous component along one near-null vector, the capacitance matrix and its right-hand side inherit that component together as a nearly rank-one term, and a backward-stable solve of the small system returns the coefficients that cancel it exactly. The ill-conditioning is concentrated, and what is concentrated can be removed.

On the squared Laplacian at a half-step twist it is not concentrated. The two samples at ±π/n see 5.8·10⁻⁶; the next pair, at ±3π/n, see 81 times that; the pair after, 625 times. The wrap’s inverse amplifies a spread of directions by factors from 10⁵ downwards, graded rather than separated, and the rank-four correction’s four columns cannot pick out a single large term to cancel because there is none. The rounding committed in the fast transform, which is backward stable in the norm and spread across every frequency, is amplified by the wrap’s inverse along all of those directions at once, and the capacitance matrix, whose condition number here is 2.1·10⁵, is ill conditioned for the same reason rather than for the one the earlier essay analysed.

So the rank-four correction does not have “fewer columns than near-null directions” in any sharp sense: at a half-step twist the nearest two samples form one conjugate pair, and a rank-four correction has room for two such pairs. What it lacks is a gap. Two near-zeros cost less than one found that a pair of samples landing on a pair of zeros gave a capacitance matrix of condition number one, because the correction absorbed exactly the directions the wrap had lost; here the wrap has lost a whole graded family of directions and the correction absorbs none of them cleanly.

Against the shift #

Through the periodic wrap, which samples the zero itself, the squared Laplacian’s corrected solve loses two digits for every decade of the shift: 1.5·10⁻¹² at s = 10⁻², 1.9·10⁻¹⁰ at 10⁻³, 2.1·10⁻⁸ at 10⁻⁴, 1.4·10⁻⁶ at 10⁻⁵, 1.0·10⁻⁴ at 10⁻⁶. That is the cancellation law at work, with the wrap’s smallest sample at $s^2$ and the intermediate solve growing as its inverse, and it is the Laplacian’s first essay again with the exponent doubled.

The half-step twist stops that growth: 2.3·10⁻¹², 1.5·10⁻¹¹, 4.6·10⁻¹¹, 1.0·10⁻¹⁰ and 7.6·10⁻¹¹ over the same shifts, a floor near 10⁻¹⁰ once s is below the nearest sample’s (π/n)⁴. Elimination sits at 1 to 4·10⁻¹² throughout. The twist removes the loss that depends on the shift and leaves the loss that depends on the order, and for this symbol the second is the larger.

What this means for the route #

The fast route’s promise was that a banded matrix could be solved in O(n log n) with a small correction and lose nothing against elimination but a constant. For symbols with a simple zero — the Laplacian, and everything the earlier essays measured — the half-step twist keeps that promise within a few times. For a zero of order p the twisted wrap’s condition number grows as $n^p$ while the band’s own grows as well, but the wrap’s is larger by a constant that grows with p, and the corrected route’s error follows the wrap. The biharmonic operator, plate bending, higher-order smoothing penalties: every operator whose symbol has a multiple zero at the frequency of the smooth modes is in this class, and on them the twisted route costs more than a constant.

A code that wants the route on such an operator has two choices, neither measured here. It can refine: one step of iterative refinement with the band matrix’s own residual, which costs a band multiplication, should remove most of a thirty-fold loss, as the correction lost to its own two-by-two solve found refinement doing for Cramer’s loss. Or it can factor the operator: the squared Laplacian is the square of a matrix the twisted route solves to within four times elimination, and two such solves in sequence would multiply those factors rather than inherit $n^4$.

The same wrap as a preconditioner #

A wrap is used far more often as a preconditioner than as a direct solver, and there the loss measured here matters less. The circulant that cannot be indefinite and two dimensions and the cluster that thins measured circulant preconditioners for banded Toeplitz systems by the steps an iteration takes, and four orders of conditioning and four steps found the count indifferent to four orders of magnitude of the matrix’s conditioning. Preconditioned by its own twisted wrap, a band is the identity plus a correction of rank four, so the preconditioned matrix has at most five distinct eigenvalues and a Krylov method would converge in five steps in exact arithmetic whatever the wrap’s condition number, with each step’s residual computed from the band itself.

That residual is what the direct route lacks. The corrected solve is a single application of the wrap’s inverse and a correction, and its error is whatever that application commits. An iteration measures its residual against the band after every step and removes what the wrap got wrong, so the κ(wrap)·u floor becomes a matter of the iteration’s stopping test rather than a limit on accuracy. The fourth-order zero makes the twisted wrap a poor direct solver and leaves it a good preconditioner, which is the use it is most often put to.

What this does not settle #

Two operators, shifts from 10⁻² to 10⁻⁸, four sizes up to 128, five right-hand sides with random entries. The fraction of κ(wrap)·u the corrected route pays, a tenth to three tenths, is measured and not explained; nor is why elimination pays a twentieth to a tenth of κ(T)·u on both bands, which is what makes the wrap’s excess over the matrix visible at all.

The explanation of the floor — normwise-stable transform rounding amplified along a graded family of near-null directions — is an argument from the measurements, not a measurement of the error’s direction. Projecting the error onto the wrap’s eigenvectors would test it.

The Laplacian at small twists falls well below the law: at a twentieth of a half-step its error is 0.008 of κ(wrap)·u. There its one near-null sample is the direction the rank-two correction removes, and the wrap’s condition number, set by that sample, is not what the error sees. That exception is the preceding essay’s mechanism, and the law drawn here is for the regime where it does not apply.

Still open: refinement, and the factorised route #

One step of refinement. A residual computed with the band matrix and a second corrected solve against it should remove most of the thirty-fold loss, if the corrected route’s error is a backward error of size κ(wrap)·u in the norm. Whether one step brings the squared Laplacian to elimination’s accuracy, and what it costs against simply eliminating, is the measurement that decides whether the route survives on high-order zeros.

Solving the square as two Laplacians. The squared Laplacian with the right boundary rows is the product of two tridiagonal matrices the half-step twist handles within a few times elimination. Two twisted solves in sequence would compound two small factors rather than inherit a wrap of order $n^4$, and whether the product form’s boundary rows match the band’s closely enough for that to be the same problem is the question the band’s corners decide.