A limit the matrix never reaches

The previous essay ends on a fragility. A circulant hands over its entire spectrum in closed form because its diagonals wrap round the matrix, and the wrapping is what makes the Fourier vectors eigenvectors. Stop the diagonals wrapping and every exact claim goes at once.

What is left is a Toeplitz matrix: constant along each diagonal, so T[i][j] depends on i − j and not on i and j separately. It is described by 2n − 1 numbers rather than n, it has no zero entries either, and no basis diagonalises the family. There is no closed form for its eigenvalues at any size.

This essay is about what survives that, and the answer is: a great deal, in the limit, and nothing at a size anybody runs. Which is a more interesting answer than either extreme, and is the reason this collection has a whole verdict category called true in a limit nobody reaches.

The test matrix, and why it is this one #

The family drawn is Kac–Murdock–Szegő: T[i][j] = ρ^|i−j| for |ρ| < 1. It is the correlation matrix of a first-order autoregressive process, which is where it comes from, and it is the test matrix of this field for three reasons that are all about what comes with it rather than about convenience.

It is symmetric and positive definite, so a condition number is a ratio of eigenvalues rather than of singular values and everything below is about one object.

Its symbol is elementary. Sum the geometric series in both directions and

f(θ)  =  Σₖ ρ^|k| e^(ikθ)  =  (1 − ρ²) / (1 − 2ρ cos θ + ρ²)

which is a smooth positive function of one variable with an obvious maximum at θ = 0 and an obvious minimum at θ = π.

And its inverse is known in closed form — this is the one that matters most here, and it is startling the first time. The inverse of a dense matrix with no zero entry anywhere is tridiagonal:

K⁻¹  =  1/(1 − ρ²) · tridiag(−ρ, [1, 1+ρ², …, 1+ρ², 1], −ρ)

Three diagonals, written down rather than computed, with the two corner entries differing from the rest. assertTheClosedFormInverseIsTheInverse multiplies the two together at three values of ρ and three sizes and finds ‖KK⁻¹ − I‖/√n below 10⁻¹⁶ every time. So this field has an exact ground truth in the same sense the elimination field has the Hilbert matrix: a family of problems whose answer is known rather than approximated.

What Szegő’s theorem says, and what it does not #

The theorem is about the family rather than about any member of it. As n grows, the eigenvalues of the n×n section fill the range of the symbol: they lie strictly inside (min f, max f) and their distribution converges to that of f sampled uniformly. Both extremes are elementary because cos θ reaches ±1:

λₘᵢₙ → (1 − ρ)/(1 + ρ)      λₘₐₓ → (1 + ρ)/(1 − ρ)      κ → ((1+ρ)/(1−ρ))²

At ρ = 0.8 that is a symbol running from 0.1111 to 9.000 and a limiting condition number of 81.000.

Now read the figure. The measured condition numbers are

Every one of them is below the limit, they climb monotonically towards it, and none of them reaches it — which the theorem also says, since the eigenvalues lie strictly inside the range and the sections approach it from within.

So the sentence “the condition number of this family is 81” is true as a statement about a limit and wrong as a statement about any matrix. At n = 8 it is out by a factor of nearly two. The assertion in lib/fft.js requires the measured κ to be below the limit at every size and to climb at every size, and the refusal is fed the claim that the 16×16 section attains its own asymptotic value and required to reject it.

The gap is the measurement #

That may read as pedantry, and the drag is what makes it not. Take ρ towards one:

The limit rises steeply — it is ((1+ρ)/(1−ρ))², which goes to infinity as ρ approaches one — and the approach to it gets slower at the same time. So the case where the asymptotic value is largest and most alarming is exactly the case where a real matrix is furthest from it. At ρ = 0.95 a reader told “this family is conditioned like 1,521” and handed a 128×128 example is being told something 14% too pessimistic — and at n = 8 it would be 74% too pessimistic.

Across the seven values the slider draws — ρ = 0.3, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95 — the limit reads 3.449, 9, 16, 32.11, 81, 361 and 1,521, the 128×128 section reads 3.447, 8.988, 15.96, 31.97, 80.14, 346.6 and 1,337, and the shortfall reads 0.0%, 0.1%, 0.2%, 0.4%, 1.1%, 4.0% and 12.1%. The error in the asymptotic statement moves by a factor of a hundred and twenty across a factor of four hundred in the statement itself, and it moves the wrong way for anybody quoting it.

The symbol’s range is where that comes from and the seven readouts give it directly. f runs 0.5385–1.857 at ρ = 0.3 and 0.02564–39 at ρ = 0.95: a factor of 3.4 between its extremes at one end and a factor of 1,521 at the other, on the same interval of θ. A finite section samples n frequencies and the ones nearest θ = 0 and θ = π are about 1/n away, so what the section misses is the height f climbs over that last 1/n — small when f is flat there and enormous when it is not. The two columns of the table are the same fact: a steep symbol has a large ratio between its extremes and a large gap between the extreme and the nearest sample.

