Symmetry is worth more than precision

Take the eight-by-eight matrix with ones on the superdiagonal and zeros everywhere else. It is nilpotent: all eight of its eigenvalues are zero.

Now put ε in the bottom-left corner. Its characteristic polynomial becomes λ⁸ − ε exactly, so its eigenvalues are the eighth roots of ε and every one of them has modulus $\varepsilon^{1/8}$.

Set ε = 10⁻¹⁶ — the size of the rounding involved in storing a number, not in computing with it. The eigenvalues now have modulus 10⁻². A perturbation at the level of representation has moved the spectrum by a hundredth.

Slope one against slope one eighth is the reading at eight unknowns, and the eighth is the size rather than a constant.

Six is the smallest block the figure will draw, for a reason the last paragraph of this section gives. One more unknown:

The exponent is the size, so every extra unknown makes the non-symmetric case worse at rounding — which is the opposite of how a reader expects a conditioning result to behave, since a larger symmetric problem is no harder here at all.

The movement at rounding is exactly 10⁻¹⁶ᐟⁿ, to every digit printed, at all five sizes. The slope is 1/n because a perturbation of size ε in the corner of an n×n Jordan block moves its eigenvalues by ε¹ᐟⁿ, and the closed form is verified against the characteristic polynomial at two values of ε rather than plotted and trusted. So “slope one eighth” is a fact about eight unknowns and the general statement is worse than it: the larger the block, the flatter the line, and the more a rounding-level perturbation costs.

The symmetric column does not move with n at all — 1.55, 4.44, 3.55, 4.44 and 6.22 × 10⁻¹⁵, no trend across a factor of 1.7 in the size, and every one of them at the eigensolver’s own rounding. Weyl’s bound has no n in it and the measurement has no n in it either.

So the gap at ε = 10⁻¹⁶ runs from 1.4·10¹² to 4.0·10¹² across these five sizes, and it widens with n because only one of the two columns notices n. That is the title as an inequality between two exponents rather than between two numbers: symmetry buys a slope, and no amount of precision buys a slope.

One boundary worth stating, because the figure enforces it. At n = 5 the movement at rounding is 6.3·10⁻⁴, and the generator refuses to draw — it asserts that a rounding-level perturbation moves the non-symmetric spectrum by more than a thousandth, and at five unknowns it does not. The effect this essay is about is a large-n effect, the assertion says where it starts, and the smallest block on which the claim survives is six.

The easy case is very easy #

For a symmetric matrix, Weyl’s inequality says that if A is perturbed by E then every eigenvalue moves by at most ‖E‖₂. Not “at most a constant times ‖E‖”, not “asymptotically” — at most ‖E‖, and the bound is attained.

That is an unusually clean statement for this subject and it has a consequence worth naming: the eigenvalue problem for a symmetric matrix is perfectly conditioned in the absolute sense, regardless of the condition number of the matrix. A symmetric matrix with κ = 10¹⁵ still surrenders its eigenvalues with absolute error at the level of ‖A‖·u.

The relative accuracy of the small eigenvalues is a separate question and a harder one — an absolute error of ‖A‖u is a relative error of ‖A‖u/λ, which for a tiny λ is large — and it is the question the best approximation there is runs into. But absolute accuracy is guaranteed, and for the symmetric case the guarantee is essentially free.

The hard case is arbitrarily hard #

For an unsymmetric matrix nothing of the sort holds. The relevant result is Bauer–Fike: an eigenvalue of A + E lies within ‖E‖·κ(V) of an eigenvalue of A, where V is the matrix of eigenvectors. When the eigenvectors are nearly parallel, κ(V) is enormous, and when the matrix is defective — when it does not have a full set of eigenvectors — κ(V) is infinite and the bound says nothing.

The Jordan block is the defective case in its purest form. Its eight eigenvalues are all zero and it has one eigenvector. There is no basis of eigenvectors to be badly conditioned; there is no basis at all.

