The circulant the problem did not contain

The matrix that is one row ended on what breaks a circulant. Change one entry, or stop the diagonals wrapping round, and the Fourier vectors no longer diagonalise the matrix, the eigenvalues are no longer a transform of anything, and the closed form is gone. It said the spectrum became “recoverable only through the correction formulae or by factorising”, and left the correction formulae there.

They are the obvious next move and they are standard. A matrix that differs from a circulant C in k entries is C plus a matrix of rank at most k, and the Sherman–Morrison–Woodbury identity solves it with k + 1 solves against C and one dense k × k system. Each solve against C is two transforms and a division. So a matrix with $n^2$ entries and a structure that is broken in two places costs a handful of transforms and a two-by-two solve — n log n, not $n^3$ — and the answer is exact algebra.

The matrix this essay uses is the simplest one with that shape. T = tridiag(−1, d, −1) is the second-difference matrix with its ends held fixed, and at d = 2 + σ it is the discrete operator −u″ + σu on an interval. Its periodic neighbour C has the same first row except that the −1 which falls off the end of the row wraps round to the far corner. So T = C + (two corner entries), a rank-two correction, and the transform route applies.

Why anyone routes a tridiagonal solve through a transform #

It is worth being plain about the model before measuring it, because a reader will object at once and be right. A tridiagonal matrix does not need a transform. Elimination on it costs about 8n operations, the transform route costs three solves at two transforms each — some 30n log₂n — and at every size drawn here elimination is the cheaper of the two by a factor of twenty or more.

The tridiagonal case is on the page because it is the smallest matrix with the shape that matters and because every quantity in it has a closed form, so each measurement can be checked against a formula rather than against another measurement. The shape it stands in for is common and expensive. A dense matrix whose entries depend only on i − j except near its boundary — a convolution with a kernel that does not wrap, an integral operator on an interval, a periodic discretisation with one interface — is a circulant plus a correction confined to a corner block, and there the comparison is $n^3$/3 against a few transforms and a k × k solve. The first essay’s table put that ratio at 546 by n = 256. The question is whether the accuracy survives the trip, and the tridiagonal matrix answers it with nothing else in the way: every digit the correction loses below, it loses to the route and not to the problem.

Two matrices a corner apart, with nothing in common #

The claim worth testing is not the cost, which is arithmetic. It is the one the cost claim silently carries: that routing a solve through C instead of through T changes the price and not the answer.

A correction cheaper than the problem already measured the general form of the objection. Every quantity in the Woodbury identity is a solve against the base matrix, so the errors are the base matrix’s, and on a construction where the updated matrix was the identity the rank-one formula lost digits in proportion to κ of the matrix that was replaced, at a fitted slope of 1.00. That essay chose its base matrix to be bad. The question here is what happens when nobody chose anything — when the base matrix is simply the circulant that sits next to the matrix in hand, picked because it is the one with a fast solve.

For T and C the two condition numbers can be written down. T’s eigenvalues are d − 2cos(πj/(n + 1)) for j = 1, …, n, the smallest is σ + 4sin²(π/(2(n + 1))), and at n = 64 that stops κ(T) growing at 1,712 however small σ becomes. C’s eigenvalues are the same function sampled at θ = 2πk/n, and one of those samples is always θ = 0, where the function is exactly σ. So κ© = (4 + σ)/σ: 4·10⁴ at σ = 10⁻⁴, 4·10⁸ at 10⁻⁸, infinite at zero.

At σ = 0, the matrix people actually meet — the plain second difference with fixed ends — has condition number 1,712 and its periodic neighbour is singular, at every size, for a reason visible in the first row: the row sums to zero, so the constant vector is in the null space. That is the circulant field’s own remark that “a circulant is singular exactly when its symbol vanishes at a root of unity”, met at the one root of unity that every size contains. It is also the vertex nobody solves for: the graph Laplacian of a ring is this C, and a Laplacian with nothing grounded has the constants in its kernel.

