The matrix that is one row

Every matrix on this site so far has been n² independent numbers, and every cost claim has followed from that. Elimination is n³ because there are n² entries to eliminate against. The sparsity field is entirely about what happens when most of those n² numbers are zero: the matrix is smaller than its shape suggests, and the argument becomes about which of the zeros survive into the factors.

This field is about the other way a matrix can be small, and the sparsity field’s machinery does not apply to it at all.

A circulant is a matrix each of whose rows is the previous one shifted right by one place. Write down its first column and the whole matrix follows. At size 16 that is sixteen numbers describing 256 entries — and not one of those entries is zero. Every ordering argument, every fill count, every elimination-graph statement in the sparsity field is inapplicable by construction. There are no zeros to preserve.

The eigenvectors are known before the matrix is #

The claim is stronger than there is a formula, and it is worth stating precisely because the strength of it is what the rest of the field is built on.

Let S be the cyclic shift — the matrix that moves each entry down one position and wraps the last one round. Every circulant is a polynomial in S:

C  =  c₀I + c₁S + c₂S² + … + c_{n−1}S^{n−1}

S has an obvious set of eigenvectors: the vectors whose entries are the powers of an n-th root of unity, since shifting such a vector multiplies it by that root. So every circulant of size n has the same eigenvectors, and they do not depend on the entries in any way. Two circulants of size 16 built from completely different numbers are diagonalised by the same 16 vectors, which are known before either matrix is written down.

What is left for the entries to decide is the eigenvalues, and they are the polynomial above evaluated at the n roots of unity — which is exactly the definition of the discrete Fourier transform of the vector c.

For the first column drawn above, whose entries are 1, ½, ⅕, ⅒ and so on down to 0.0154 and back up again, the transform gives

2.9040  2.1633  1.4373  0.9937  0.6790  0.4863  0.3582  0.2952
0.2702  0.2952  0.3582  0.4863  0.6790  0.9937  1.4373  2.1633

and the eigensolver, given the assembled 16×16 matrix and told nothing about its structure, returns the same sixteen numbers to a worst relative disagreement of 1.2·10⁻¹⁵.

The same matrix, read as a function #

There is a second way to hold a circulant that makes the next essay possible, and it is worth setting up here where the algebra is exact.

The eigenvalues above are a transform, and a transform is a function evaluated at n points. Define

f(θ)  =  c₀ + c₁e^(iθ) + c₂e^(2iθ) + … + c_{n−1}e^{(n−1)iθ}

— the symbol of the matrix. The eigenvalues are f evaluated at the n roots of unity, θ = 2πk/n. So a circulant is not really a matrix with a formula attached; it is a function on the circle, sampled. Two circulants of different sizes built from the same coefficient sequence are two samplings of the same function, and everything that is true of them in the limit is a statement about f.

That reframing does three things at once. It says why the eigenvectors do not depend on the entries — they are the sampling points, and the entries are the function. It says immediately which circulants are singular, since a zero eigenvalue is f vanishing at a sampled point. And it is the only reason anything can be said about the Toeplitz matrices of the next essay, which have no exact spectrum at all and whose extreme eigenvalues are governed, in the limit, by the range of exactly this f.

A circulant multiply is a convolution #

The other thing worth noticing before leaving the exact case: multiplying by a circulant is the same operation as convolving cyclically with its first column. That is not an analogy — write out the i-th entry of Cx and it is Σⱼ c_{i−j mod n} xⱼ, which is the definition.

So the convolution theorem and the diagonalisation above are one fact stated twice. Convolution becomes multiplication under the Fourier transform and a circulant is diagonalised by the Fourier matrix say the same thing about the same object, and a reader who has met one has met the other.

It also explains where circulants come from. Nobody assembles one by writing down a first column. They arise wherever a computation applies the same local rule at every position of a periodic domain — a blur, a difference stencil on a ring, a filter — and the periodicity is the assumption that turns a Toeplitz operator into a circulant one. Which is why the field’s second essay is about what happens when the domain is not periodic, and the answer is that most of this survives approximately and none of it survives exactly.

Why that comparison is worth making #

Because the two routes share nothing. One is a transform of sixteen numbers, log₂16 = 4 stages of butterflies, about 64 multiply-adds. The other is jacobiEigSym sweeping over 256 entries, annihilating off-diagonal pairs by rotation until the off-diagonal norm falls below a tolerance, which on this matrix takes several sweeps and several thousand operations.

