The spectrum that predicts nothing

The three essays before this one are about conjugate gradients, and every quantitative claim in them depends on a hypothesis stated once and then relied on throughout: the matrix is symmetric and positive definite.

That hypothesis is doing far more work than it appears to. It supplies real eigenvalues, so a spectrum can be summarised by its two ends. It supplies a norm in which the method is optimal, so convergence can be bounded. And it supplies a short recurrence, so each step costs the same as the last.

Drop it, and all three go. The first two go in a way that is genuinely surprising, and this essay is about that.

“Always exactly” is a strong enough claim to be worth the slider rather than the sentence, because it is the kind of statement that is usually true of the sizes somebody looked at.

The word doing the work is exactly. Not 0.997, not falling by a per cent a step — the residual is the same number it started at, and then it is zero.

Five sizes, and the stall is n − 1 at every one of them. Six, fourteen, seventeen and twenty give five, thirteen, sixteen and nineteen steps of no progress whatever, and the twelve of the hero gives eleven. There is no size at which the method gets an early foothold and none at which it does worse than the theoretical worst.

That is what makes this matrix the counterexample rather than an anecdote. Its eigenvalues are the n-th roots of unity — all on the unit circle, all simple, evenly spaced, and by every spectral reading a benign spectrum. A convergence estimate built from that spectrum predicts steady progress at every n. The measurement is zero progress at every n, for n − 1 steps, exactly.

And the residual is exactly 1 rather than nearly 1, which rules out the reading that the method is converging slowly. Slow convergence is a rate; this is the absence of one. Until the Krylov space has dimension n it contains no component of the answer at all, so there is nothing for a least-squares step to reduce, and the first moment it does contain one the residual goes to zero in a single step.

What GMRES does, and the one thing it guarantees #

The construction is the natural generalisation. Build the Krylov subspace as before — b, Ab, A²b, and so on — and at step k return the vector in that subspace which minimises the ordinary 2-norm of the residual ‖b − Ax‖.

One consequence is immediate and it is the method’s whole reliability story. The residual cannot increase. The subspace at step k + 1 contains the subspace at step k, so the minimum over the larger set is no larger than the minimum over the smaller. Monotonicity is not a theorem about convergence; it is a property of minimising over a growing set, and it holds for every matrix, at every step, with no hypothesis at all.

That is asserted on every frame of the figure: at each of up to twenty steps, the residual is checked against the one before. It is the only property in this essay that survives dropping symmetry, and it is worth knowing that it survives because it never depended on anything.

What it costs is the recurrence. Minimising over the whole subspace requires an orthonormal basis for the whole subspace, so GMRES runs Arnoldi — modified Gram–Schmidt against every previous basis vector — and step k costs k inner products. The work per step grows without bound and so does the storage. That is why real implementations restart: run m steps, take the answer, throw the basis away, start again. Restarting caps the cost and loses the optimality, and a restarted run can stagnate where the full one would have converged.

How the residual is known without computing the answer #

Worth a paragraph of machinery, because it explains how the figure can plot a residual at every step cheaply, and because the mechanism reappears elsewhere on this site.

Arnoldi produces, alongside the orthonormal basis V, a small matrix H of shape (k+1)×k that is upper Hessenberg — zero everywhere below the first subdiagonal. The relation AVₖ = $V_{k+1}H_k$ holds exactly, and it converts the minimisation over the Krylov subspace into a least-squares problem of size (k+1)×k, which is tiny.

That small problem is solved by Givens rotations applied one column at a time, each rotation zeroing the single subdiagonal entry the new column brought. The rotations are accumulated, so adding a step costs one rotation rather than a fresh factorisation — and the residual norm falls out as the magnitude of the last entry of the rotated right-hand side, available before the solution vector is ever assembled.

So GMRES knows its residual at every step for free, and computes x only when the run stops. The Hessenberg form doing the work here is the same structure the form that makes it affordable is about in the eigenvalue problem: one subdiagonal is what makes an O(n³) step into an O(n²) one, in both places, for the same reason.

The matrix that settles it #

The cyclic shift: the permutation that moves each coordinate to the next and wraps the last around to the first. Its matrix is all zeros except for a single 1 in each row and column.

Consider what it does to the Krylov subspace generated by e₁, the first coordinate vector. Se₁ = e₂. S²e₁ = e₃. The subspace after k steps is spanned by the first k coordinate vectors, and it does not contain e₁ back again until step n.

So the best approximation to the solution available at step k is, for every k below n, no better than the zero vector. The residual is exactly 1 — checked to 10⁻⁹, at every step, at every size on the slider — and then at step n the space closes up and the answer is exact.