The mechanism is not subtle once the symbol is in view. The extreme eigenvalues are approaching the extreme values of f, which occur at θ = 0 and θ = π; a finite section samples n frequencies and the sampled points nearest those extremes are a distance of order 1/n away. The steeper f is near its extremes, the more that 1/n costs — and f is steep near its extremes exactly when ρ is near one. The approach and the limit are governed by the same feature of the same function.

A Toeplitz product, exactly, through a circulant twice the size #

The previous essay’s machinery is not wasted, and the way it is recovered is the neatest thing in this field.

A Toeplitz matrix is not a circulant. But it embeds in one: build the 2m × 2m circulant whose first column is the Toeplitz first column, then zeros, then the first row reversed, and T sits in its top-left block. Multiply by that circulant with a transform — which is exact, and costs O(m log m) — and the first n entries of the answer are the Toeplitz product. The rest is the wrap-around, discarded.

Nothing is approximated and nothing is truncated. assertTheEmbeddingIsExact compares the embedded product against the dense one at n = 5, 8, 13 and 32, and the worst relative disagreement is 1.1·10⁻¹⁵. The padding to a power of two is deliberate rather than incidental, and the return value says how much padding was used — at n = 13 the embedding needs 26 and pads to 32 — because fft refuses a length that is not a power of two rather than padding silently. A transform routine that padded on its own initiative would be computing the transform of a longer signal without saying so, which is a different quantity.

So a Toeplitz matrix–vector product costs O(n log n) and touches 2n − 1 numbers. That is the whole of what the structure buys computationally, and it is enough to build an iterative solver on, which is a preconditioner that changes sign.

The structure says nothing about the conditioning #

The available mistake in this field, and the one the refusal at the end of lib/fft.js is aimed at, is to hear “described by 2n − 1 numbers” as “easy”.

At ρ = 0.98 and n = 64 the matrix is described by 64 numbers, has no zero entry, and has a condition number near 10⁴. At ρ = 0.999 it is worse. The exponential decay of the entries is what makes the matrix compressible and it is also, read the other way, what makes the far corners of the matrix nearly independent of the near ones — which is precisely a small eigenvalue.

assertTheStructureAssertionsReject feeds the claim that a matrix described by few numbers is a well-conditioned matrix a 64×64 KMS at ρ = 0.98 and requires it to fail. That refusal is doing more work than most: the two properties feel connected and are not connected at all, and the reason a Toeplitz solver is worth building is that the problem is genuinely large rather than that it is genuinely easy.

Two things the exact inverse settles #

A tridiagonal inverse is not a sparse inverse in any useful sense. It is tempting to read K⁻¹ being tridiagonal as saying the problem is really a tridiagonal one wearing a disguise, and solvable in O(n) by looking at it correctly. That is true for this family and it is a statement about this family, not about Toeplitz matrices — a general symmetric positive-definite Toeplitz matrix has a completely dense inverse. What makes KMS special is that it is the correlation matrix of a Markov process, and the inverse of a correlation matrix is a matrix of conditional independences: ρ^|i−j| says the process forgets, and the tridiagonal inverse says it forgets after one step.

And it makes the whole field checkable. Every claim about a Toeplitz solve on this site is measured against kmsSolveExact, which applies the closed-form inverse. That is the difference between reporting “the fast solver agrees with the dense one” — two floating-point computations agreeing, which they might both be wrong about — and reporting a forward error against a known vector. The site has said since its foundation phase that a small residual is not a small error; a field without a known answer can only ever print the first.

Where the exactness went #

It is worth being precise about what was lost between the two essays, because “no closed form” covers several different situations and this is the mild one.

A circulant is diagonalised by a basis that does not depend on its entries. A Toeplitz matrix is not diagonalised by any fixed basis, and its eigenvectors depend on the entries in a way with no useful description. So the exact statements are gone: there is no formula for λₖ, no way to read off singularity from four numbers, and no solve by dividing.

What survives is everything that was really a statement about the symbol rather than about the sampling. The extreme eigenvalues are governed by the range of f; the distribution of the spectrum is governed by the distribution of f; and — the practical one — the circulant whose symbol agrees with f is a good approximation to T in a sense the preconditioning essay makes precise. The function was always the object and the circulant was one especially convenient sampling of it.

That is why the essay before it spends a section reframing a circulant as a sampled function before this one needs it. Read as a matrix with a formula, a circulant is a curiosity and a Toeplitz matrix is a different curiosity. Read as a function on the circle, the two are the same object under two boundary conditions, and everything asymptotic transfers.

And a circulant whose condition number is the square root of the limit #

cond2 in lib/matrix.js obtains a condition number by computing the SVD — a one-sided Jacobi sweep over n² entries. Every number in the two tables above came from it.

The symbol gives a second computation, and it does not give the same answer, which is the informative part. Build the circulant from the KMS matrix’s first column — (1, ρ, ρ², …) — and its eigenvalues are the transform of that column, which is the one-sided geometric sum

g(θ)  =  1 / (1 − ρe^(−iθ))

whose modulus runs from 1/(1 + ρ) to 1/(1 − ρ). So the condition number of that circulant is

κ  =  (1 + ρ)/(1 − ρ)