This site’s two routes to a number habit exists because a single computation confirms only that it is self-consistent. Here the two routes are about as far apart as two computations of the same quantity can be, and the agreement is at the level of the arithmetic rather than at the level of a tolerance somebody chose.

It is also the site’s exact ground truth thread reaching somewhere it usually cannot. The Hilbert matrix has an exact rational inverse, so a Hilbert solve has a known answer; that is one family. This is a whole class of matrices whose spectrum is available in closed form for any entries at all, which is unusual enough in this subject to be worth a field.

The factorisation nobody forms #

C = F*ΛF is a factorisation. It is not a metaphor for one and it is not an approximation: F is the Fourier matrix, Λ is the diagonal of eigenvalues, and the product is exactly C.

So the rule this site is named for applies to it, and the figure above carries a residual badge like every other factorisation figure here. Solving Cx = b through the transform — transform b, divide by the eigenvalues, transform back — gives a solution whose relative residual is 4.7·10⁻¹⁶, which is the same quantity an elimination’s badge prints and is of the same size.

What is different is that no factor was ever formed. An LU factorisation stores n² numbers; a QR stores two matrices. This one stores the n eigenvalues and refers to F, which is not stored anywhere because it is the same F for every circulant of that size. The residualcheck gate lists circulantSolve alongside luFactor and qrHouseholder in its factoriser list for exactly this reason: a decomposition that costs nothing to store is still a decomposition, and a figure drawing one is still obliged to print what it cost in accuracy.

What is real is asserted, and what is real is checked #

Two smaller claims sit on the figure and both would fail silently without a check.

The spectrum is real. A symmetric circulant — one whose first column reads the same forwards as backwards from the second entry on — has real eigenvalues, and the transform of such a column has zero imaginary part. That is a statement about the arithmetic as well as the algebra: a sign error in one butterfly of the fast transform would leave a small imaginary part and nothing else would show, because the real parts would still be close to right and the picture would look identical. The measured imaginary part is 9.7·10⁻¹⁷, and the assertion requires it below 10⁻¹² times the scale.

The matrix is not singular. circulantSolve refuses a circulant whose smallest eigenvalue is at the level of rounding rather than returning infinities, and this is the one place in the field where that refusal is informative rather than defensive. A circulant is singular exactly when its symbol vanishes at a root of unity — the vector (1, −1, 1, −1) has a first column summing to zero, so λ₀ = 0 exactly, and the matrix is singular for a reason that is visible in four numbers before anything is computed. Every other class of matrix on this site has to be factorised before that can be known.

The transform is written out, and checked against the definition #

lib/fft.js contains a radix-2 Cooley–Tukey transform in about forty lines: bit-reverse the input, then log₂n stages of butterflies. It is written out rather than imported, for the same reason matrix.js writes out LU and QR — the content of this site is how these routines behave, and a library call hides it.

Forty lines of index arithmetic is also forty lines in which an off-by-one produces a plausible answer, so it is compared against a direct transform written from the definition,

X[k] = Σⱼ x[j] · e^(−2πi jk/n)

which is O(n²), shares no line of code with the fast version, and exists on this site for no other purpose. Over random complex input at n = 2, 4, 8, 16 and 64 the two agree to a worst relative difference of 1.5·10⁻¹⁴, and the round trip through the inverse returns the input to 3.4·10⁻¹⁵. Nothing else in this field is believed until that comparison passes.

And a second reason to prefer it that nobody mentions #

The fast transform is the fast one. It is also, on the same input, the more accurate one, and that has nothing to do with speed.

The direct transform sums n terms in index order. The order they are added in establishes what that costs: a sum taken sequentially accumulates error like √n in the best case and like n when the terms are biased, because each partial sum is rounded and the roundings compound down a chain of length n. The fast transform is not a chain. It is a binary tree of depth log₂n, so the error compounds down a path of length four rather than sixteen at n = 16, and like √(log n) rather than √n as the size grows.

That is the summation field’s argument applied to a transform, and it is the reason a numerical analyst would prefer the fast transform even on a machine where both took the same time. It almost never gets said, because the cost argument is so overwhelming that nobody looks for a second one.