Now the spectrum. S is an orthogonal matrix, which is asserted on the figure rather than assumed: ‖SᵀS − I‖ is checked below 10⁻¹³. Every eigenvalue of an orthogonal matrix has modulus one, and these are precisely the n-th roots of unity — real and imaginary parts spread evenly around the circle, none clustered, none near zero, condition number exactly 1.

A perfectly conditioned, beautifully distributed spectrum, and an iteration that achieves nothing for n − 1 steps.

Why this is not a curiosity #

It would be comfortable to file the shift matrix as a pathological example, and the theorem that covers it forbids that.

Greenbaum, Pták and Strakoš proved that any non-increasing residual curve can be realised by a matrix with any prescribed spectrum. Pick a convergence history — stagnate for fifty steps, then plunge; descend smoothly; anything non-increasing — and pick a set of eigenvalues. There is a matrix with those eigenvalues on which GMRES produces exactly that history.

That is not a statement about rare matrices. It says the spectrum and the convergence behaviour are independent: knowing one places no constraint whatsoever on the other. Every heuristic of the form “cluster the eigenvalues and GMRES will converge quickly” is, as a general claim, false — and the site’s refusal for this figure is exactly that sentence, fed the shift matrix’s own numbers and required to throw.

What is actually governing it #

If not the eigenvalues, then what — and the honest answer is that there is no single replacement, which is why the non-symmetric case is a harder field than the symmetric one.

The missing ingredient has a name: normality. A matrix is normal when AᵀA = AAᵀ, equivalently when it has a complete orthonormal set of eigenvectors. Symmetric matrices are normal. For a normal matrix the eigenvalues do govern GMRES convergence, and the whole intuition holds.

The shift matrix is orthogonal and therefore also normal, which is why it is such a sharp example: it is not even a departure from normality that produces the stagnation. What produces it is that its eigenvectors, though orthogonal, are the Fourier modes — and e₁ has equal components along all of them, so no small set of eigenvalues can be annihilated by a low-degree polynomial. The right-hand side matters, which is a dependency conjugate gradients’ bound does not have.

For a non-normal matrix there is a further problem. The eigenvector matrix can be badly conditioned, and everything derived from a diagonalisation then carries that condition number as a factor. This is the same fact symmetry is worth more than precision measures from the other side: a symmetric matrix’s eigenvalues move by at most the size of a perturbation, and a defective matrix’s can move by its eighth root.

The pseudospectrum, in one paragraph #

The tool that replaces the spectrum, since something has to.

The ε-pseudospectrum of A is the set of numbers that are eigenvalues of some matrix within ε of A. For a normal matrix it is just the eigenvalues fattened by ε, and it says nothing new. For a non-normal one it can be enormous — a matrix whose eigenvalues are all at 0 can have a pseudospectrum, at ε = 10⁻¹⁰, that reaches out to modulus 1.

That is the object convergence bounds for non-symmetric Krylov methods are stated in terms of, and it is the honest reason this essay does not offer a rate. There is no formula here of the shape the rate the condition number predicts provides, because for the non-symmetric case the corresponding quantity is a region in the complex plane rather than a number, and computing it costs far more than the solve.

Restarting, and the stagnation it manufactures #

The growing cost per step is not a theoretical inconvenience; it is why the method as described is almost never what runs.

GMRES(m) runs m steps, forms the solution, discards the basis, and starts again from the new residual. Storage is capped at m vectors and the cost per step is capped with it. What is lost is the property the method was built on: the minimisation is now over a subspace that has forgotten everything before the restart, so the guarantee of optimality applies only within the current cycle.

The consequence is a failure mode with no analogue in the full method. GMRES(m) can stagnate completely — each cycle reducing the residual by a factor so close to 1 that the run never terminates — on a problem the full method solves in 2m steps. The residual is still monotone, because each cycle starts from the previous residual and can only improve it, so the plot looks like a method that is working slowly rather than one that has stopped.

And the choice of m is not monotone in the way anybody expects. Increasing the restart length usually helps and sometimes hurts, and that is a genuinely uncomfortable fact about a parameter most codes leave at a default — so it is worth measuring rather than citing.

A cyclic shift on twenty unknowns with the wrap-around entry set to θ, so θ = 1 is the shift whose spectrum predicts nothing and smaller θ pulls it in. Total inner steps to a relative residual of 10⁻⁸, stall meaning four hundred cycles without one:

A longer restart costs more work at four of the six θ. At θ = 0.7, m = 9 takes 738 steps and m = 10 takes 860 — seventeen per cent more for storing one extra vector — and m = 12 takes 552 where m = 14 takes 728, thirty-two per cent more. These are not marginal reversals: a parameter left at a default, moved in the direction everyone assumes is safe, costing a third of the run.

And the stagnation threshold is a number that moves with the problem. Every m below 8 stalls at θ = 0.9, below 4 at θ = 0.8, below 3 at θ = 0.7, and nothing stalls at or below θ = 0.6. There is a minimum restart length, it is not knowable from the matrix’s size, and below it the residual is still monotone — the plot shows a method working slowly rather than one that has stopped.