The picture makes the asymmetry concrete at n = 16. The large dots are C’s sixteen eigenvalues and one of them sits on the curve’s minimum, at |λ| = 10⁻⁴. The small dots are T’s, mirrored to fill the turn, and the nearest of them to zero is 0.0342 — the first sample is a sixteenth of a turn in. Same curve, same entries in every interior row, and condition numbers of 4·10⁴ against 116.

The correction follows the wrap, and faster than it #

Now solve. Take a right answer x from a seeded Gaussian draw, form b = Tx, and recover x twice. Five draws at each σ, and the median is what is plotted.

Elimination on T returns 1.1·10⁻¹⁴ at σ = 10⁻⁸ and never gets worse than 1.8·10⁻¹⁴ across the whole range — κ(T)u is 1.9·10⁻¹³, and elimination does better than its bound, as it usually does. The correction through C returns:

The first column that matters is the fourth, read against the third. At σ = 10⁻⁴ the correction does better than κ©·u. By σ = 10⁻⁶ it is 25 times worse than κ©·u, by 10⁻⁸ it is 2,700 times worse, and by 10⁻¹⁰ the answer has no correct digit in it at all, on a matrix whose condition number is 1,712. The slope of the correction’s error is very nearly two decades per decade of σ where κ©'s is one.

So the rank-one result does not transfer at its own slope. The loss here is not κ of the replaced matrix; it grows like the square of it.

Where the digits go #

The identity is short enough to follow the arithmetic through it. With U the two unit vectors at the ends and V the two corner couplings, T = C + UVᵀ, and

Two things in that line can amplify rounding, and both can be measured separately.

The first is S, the two-by-two capacitance matrix. It is built from the corner entries of the two columns $C^-$¹e₀ and $C^-$¹eₙ₋₁, and near singularity both columns are dominated by the same thing — the constant vector divided by σ — the condition number is an amplifier, and this is the direction it amplifies. So S is a two-by-two matrix whose four entries are each about 1/(nσ) plus something of order one, and it is nearly rank one. Its condition number is measured at 6.37·10², 6.18·10³, 6.16·10⁴, 6.15·10⁵ and 6.15·10⁶ for σ from 10⁻⁴ to 10⁻⁸ at n = 64. Multiply each by nσ and the answer is 4.08, 3.96, 3.94, 3.94, 3.94: κ(S) = 4/(nσ), which is κ©/n, to three figures.

The second is the last subtraction. The intermediate y = $C^-$¹b is enormous for the same reason — b has a component along the constant vector, C divides it by σ — and the answer x is of order one. So x is formed as the difference of two vectors each about ‖y‖ in size, and whatever relative error those vectors carry is multiplied by ‖y‖/‖x‖. That ratio reads 1.6·10², 1.6·10³, 1.6·10⁴, 1.6·10⁵, 1.6·10⁶ over the same range. Cancellation takes the answer, not a digit is the general statement: the subtraction itself is exact, and it exposes error the operands were already carrying.

Put the two together and the forward error divided by the unit roundoff should be of the order of their product. It is — the ratio of the one to the other is 0.22, 0.11, 0.07, 0.11 and 0.07 between σ = 10⁻⁵ and 10⁻⁹, a constant of about a tenth that does not trend. Both factors are proportional to 1/σ, so the error is proportional to 1/σ², and that is the slope of two.

The mechanism is worth saying in words, because it is not specific to circulants. A rank-k correction solves a k × k system whose entries are values of the base matrix’s inverse. If the base matrix has a near-null direction that the correction’s columns can see, every entry of that system carries the same huge component and the system is nearly singular in exactly the proportion the base matrix is. At rank one the “system” is a scalar, and a scalar has condition number one whatever its size — which is why the rank-one measurement found a slope of one and not two. At rank two the capacitance matrix has a condition number of its own, and it is the base matrix’s divided by n.

The size of the matrix enters twice, in the right direction #

