Numerical Linear Algebra
where the algebra stops being the arithmetic
Iterating, instead of factorising

The spectrum that predicts nothing

For a symmetric matrix the eigenvalues govern how fast an iteration converges. Drop symmetry and they stop governing anything — there is a matrix whose eigenvalues are as evenly spread as eigenvalues can be, on which GMRES makes no progress at all until the last possible step.
15 min read 7 figures Where the algebra stops being the arithmeticAssertions that rejectMeasured, not assertedWhose fault is it

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.

GMRES on the Laplacian and on the cyclic shift, both 12×12A semi-logarithmic plot of relative residual against iteration. One curve falls steadily; the other is flat at one for every step until the last, where it drops to zero.02468101210⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1iteration‖r‖ / ‖b‖Laplaciancyclic shiftno progress at allevery eigenvalue of the shift is on the unit circleand it predicts nothing
Fig. 1 Two GMRES runs. One is the model problem, converging steadily. The other is the cyclic shift matrix, whose residual is exactly 1 at every step until the last, where it drops to zero. Every eigenvalue of that matrix lies on the unit circle. Drag the size and the stall is always exactly n − 1 steps.

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_k = Vk+1HkV_{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.

How far a perturbation of size ε moves an eigenvalue, 8×8A log–log plot of eigenvalue movement against perturbation size. The symmetric case lies on a line of slope one; the non-symmetric case lies on a line of slope one eighth, and at a perturbation of ten to the minus sixteen it has already moved by a hundredth.10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/8)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision
Fig. 2 What non-normality costs in the eigenvalue problem, which is the same underlying fact. A symmetric matrix’s spectrum is stable under perturbation; a defective matrix’s moves by a fractional power. In both cases the eigenvalues of a non-normal matrix are a less informative object than they appear.

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: there are documented problems on which GMRES(10) converges and GMRES(20) stagnates. That is a genuinely uncomfortable fact about a parameter most codes leave at a default.

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 previous three essays

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.

Conjugate gradients at κ = 104, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417
Fig. 3 The symmetric case, for comparison. A rate, a bound, a prediction from a single number computable in advance — none of which has a counterpart on the other side of the hypothesis.
The residual basis of conjugate gradients, at κ = 106A semi-logarithmic plot against iteration count showing the loss of orthogonality among the residual vectors rising while the relative residual falls.05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here
Fig. 4 And the failure that does carry over, worse rather than better. GMRES orthogonalises explicitly against a stored basis, so what conjugate gradients lose implicitly it can at least measure — and must, because the ghost eigenvalues that follow from an unmeasured drift are indistinguishable from real ones in the answer.
Least squares by QR and by the normal equations in binary32Relative error of the computed coefficients against epsilon, on log axes. The QR route is a flat line near the bottom; the normal-equations route climbs and then stops, at the epsilon where the cross-product matrix becomes exactly singular.10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴
Fig. 5 The cost of the route back to symmetry. Forming AᵀA squares the condition number, and since the convergence rate goes as its square root, the square root of the square is the original — which is exactly the factor symmetry was supposed to buy.
Householder reduction to Hessenberg form, on a symmetric matrixTwo matrices side by side. The first is full; the second has zeros everywhere below the first subdiagonal, and being symmetric is tridiagonal.-0.16-0.63-0.493.43-0.025-0.632.3-0.69-0.95-1.70.076-0.49-0.694.3-1.6-1.41.93.4-0.95-1.65.40.0141.63-1.7-1.40.0144.60.055-0.0250.0761.91.60.0553.1A, symmetric→-0.164.600004.65.92.500002.51-1.20000-1.24.90.3800000.385.80.2400000.242H = QᵀAQ, tridiagonal‖A − QHQᵀ‖/‖A‖1.1·10⁻¹⁵below the subdiagonal0worst eigenvalue movement7.1·10⁻¹⁵a similarity, so the spectrum is untouched — and every later step is O(n²) rather than O(n³)one reduction, then every iteration is cheapthe eigenvalues did not move
Fig. 6 The form Arnoldi produces whether or not anybody asked for it. One nonzero subdiagonal is what turns the least-squares problem inside GMRES into a sequence of Givens rotations, and it is the same structure the eigenvalue algorithm reduces to deliberately.
Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 7 And the failure GMRES inherits in full. Arnoldi is modified Gram–Schmidt against a stored basis, so it has every opportunity to lose orthogonality and every reason to re-orthogonalise — which is why serious implementations do and conjugate gradients can afford not to.

Where this goes next

0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417 The rate the condition number predicts Conjugate gradients converge at a rate governed by the square root of the condition number. That is a bound rather than an estimate, it is provable, and it is loose enough that provisioning iterations from it wastes nine out of ten. 10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/8)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision Symmetry is worth more than precision A symmetric matrix gives up its eigenvalues to full accuracy however ill-conditioned it is. An unsymmetric one can move them by the eighth root of a perturbation, so the rounding involved in merely storing the matrix shifts the spectrum by a hundredth. 05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here An orthogonalisation nobody calls one Conjugate gradients are derived as a minimisation and behave as an orthogonalisation, which is why the finite-termination property in every textbook is not a property the method has in floating point.