What replaces the bound is the fractional-power behaviour: a perturbation of size ε moves the eigenvalues by $\varepsilon^{1/n}$ for an n×n Jordan block. That is not a pathology of the perturbation being placed in the corner — a random perturbation gives the same rate — and it gets worse as the matrix gets larger, which is the opposite of intuition.

The ground truth here is a theorem #

The Jordan curve in the figure is not measured, and that is deliberate.

Computing the eigenvalues of a defective matrix numerically is precisely the thing that is ill-conditioned, so a measured curve would be measuring the solver’s own error as much as the phenomenon. Instead the closed form is used: λ⁸ = ε, so $|\lambda| = \varepsilon^{1/8}$, exactly.

And the closed form is checked. At two values of ε the claimed root is substituted into the characteristic polynomial — computed by cofactor expansion of λI − M, which shares no arithmetic with taking an eighth root — and required to give zero to twelve digits. That is two routes to a number in the one form available for a quantity no solver can be trusted with.

The symmetric curve is measured, over six seeded perturbations at each of thirty-one sizes, because there the solver can be trusted and a measurement is the stronger evidence.

The measurement has a floor, and the assertion says so #

Weyl’s bound is an inequality with no constant, so the natural assertion is that the measured movement never exceeds ‖E‖. Below ε ≈ 10⁻¹⁴ that assertion fails, and the reason is not that Weyl is wrong.

At those perturbation sizes the movement being measured is the eigensolver’s own rounding rather than the perturbation’s effect. The Jacobi method returns eigenvalues accurate to about ‖A‖u, so comparing two runs with a 10⁻¹⁶ perturbation between them compares two numbers whose difference is dominated by their own noise.

The assertion therefore carries the measurement floor explicitly: the movement must be below ‖E‖ + 20u‖A‖. That is not a fudge to make the test pass — it is the correct statement of what the measurement can support, and writing it down is the difference between a check that means something and one that has been tuned until it stopped complaining. The alternative, restricting the range to where the bound holds cleanly, would have hidden the floor rather than reporting it.

What Jacobi buys, and why it is used here #

Every eigenvalue in the figures on this site comes from the cyclic Jacobi method, which is not the algorithm any library runs. LAPACK uses a tridiagonal reduction followed by QR iteration with shifts, which is several times faster.

Jacobi is used here for three reasons, in increasing order of importance. It is provably convergent — the sum of the squares of the off-diagonal entries decreases at every rotation, monotonically, with no shift strategy to get wrong. Each step is a plane rotation, so orthogonality is structural in the sense of a reflection cannot stop being one. And it computes the small eigenvalues to high relative accuracy, which the tridiagonal route does not, because reduction to tridiagonal form mixes the small eigenvalues with the large ones.

A site about accuracy should use the accurate method and say what it costs.

Three routes to the same spectrum #

The Jacobi figure carries three independent checks, and they are worth listing because together they are stronger than any one of them.

Every eigenpair satisfies its own equation. Av = λv to 10⁻¹², for each of them. This is the direct check and it is the weakest, because it would pass for a subset of the eigenvalues found twice.

The sum of the eigenvalues is the trace. The trace is invariant under every similarity transformation, so it is conserved by every rotation Jacobi performs, and the sum of the diagonal at the end must equal the sum of the diagonal at the start. This catches a lost or duplicated eigenvalue, which the first check does not.

The product of the eigenvalues is the determinant, computed by cofactor expansion. Another invariant, another arithmetic, and sensitive to a sign error in a way the trace is not.

Plus the eigenvector matrix’s orthogonality, ‖VᵀV − I‖, which for the six-by-six Hilbert matrix comes out at 10⁻¹⁵. A method that produced the right eigenvalues with a non-orthogonal V would pass everything above and be useless for anything downstream.

What this means in practice #

Three consequences, and the first is the one worth carrying furthest.

Exploit symmetry. A matrix that is symmetric, in code that does not know it, is being handed to a method with no guarantee on a problem that has one. Enforcing symmetry explicitly — replacing A with (A + Aᵀ)/2 when it should be symmetric and is not, because it was assembled by a process with rounding in it — costs nothing and buys the whole of Weyl. Costs nothing is doing more work in that sentence than it looks, and the next section is what it means.