which at ρ = 0.8 is exactly 9.00 — and is exactly the square root of Szegő’s limit for the Toeplitz family. Measured at n = 8, 16 and 128 it is 9.00 to every digit at all three sizes, because the transform samples the same function however many points it uses and the extremes of |g| are attained at θ = 0 and θ = π, which are sampled at every even n.

That is worth having in view for two reasons.

A circulant does not approach anything, and a Toeplitz section does. The circulant’s condition number is a value rather than a limit, at every size, because it samples its symbol exactly. Which makes the two objects a matched pair for this essay’s argument: the same coefficients, one boundary condition that reaches its asymptotic value at n = 2 and one that does not reach it at n = 128.

And the factor of two in the exponent is the two-sidedness. The KMS matrix’s symbol is the one-sided sum in both directions, so f = |g|² up to the normalising constant, and its range is the square of |g|'s. Szegő’s ((1+ρ)/(1−ρ))² and the circulant’s (1+ρ)/(1−ρ) are one relationship, and seeing it written out is the difference between remembering the limit and knowing where it comes from.

What a Toeplitz solver has to be #

Everything above is descriptive, and the descriptions decide the method.

There is no diagonalisation, so there is no solve-by-dividing. There is an exact O(n log n) product, so a method that only needs products is affordable. The matrix is symmetric positive definite, so conjugate gradients applies. And κ has a limit, so the step count has one too — which is not obviously true and is the reason this family is a fair test rather than a rigged one.

That last point is worth spelling out, because it is where this field meets the iterative one. Conjugate gradients converges at a rate governed asymptotically by √κ. If κ grew with n the step count would grow with it and an iterative Toeplitz solver would be no better than a direct one at large sizes. It does not grow: it climbs to ((1+ρ)/(1−ρ))² and stops.

So an unpreconditioned conjugate gradient solve on this family should take a number of steps that climbs and then plateaus, and it does — 16, 24, 31, 34, 34 across n = 16 to 256 at ρ = 0.5. The plateau is not the method improving. It is the condition number reaching Szegő’s limit, and the two figures are the same figure read twice.

At ρ = 0.9 the limit is 361 and the sections have not reached it by n = 256, so the count is still climbing there — 37, 59, 94, 133 — and the plateau is off the right-hand edge. Everything about this family’s iterative behaviour is downstream of how far along that approach a given size sits, which is what makes an essay about an asymptotic condition number a practical essay rather than a pedantic one.

The direct route, which costs no accuracy here #

The other family of Toeplitz solvers — Levinson, Durbin, Schur — solves the system in n² operations rather than n³, by exploiting the structure with no transform at all. This field takes the iterative route throughout, and the usual reason for setting the direct one aside is that its stability is a delicate subject with a literature of its own. That is fair about the literature and it leaves a reader with the impression that the exponent is bought with accuracy. On the family measured here it is not.

Levinson grows the solution one row at a time, carrying a forward vector that a reflection coefficient extends. For a symmetric T the backward vector is the forward one reversed, which is what makes the recurrence two vectors rather than three. Against the closed-form answer, n = 64:

Both methods sit at a small and roughly constant fraction of κu — Levinson between 1.8% and 4.5% across four decades of conditioning — which is what a backward-stable solve looks like from the forward side. Neither is losing digits beyond what the conditioning costs. And Levinson is the more accurate of the two at every ρ up to 0.99, by a factor of two to five, at n² operations against n³.

The delicacy in the literature is real and it is about a different question. The classical recurrence can break down on an indefinite or non-symmetric Toeplitz matrix, where a leading principal minor vanishes and the algorithm divides by a zero it had no reason to expect. Every matrix in this field is symmetric positive definite, so every leading minor is positive by construction, and there the route is weakly stable and behaves accordingly.

So the honest form of the caution is narrower than the stability is delicate: the O(n²) route is safe on the class measured here and delicate outside it, which is a statement about which matrices a code will meet rather than about the algorithm. That is the same distinction a preconditioner that changes sign draws about the circulant preconditioner, where what fails is a construction’s hypothesis rather than its arithmetic.

assertTheDirectRouteCostsNoAccuracyHere measures the whole table and requires Levinson’s error to stay under the conditioning bound, to use a roughly constant fraction of it, and to be the more accurate of the two over most of the range.

What is left #

The preconditioner, which is a preconditioner that changes sign’s subject and is where the field’s finding is. Conjugate gradients on a Toeplitz system costs one embedded transform per step, so the whole question is the step count — and the standard way of controlling it is to precondition with the nearest circulant, which turns out to be indefinite at exactly the sizes anybody would run.

Non-symmetric and block Toeplitz operators, which is where the subject actually lives: a two-dimensional convolution is block-Toeplitz-with-Toeplitz-blocks, the symbol becomes a function of two variables, and Szegő’s theorem still holds with the range of that function. Nothing here reaches it, and every claim above would have to be re-measured there — which is the same shape of deferral the multigrid field made about the second dimension and which turned out to change one of its answers.

Non-symmetric and block Toeplitz operators remain the largest gap — see above.

A limit that is a defective matrix #

A family approaching a limit it never reaches is the shape of several results on this site. The sharpest is a matrix walking up to a defective one, where the answer has a limit and the standard formula for it does not.