What this sweep does not find is the strongest form of the claim: no row here stalls at a larger m having converged at a smaller one. On this family the stagnation threshold is monotone even though the work is not, and that is worth saying because it is the weaker of the two claims and the one the measurement supports.

assertALongerRestartIsSometimesMoreWork runs the sweep and requires at least one reversal, the unrestarted method to converge on the same problem, and a long enough restart to converge — so the threshold is a restart length rather than a property of the matrix alone.

The short-recurrence alternatives, and what they trade #

Since the growing cost is the problem, the obvious question is whether the short recurrence can be recovered without symmetry, and there is a family of methods that answer yes at a price.

BiCG, CGS, BiCGSTAB and QMR all build a Krylov basis with a fixed-length recurrence, so the cost per step is constant and the storage is a handful of vectors. What they give up is the minimisation: none of them minimises the residual over the subspace, and consequently none of them has a monotone residual. A BiCGSTAB run can and does show the residual jumping up by orders of magnitude before coming down again, and it can break down outright — a division by a quantity that has become zero for reasons unconnected to convergence.

So the non-symmetric case offers a choice between a method whose cost grows and whose residual is guaranteed, and methods whose cost is fixed and whose residual is not. The symmetric case offers both at once, which is what symmetry is worth more than precision is about and is the strongest available argument for symmetrising a problem when the option exists.

What is asserted here #

The residual never rises, at every step of the model problem’s run and at every size on the slider. The one guarantee that survives is checked rather than cited.

The shift makes no progress, with its residual equal to 1 to within 10⁻⁹ at every step from 1 to n − 1.

And then solves exactly at step n, below 10⁻¹².

The shift is orthogonal, ‖SᵀS − I‖ below 10⁻¹³ — which is what entitles the essay to say its eigenvalues are on the unit circle without computing them.

The refusal: the claim that a well-spread spectrum implies fast GMRES convergence is fed the shift’s residual at step n/2 and must throw. It does. That refusal is the essay.

The tempting way back to symmetry, and what it costs #

There is an obvious route out of all of this, and it is worth following far enough to see the bill.

AᵀA is symmetric positive definite for any A of full column rank. So conjugate gradients can be applied to AᵀAx = Aᵀb, with all the guarantees restored: a short recurrence, a rate governed by a condition number, monotone error in a norm. The method is called CGNR and it takes four lines to write given a CG implementation.

The bill is on this site already. κ(AᵀA) = κ(A)², which is the road that squares the problem — and since the convergence rate goes as the square root of the condition number, the square root of the square is the original. The iteration count for CGNR is governed by κ(A) rather than √κ(A), which is precisely the factor symmetry was supposed to buy.

So the trade is exact and it is not free: symmetry is recovered, the guarantees come back, and the rate reverts to what a non-symmetric method would have needed anyway. Worse, forming AᵀA explicitly loses accuracy at the small end of the spectrum in the manner rank is a decision measures, so the usual advice is to apply the products implicitly — one multiplication by A and one by Aᵀ per step — which is what LSQR does, and LSQR is CGNR with the arithmetic arranged so the square is never formed.

That is the honest summary of the non-symmetric situation. Every route out of it pays somewhere: in storage with full GMRES, in guarantees with a restart, in monotonicity with the short-recurrence family, or in the rate with the normal equations.

What this changes about the essays before it #

Not their correctness — their scope, and the scope was stated in each of them as a hypothesis that looked like a technical precondition.

The convergence bound in the rate the condition number predicts is a symmetric-case result. Its loose factor of ten, its dependence on √κ, its behaviour under the slider: all of it lives inside the assumption that the eigenvalues are real and the eigenvectors orthogonal.

Preconditioning as changing κ, in changing the condition number on purpose, is also symmetric-case. Preconditioning non-symmetric systems is still the right thing to do and still helps enormously in practice, but “design the preconditioner to reduce κ” has no comparable justification, because κ does not control the rate.

The loss of orthogonality in an orthogonalisation nobody calls one is the one that transfers directly, and it transfers worse. Arnoldi is Gram–Schmidt written out, and GMRES stores the whole basis, so it has every opportunity to lose orthogonality and every reason to care — which is why serious implementations reorthogonalise there and can afford not to bother in conjugate gradients.

The pattern is the site’s own, arriving one level up. A derivation that assumes an invariant cannot see the cost of the invariant failing; a field built on a hypothesis cannot see which of its results were consequences of the hypothesis until something drops it. The shift matrix is what dropping it looks like.

The same failure, for an exponential #

A non-normal matrix’s spectrum does not predict how a Krylov method behaves on it. It does not predict the size of the matrix’s exponential either, and there the discrepancy has a closed form.