Both amplifications fall with n. κ(S) = 4/(nσ) says so directly, and the cancellation falls too, because the corner entries through which b reaches the constant vector are a smaller share of a longer vector.

At n = 16 the correction returns 2.3·10⁻⁷ at σ = 10⁻⁶, 3.2·10⁻³ at 10⁻⁸ and 40 at 10⁻¹⁰. At n = 256 the same three read 5.9·10⁻¹⁰, 3.2·10⁻⁶ and 0.044. Sixteen times the size buys between four hundred and a thousand times the accuracy at every σ.

That dependence is the opposite of the usual one. T gets harder as it grows — κ(T) goes 116, 1,712, 26,768 across the three sizes, the familiar $n^2$ of a second difference — and elimination’s error rises from 10⁻¹⁵ to 10⁻¹³ with it. The correction gets easier, because what it is exposed to is not T’s conditioning but the size of C’s near-null component as seen through two corners, and that shrinks. At a large enough n the two lines would cross. They do not cross in the range drawn, and at σ = 10⁻⁸ the gap between them is still 10¹² at n = 16, 10¹⁰ at n = 64 and 3·10⁷ at n = 256 — closing by a factor of a hundred or more for each factor of four in size.

Refinement repairs it, until it cannot start #

The rank-one measurement had a repair and it is worth trying here because the residual is cheap. Compute r = b − Tx̂ with T itself, which is a three-term stencil and involves no inverse of anything, then solve for the correction by the same route, and add it. The residual is exact to rounding; only the correction step is exposed. It is the same division of labour a small residual is not a small error insists on: the residual is the trustworthy number and the answer is the thing it is used to repair.

At σ = 10⁻⁸ the error goes 1.2·10⁻⁴, 1.2·10⁻⁸, 4.2·10⁻¹², 6.3·10⁻¹⁵ — rounding in three passes. At 10⁻⁹ it starts at 7.6·10⁻³ and needs six. At 10⁻¹⁰, where the unrefined answer is already wrong by 1.3, it diverges: 1.3, 10, 39, 230, 1,900, 12,000, 38,000.

The contraction factor per pass is the error of the solver, because the correction to the correction is computed by the same inaccurate route — which is the rank-one essay’s rule, that a pass reduces the error by the solver’s own relative accuracy. What changes is where that rule stops working. There, the solver’s accuracy was κ(A)·u, so refinement converged for any base matrix with κ below 1/u, about 10¹⁶. Here it is roughly κ©²u/n² times a tenth, and it reaches one at κ© near 10¹⁰ at n = 64. The route stops being repairable six orders of magnitude earlier than the rank-one rule says, and the reason is the capacitance matrix the rank-one case does not have.

The gap refinement can close drew the line between a bad route and a bad question, and this is squarely a bad route: T is a perfectly good question at every σ on the axis. It is a route bad enough that the repair for bad routes fails on it.

The repair that changes the route #

Refinement treats the symptom. The cause is that the periodic wrap samples its symbol at θ = 0, where the symbol is σ, and nothing forces that. The two corner entries were chosen to make C a circulant; they could equally be chosen to make it a different matrix with a fast solve.

Put −1 in both corners with the opposite sign — +1 where C has −1 — and the matrix is anti-periodic: shifting a vector round the end flips its sign. It is still diagonalised by Fourier vectors, modulated by a half-sample phase, and its eigenvalues are the same symbol sampled at θ = (2k + 1)π/n. That grid is the periodic one moved half a step, and at d = 2 its nearest sample to the minimum is π/n away. The correction back to T is still rank two, and the solve still costs two transforms and a modulation.

The anti-periodic correction’s error stays between 2.9·10⁻¹⁵ and 6.8·10⁻¹⁴ across the whole axis, three to eight times elimination’s, and it has no trend at all as σ goes to zero. Its condition number at σ = 10⁻¹⁰ is 1,659, three per cent under T’s own — so the wrap is as well conditioned as the matrix it approximates, the capacitance matrix is well conditioned with it, and there is nothing left to amplify.

