Eigenvalues, singular values, rank

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.

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 ε^(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.

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. 1 How far the eigenvalues move, against the size of the perturbation, for two matrices. The symmetric case follows a line of slope one — an eigenvalue moves by at most the size of the perturbation, which is Weyl’s inequality and has no constant in it. The Jordan block follows a line of slope one eighth.

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 ε^(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 |λ| = ε^(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.

Jacobi sweeps on the 6×6 Hilbert matrixA plot of the off-diagonal norm against sweep number, falling from about one to ten to the minus seventeen in five sweeps, with the matrix shown at three stages beneath it.01234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweep‖off-diagonal‖eigenvalues, largest first1.6190.24240.016326.157·10⁻⁴1.257·10⁻⁵1.083·10⁻⁷κ = 1.5·10⁷‖VᵀV − I‖ = 1.9·10⁻¹⁵1.000.500.330.250.200.170.500.330.250.200.170.140.330.250.200.170.140.130.250.200.170.140.130.110.200.170.140.130.110.100.170.140.130.110.100.09as given1.620.050.00-0.040.00-0.000.050.240.000.010.00-0.000.000.000.000.00-0.00-0.00-0.040.010.000.020.000.000.000.00-0.000.000.00·-0.00-0.00-0.000.00·0.00after one sweep1.62······0.24······0.00······0.02······0.00······0.00after 8 sweepseach sweep is n(n−1)/2 rotationsthe trace is conserved
Fig. 2 Jacobi on the six-by-six Hilbert matrix. The off-diagonal norm falls from about one to below 10⁻¹⁷ in four sweeps, and the fall is quadratic once the rotations are small. Beside the plot: the six eigenvalues, spanning seven orders of magnitude, and the orthogonality of the eigenvector matrix. The trace is conserved by every rotation and is checked to 10⁻¹³ against the trace of the original.

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.

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 ε^(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.

How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 3 The comparison worth making. For a linear system, the amplification is κ and it is a fixed factor — the distribution here slides one decade per decade and keeps its shape. For a defective eigenvalue problem there is no such factor: the relationship is a fractional power, and no single number describes it.

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 ε^(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.

Singular values of a rank-4 matrix with noise of relative size 10⁻⁸Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 4 A question that is well conditioned for any matrix, symmetric or not. Singular values are as stable as symmetric eigenvalues — a perturbation moves them by at most its own size — which is why they are the safe tool when the eigenvalues are not.

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.

Error of the best rank-k approximation to a 10×10 matrixApproximation error against k on a logarithmic axis, with the measured error and the next singular value drawn as separate curves lying exactly on top of one another.12345678910⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 4.3·10⁻⁹worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁4.3·10⁻⁹worst Frobenius discrepancy4.3·10⁻⁹κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁
Fig. 5 The related well-conditioned question. The best rank-k approximation’s error is exactly σₖ₊₁ — an equality rather than a bound, and a quantity as stable as the singular values are. Where an eigenvalue question is ill-posed, the corresponding singular value question is usually 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.

Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 6 The framework this field departs from. For a linear system the algorithm’s contribution is flat and the problem’s rises with κ. For a defective eigenvalue problem there is no κ to plot against: the relationship is a fractional power, and the perturbation that does the damage arrived before any algorithm ran. See the exact answer to a nearby problem for the vocabulary and the valley with no bottom for the least-squares analogue.

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.