What the saving actually is #

The cost claim is usually written as “n log n instead of n³” and left there, which understates it by enough to be worth writing out. At the sizes this site’s other fields work at:

The ratio itself grows like n²/log n, so the two methods are not competing algorithms for the same job at any interesting size — they are a computation that happens and a computation that does not. The bottom row is a matrix of about four billion entries, which cannot be stored; the transform of the 65,536 numbers describing it is a fraction of a second.

And storage moves the same way and further. The elimination needs n² numbers for the factors, which is 4.3·10⁹ at the bottom row. The transform needs the n numbers it started with.

None of that is a statement about this site’s arithmetic, which counts operations rather than seconds — see the cost field on why a duration is not a quantity this collection reports. It is arithmetic about arithmetic, and it is the reason structure is worth detecting even when detecting it is awkward.

What breaks it #

The structure is fragile in a way worth stating, because it is what makes the rest of the field necessary.

Change one entry and the matrix is no longer a circulant. Not approximately not a circulant — its eigenvectors are no longer the Fourier vectors, its eigenvalues are no longer a transform of anything, and the closed form is gone entirely. A matrix that is a circulant plus a rank-one correction is a perfectly ordinary object with a perfectly ordinary spectrum, recoverable only through the correction formulae or by factorising.

Stop the diagonals wrapping and the same thing happens. The wrap is what makes S a shift with Fourier eigenvectors; without it S is nilpotent, its only eigenvalue is zero, and no basis diagonalises the family.

That fragility is why the field is three essays rather than one. If circulants were robust to small changes there would be nothing more to say. They are not, and what happens to each of the exact claims above under the mildest possible perturbation of the structure — losing the wrap — turns out to be a different answer for each of them: the spectrum becomes approximate, the condition number gets a limit rather than a value, and the fast solve becomes a preconditioned iteration whose behaviour depends on a sign nobody mentions.

The refusal #

The claim this essay is fed, and must refuse, is that a dense matrix costs n³ to factorise.

It does if the only thing known about it is that it is dense. rejects() is handed a circulant whose first column is (1, −1, 1, −1) and required to refuse it — a matrix with no zero entries, no obvious structure to the eye, and a spectrum available in four operations, one of which is zero. A routine that returned an answer for it would be a routine that had not looked at what it was given.

One more thing the exactness buys #

A closed-form spectrum is usually presented as a way of avoiding a computation. On this site it is worth more as a way of testing one, and the figure above is the demonstration.

jacobiEigSym is written out in lib/matrix.js and is used by thirty-odd figures across four fields. Nothing in the site can check it against a better eigensolver, because there is no better eigensolver here; the site’s own is the only one. What the site’s own routines are usually checked against is reconstruction — does VΛVᵀ come back to A — and reconstruction is a real test that misses one whole class of error. A routine that returned the right eigenvalues in the wrong order, or that converged to a nearby symmetric matrix rather than the one it was given, would reconstruct beautifully.

A circulant closes that gap. The eigenvalues are known independently, so an eigensolver run on one is being compared against an answer rather than against itself. That the comparison passes at 1.2·10⁻¹⁵ is a fact about jacobiEigSym first and about circulants second, and it is the sort of thing this collection is built to notice: a family of matrices whose only unusual property is that they can report whether the instrument works.

What is left #

The structure that does not wrap. A circulant’s diagonals wrap round the matrix, which is what makes the Fourier basis diagonalise it and is a strong assumption about the problem. Drop the wrapping and the matrix is Toeplitz — constant along each diagonal, 2n − 1 numbers instead of n — and no transform diagonalises it. The next essay is about what survives that loss, which is more than seems likely and less than would be convenient.

The conditioning, which the structure says nothing about. Everything above is a statement about the spectrum being available, and none of it is a statement about the spectrum being good. A circulant can be arbitrarily ill conditioned and its condition number is as available as everything else: for the matrix drawn here it is 10.75, which is the ratio of the largest eigenvalue to the smallest and took one transform to obtain.

And the cost of the transform itself, which is n log n multiply-adds and a memory access pattern that is nothing like a matrix multiply’s. The cost field is about exactly that distinction and does not reach the transform; the bit-reversal permutation at the start of every fast transform is one of the least cache-friendly operations in numerical computing, and it is a subject rather than a footnote.