It is the circulant field’s recurring lesson in a new place. The circulant that cannot be indefinite found that which circulant stands in for a Toeplitz matrix decides whether a preconditioner is a preconditioner at all; here the same freedom decides whether a direct solve through it has any digits.

That is the practical conclusion, and it is sharper than “check the base matrix”. The circulant in a circulant-plus- correction solve is not given by the problem. It is a choice of what to put in two corners, and the choice decides which points of the symbol are sampled. A choice that samples the symbol’s zero hands the solve a singular matrix the problem never had; a choice half a step away hands it a matrix no worse than T.

What this settles, and what it measures on #

Everything here is one family, tridiag(−1, 2 + σ, −1), at three sizes, with five seeded right answers per point and the median drawn. The correction uses the fast transform for every solve against the wrap and forms the two-by-two capacitance matrix explicitly; a code that factorised S with pivoting instead of by Cramer’s rule would change the constant of a tenth and not the slope, since S is two by two and the factorisation is the same arithmetic either way.

The error constant of about 0.1 is measured, not derived, and it moves between 0.07 and 0.25 over the sizes and σ drawn. The claim is the product law and the slope of two, both of which hold at every size; the constant is what the five draws happened to give.

The cancellation ratio depends on the right-hand side. A b with no component along the constant vector — one whose end entries happen to cancel — would pass through C without being amplified, and the error would fall to κ(S)·u alone. The seeded draws are ordinary vectors, and ordinary vectors have that component.

The claim this has to refuse #

The claim is the one a reader of the cost argument would naturally make: that a solve routed through a nearby circulant and a low-rank correction is as accurate as a solve of the matrix, provided the matrix is well conditioned. Fed tridiag(−1, 2 + 10⁻⁸, −1) at n = 64, with κ = 1,712, and required to show the correction within 10⁻¹⁰ of the right answer, it must fail — and it does, at 1.2·10⁻⁴. And the anti-periodic wrap is measured beside it for the reason a refusal needs a control: without it the result could be a defect of the correction formula rather than of the wrap it was pointed at.

Whose fault is it is usually answered with two numbers, a condition number for the problem and a backward error for the algorithm. Here the problem’s condition number is small, the algorithm is exact algebra, and the fault belongs to a third party — a matrix introduced by the method for its own convenience, which neither of the two usual numbers describes. The usual pair is built to say whether an answer is the exact answer to a nearby problem and whether that problem is sensitive; neither was designed to notice a route that manufactures its own sensitivity. The number that does notice is the capacitance matrix’s condition number, and it is available for nothing, since S is formed anyway and is two by two.

Still open: which wrap, and why a twist sometimes costs nothing #

The choice of corner. Anti-periodic is right for this symbol because its zero sits at θ = 0. A symbol whose zeros sit elsewhere — an indefinite operator, d below 2, whose zeros are at ±arccos(d/2) — needs a different half-step, and the general wrap puts a phase e^(iφ) on the corners so that its samples sit at (2πk + φ)/n. Which φ to choose is then readable from the symbol, and the measurement is whether that reading is the whole of the choice.

Two near-zeros at once. Everything above had one near-null direction, the constant vector, and the rank-two capacitance matrix inherited it. An indefinite symbol has two zeros, and a real wrap can land on both at once. Whether two near-null directions against a rank-two correction behave like one near-null direction twice — or like none — is the question two near-zeros cost less than one takes up.

Rank beyond two. A pentadiagonal matrix is a circulant plus four corner entries, and a two-dimensional operator on a grid with fixed edges is a circulant along each axis plus a correction whose rank is the length of an edge. The capacitance matrix grows with it, and whether κ(S) stays at κ©/n or picks up a further factor from the correction’s own geometry is the measurement that decides whether the route is usable on the problems it is most often proposed for.