What symmetrising actually costs #

The natural reading of “costs nothing” is that the eigenvalues move by about the asymmetry that was thrown away. They move by very much less than that, and the reason is exact rather than statistical.

Write A = S + K with S symmetric and K antisymmetric — which is what (A + Aᵀ)/2 and (A − Aᵀ)/2 are. The first-order shift of a simple eigenvalue of S under a perturbation K is vᵀKv for the real unit eigenvector v, and for antisymmetric K that quantity equals its own negative and is therefore zero. Every eigenvalue’s first-order term vanishes identically. The movement is second order.

Measured on a 6×6 symmetric matrix with a well-separated spectrum, with an antisymmetric part of Frobenius norm δ added:

A hundredfold in δ is a ten-thousandfold in the movement — δ² exactly, with the constant unchanged to three digits. At δ = 10⁻⁴ the discarded asymmetry moves the spectrum by 10⁻⁹, which is five orders below the thing discarded. So the advice is right and understated: symmetrising a matrix that is nearly symmetric is nearly free in a strong sense.

The exception is a near-degeneracy, and there the price is not accuracy but kind. Second-order perturbation theory divides by the gap between eigenvalues, so the cheapness above is bought against separation. Push two eigenvalues to within 10⁻³ of each other and apply δ = 10⁻², and the pair leaves the real axis: the eigenvalues of A have |Im λ| = 1.1·10⁻³, and the symmetrised matrix reports two real numbers.

That is not an error of 10⁻³ in a number. An oscillating mode has been reported as two decaying ones, and no statement about eigenvalue accuracy describes the difference — which is the same shape as the exact answer to a nearby problem, one level up: the computation is a correct answer to a question that has quietly become a different question. The rule that follows is narrow and checkable: symmetrise when the asymmetry is smaller than the smallest gap in the spectrum, and look before assuming it is.

For unsymmetric matrices, ask whether eigenvalues are needed at all. A great many uses of the spectrum — stability of a linear system, the behaviour of a matrix power, the response to a perturbation — are better answered by the singular values or by the pseudospectrum, both of which are well conditioned. Eigenvalues of a non-normal matrix can be a badly conditioned answer to a question that had a well-conditioned one.

No amount of precision helps. At $\varepsilon^{1/8}$, gaining sixteen more digits of precision buys two digits of eigenvalue accuracy. This is the clearest case on the site of a problem’s fault rather than an algorithm’s: the perturbation that does the damage is the one introduced by storing the matrix, before any algorithm has run.

What is asserted here #

The closed form is verified by substitution into the characteristic polynomial. Weyl’s bound is checked at all thirty-one perturbation sizes with the measurement floor written into the assertion. The Jordan curve at ε = 10⁻¹⁶ must exceed 10⁻³, which is the headline claim in checkable form. And Jacobi is required to refuse a matrix that is not symmetric rather than returning plausible numbers for it — a precondition, tested on every build by handing it an upper triangular matrix and requiring it to throw.

That last one matters more than it looks. A symmetric eigensolver applied to a nearly symmetric matrix returns numbers that are nearly right, and applied to an unsymmetric one returns numbers that are wrong in a way nothing else would flag. The precondition is the only thing standing between the site’s figures and a whole class of silent errors.

Pseudospectra, in one section #

If eigenvalues of a non-normal matrix move so much under perturbation, the natural question is what to use instead, and there is a good answer.

The ε-pseudospectrum is the set of complex numbers that are eigenvalues of some matrix within ε of A. For a symmetric matrix it is a small disc around each eigenvalue, of radius ε — that is Weyl’s inequality restated. For the Jordan block it is a disc of radius $\varepsilon^{1/n}$, which is enormous.

The point is that the pseudospectrum is the honest object. It answers “where could the eigenvalues be, given that my matrix is only known to within ε”, and for a non-normal matrix the answer is a large region that the eigenvalues themselves do not describe. Plotting the eigenvalues alone is plotting one point from a region and treating it as the region.

It also predicts behaviour the eigenvalues get wrong. A matrix whose eigenvalues are all inside the unit disc has powers that eventually decay, so the eigenvalues say the system is stable. For a non-normal matrix the powers can grow by many orders of magnitude first — transient growth, which the spectrum cannot see and the pseudospectrum can, and which is the reason fluid stability analysis based on eigenvalues alone gave wrong answers for decades.

Where symmetry comes from, and how it is lost #

Practically, most matrices that ought to be symmetric are, and the exceptions are worth recognising.

Symmetric by construction: Gram matrices AᵀA, covariance matrices, stiffness matrices from a variational formulation, graph Laplacians, Hessians of smooth functions. For all of these the symmetry is a theorem, and any asymmetry in the stored matrix is rounding from the assembly.

Symmetric in principle and not in storage: a matrix assembled entry by entry with floating-point arithmetic can be asymmetric in its last bits, and a symmetric eigensolver will refuse it or, worse, silently use only the upper triangle. Replacing A with (A + Aᵀ)/2 costs one pass and removes the question, and it is what careful code does.

Genuinely unsymmetric: transfer matrices, Markov chains, discretised advection, any system where the coupling from i to j differs from the coupling from j to i. These are the cases where the material in this essay is not a caution but the governing reality, and where the eigenvalues may not be the right question at all.

Two ways the eigenvalue problem is easier than it looks #

Balance against the pessimism, because for the symmetric case the news is unusually good.

Absolute accuracy is free. Any backward-stable symmetric eigensolver returns eigenvalues with absolute error at most about ‖A‖u. Not relative to the eigenvalue — relative to the norm of the matrix — which means the large eigenvalues come out to full precision regardless of how badly conditioned the matrix is.

The eigenvectors are only as good as the gaps. The sensitivity of an eigenvector is governed by its distance from the nearest other eigenvalue: a well-separated eigenvalue has a well-determined eigenvector, and two nearly equal eigenvalues have eigenvectors that can rotate freely in the plane they span. The subspace is still well determined; the individual vectors are not.

That second point is the eigenvector analogue of the valley with no bottom, and it has the same resolution: report the invariant subspace, which is determined, rather than the individual vectors, which are not.

What this field is about #

Three essays, and a single organising fact: the sensitivity of a spectral quantity has nothing to do with the condition number of the matrix.

This essay establishes that symmetric eigenvalues are perfectly conditioned in the absolute sense whatever κ is, and that non-symmetric ones can be arbitrarily worse. Rank is a decision establishes that rank is not a property but a judgement, and that the evidence for it degrades exactly in proportion to the noise. The best approximation there is establishes the one sharp equality in the subject and uses it as a check on the computation.

All three are about the singular value decomposition or its symmetric cousin, and all three come back to the same small quantity: the smallest singular value, or the smallest gap. Everything difficult in this field is that number being small.

The practical checklist #

Four questions, in the order worth asking them, when a spectrum is about to be computed.

Is the matrix symmetric? If it is by construction, enforce it in storage with (A + Aᵀ)/2 and use a symmetric solver. Everything gets easier and the guarantees are unconditional.

If not, are eigenvalues actually the question? Stability of a linear system, the behaviour of matrix powers, and the response to a perturbation are all better answered by singular values or by the pseudospectrum, and both are well conditioned for any matrix.

If eigenvalues are the question, how non-normal is the matrix? ‖AᵀA − AAᵀ‖ relative to ‖A‖² is a cheap measure, and a matrix that is nearly normal has nearly the symmetric case’s guarantees. A matrix that is far from normal does not, and the eigenvalues should be reported with the caveat.

How far apart are the eigenvalues? Eigenvectors are governed by the gaps rather than by symmetry, so even a symmetric matrix with two nearly equal eigenvalues has eigenvectors that are not individually determined — the invariant subspace is, and it is the thing to report.

None of those four costs more than a few lines, and between them they distinguish the case where the spectrum is a fact from the case where it is an artefact of the last bit of the input.