The eigenvalues that are not there

This site has measured eigenvalue sensitivity twice. A condition number for one eigenvalue is the reciprocal of the inner product of its left and right eigenvectors, and symmetry is worth more than precision found the 8×8 Jordan block whose eigenvalues move like ε^{1/8}, so a perturbation at the level of storing the matrix shifts them by a hundredth.

Both of those are statements about a limit. They describe how far an eigenvalue moves as the perturbation goes to zero, and the second one is careful to say that the exponent means the limit is reached slowly.

The pseudospectrum is the same question asked at a perturbation size that is actually present:

Λ_ε(A)  =  { z : σ_min(zI − A) ≤ ε }
        =  { z : z is an eigenvalue of A + E for some ‖E‖₂ ≤ ε }

Those two definitions describe the same set — the second follows from the first by taking E to be the rank-one matrix built from the singular vectors — and the second is the one that says why it matters. It is the set of places an eigenvalue can be put by a perturbation the size of the rounding in the entries.

The control, and it is an exact identity #

Set the superdiagonal to zero and the matrix is normal — AᵀA = AAᵀ, diagonalisable by an orthogonal matrix. For a normal matrix

σ_min(zI − A)  =  dist(z, Λ(A))

exactly, with no constant in front, because the singular values of zI − A are the moduli |z − λ_i|. So Λ_ε is the union of ε-discs around the eigenvalues, the picture carries nothing the spectrum did not, and there is no subject.

assertANormalMatrixHasNoPseudospectrum checks it at forty points scattered over the plane, by two routes that share no code: the smallest singular value of a 2n×2n real embedding, and the distance to the nearest eigenvalue computed from the site’s real Schur form. The worst relative departure over those forty points is 1.6·10⁻¹³.

That is the control the rest of the essay is measured as a departure from, and it is why the first stop of the slider is a normal matrix rather than a mild one.

What the shading is of, and how it is computed #

Each cell carries log₁₀ σ_min(zI − A) for z at its centre. Dark means small, which means near a place an eigenvalue could be moved to.

zI − A is complex and this site has no complex arithmetic. It does not need any: a complex matrix X + iY has the same singular values, each twice, as the real [[X, −Y], [Y, X]], so the existing real one-sided Jacobi SVD computes a complex singular value with four lines of index arithmetic and no second implementation to be wrong.

That is a deliberate choice rather than an economy. A second SVD is a second thing that can be wrong, and the site has spent two phases on findings that turned out to be defects in a measuring instrument rather than in the object.

The departure, measured along an axis #

At every stop past the first, the spectrum is identical: one eigenvalue at 0.8, with multiplicity six. No eigenvalue solver would report any difference between the five matrices the slider draws.

What changes is how far the shaded region reaches. assertANonNormalMatrixHasOne measures the ratio of the distance-to-spectrum to σ_min at scattered points, and finds it growing with the superdiagonal: 3.6·10⁴, 1.2·10⁶, 8.7·10⁶ and 3.7·10⁷ at m = 1, 2, 3, 4.

So the resolvent norm is not one over the distance; it is up to 3.7·10⁷ times that, and the factor is a monotone function of the one parameter the spectrum cannot see.

The same statement made by perturbing instead #

A resolvent norm is a computation about a matrix. The definition it comes from is about perturbations, and the two can be checked against each other.

The normal matrix’s eigenvalues move by 8.9·10⁻⁹ at the worst — the size of the perturbation, which is Bauer–Fike with an eigenvector condition number of exactly one. The bidiagonal’s move by 0.070, a factor of 7.9·10⁶ further on the same perturbations.

And the two computations agree, which is the check that neither is wrong: every perturbed eigenvalue λ̃ is asserted to satisfy σ_min(λ̃ I − A) ≤ ε. That is the definition rather than a finding, and it is exactly why it is worth asserting — if it failed, one of the two routes has a bug in it.

Why it moves like the n-th root #

The cloud grows as ε^{1/n} rather than as ε, which is why four decades of ε move it by about one decade rather than four.

The reason is the characteristic polynomial. For the bidiagonal with ρ on the diagonal and m above it, det(zI − A) = (z − ρ)ⁿ, and perturbing the bottom-left corner by δ changes it to (z − ρ)ⁿ − δ·mⁿ⁻¹. The roots are ρ + (δ mⁿ⁻¹)^{1/n} times the n-th roots of unity: a ring of n points at radius proportional to δ^{1/n}, which is what the figure draws.

At n = 6 and δ = 10⁻⁸ that radius is 10⁻⁸/⁶ ≈ 0.046 times m^{5/6}, and the measured 0.070 at m = 2 is exactly that. The site’s Jordan-block measurement is the same formula at n = 8.

A fractional power is what makes this a statement about real perturbations rather than about a limit. ε^{1/6} at ε = 10⁻¹⁶ is 2·10⁻³ — so even a perturbation at the level of double-precision rounding moves these eigenvalues by two parts in a thousand, and no amount of extra precision in the solver changes that, because the perturbation is in the stored matrix.

What a solver reports, and what it means #

An eigenvalue solver run on the bidiagonal returns 0.8 six times, to about 10⁻³. It is not making a mistake: it computed the exact eigenvalues of a matrix within rounding of the one it was given, which is what backward stability promises, and the eigenvalues of that nearby matrix genuinely are spread over a ring of radius 10⁻³.

So the output is correct, the algorithm is stable, and the six numbers reported are not six eigenvalues of anything in particular. The site’s identity accounts for it — forward error ⪅ condition number × backward error, with the condition number of a defective eigenvalue infinite and the ε^{1/n} law being what replaces the linear bound — and a condition number for one eigenvalue is where that number is measured.

Three ways to compute one, and why the picture is a grid #

Computing a pseudospectrum is computing σ_min(zI − A) at a lot of points, and how many points is the whole cost.

The direct way, which is what this figure does: one SVD per grid point. At 34×34 points and a 6×6 matrix that is 1,156 SVDs of a 12×12 real embedding, which is a fraction of a second. At a 1,000×1,000 matrix and a 200×200 grid it is forty thousand SVDs of a 2,000×2,000, which is not a computation anybody performs.

The Schur way, which is what any serious code does: reduce A to triangular form once, and note that σ_min(zI − A) = σ_min(zI − T) because a unitary similarity does not change singular values. Then each grid point is a triangular matrix, and its smallest singular value comes from a few steps of inverse iteration at O(n²) rather than O(n³). That is the algorithm behind EigTool and behind pseudospectra in every package that has one.

And the projected way, for matrices too large even for that: run a few dozen steps of Arnoldi, compute the pseudospectrum of the small Hessenberg matrix, and hope. It is a lower bound on the real one in the sense that the projected resolvent norm is smaller, so the region drawn is inside the true one — one-sided again, and again in the direction that under-reports the trouble.

The choice made here is the direct one, at a small size, because the argument is about what the object is and not about how to compute it at scale. The cost is stated at the call site rather than hidden, which is why nx and ny are arguments the figure’s assertions bound.

What the contours are not #

Two readings of this picture are wrong and both are natural.

It is not a plot of where the eigenvalues “really” are. The eigenvalues are at 0.8, exactly, six times. The region says where they could be moved to, which is a statement about a family of matrices rather than about this one.

And the shaded bands are not error bars. An error bar is symmetric and describes a distribution; this region is a set, with a hard boundary, and every point inside it is achievable by some perturbation of norm ε. Whether a particular perturbation reaches a particular point is a different question, and the cloud in the perturbation figure is what its answer looks like: the achievable set is a disc and the achieved set is a ring, because a random perturbation almost never picks the rank-one direction that reaches the interior.

That gap between what is achievable and what is achieved is the same distinction the condition number is an amplifier measured for condition numbers, where a random perturbation direction gets 0.29 of the available amplification at the median.

Where this changes an answer #

Three places, and none of them is about eigenvalues being printed.

Convergence of an iteration. x_{k+1} = Ax_k decays if ρ(A) < 1 and the spectrum says nothing about what it does first. That is the next essay, a spectral radius that grows first.

Krylov methods. The spectrum that predicts nothing already measured GMRES on a matrix whose eigenvalues are all on the unit circle and which makes no progress for n − 1 steps. The pseudospectrum is the object that does predict it: GMRES’s residual is bounded by the length of the ε-pseudospectrum’s boundary divided by ε, and for a normal matrix that reduces to the eigenvalue bound everyone quotes.

Stability of a discretisation. A time-stepping scheme is stable if the amplification matrix has spectral radius below one, in the sense that solutions eventually decay. Whether they decay monotonically — which is what a physical argument usually wants — is a question about the pseudospectrum, and the Kreiss matrix theorem is the statement that connects them.

What a backward-stable eigenvalue solver promises #

The sentence a reader carries away from every stability discussion is that a good algorithm returns the exact answer to a nearby problem, and it is worth watching that sentence do its work here, because it does exactly what it says and the result is unhelpful.

An eigenvalue solver run on the bidiagonal returns six numbers near 0.8. What backward stability promises is that those numbers are the exact eigenvalues of A + E for some E with ‖E‖ ≤ c(n)·u·‖A‖. That promise is kept.

But this essay’s whole content is that Λ_ε(A) for ε at that size is a ring of radius 10⁻³. So the set of matrices the promise covers has eigenvalues spread over that ring, and any six numbers in it are a correct output. The algorithm is not choosing badly among them; every choice is equally correct.

Two consequences follow, and the second is the one worth carrying.

The output is not reproducible in a strong sense. Two backward-stable solvers, or one solver on two machines with different summation orders, can return numbers differing in the third digit and both be right. That is not a bug report; it is what the promise permits.

And no amount of precision fixes it. Doubling the working precision shrinks ε by 10⁻¹⁶, and the ring’s radius by the sixth root of that — a factor of 400. Which is real, and is not the sixteen digits a reader would expect from sixteen more digits of arithmetic.

Why the definition uses the 2-norm #

A detail that changes the picture and is usually left implicit.

Λ_ε(A) = { z : ‖(zI − A)⁻¹‖ ≥ 1/ε } needs a norm, and the two definitions given at the top of this essay coincide only in the 2-norm. The reason is that the perturbation achieving the bound is the rank-one matrix built from the singular vectors of zI − A, and its 2-norm is exactly σ_min while its Frobenius norm and its 1-norm are not.

In another norm the set is still defined and still meaningful, and the equivalence between “the resolvent is large here” and “a perturbation of this size puts an eigenvalue here” acquires a dimensional constant. Every picture in this essay is in the 2-norm for that reason, which is also why the computation is a smallest singular value rather than anything cheaper.

What is worth carrying #

For a normal matrix the pseudospectrum is a set of discs and carries nothing. That is not a caveat; it is why the object exists, and it is asserted here as an exact identity so everything else is a measured departure from it.

For a non-normal one, a perturbation the size of rounding moves the eigenvalues by a fractional power of it — the n-th root for a defective eigenvalue of multiplicity n — so the movement is enormous compared to the perturbation and gets relatively larger as the perturbation shrinks.

Two matrices with identical spectra can behave completely differently, and the spectrum is the thing every library reports.

And the eigenvalues a solver returns for such a matrix are correct and not informative. The algorithm is backward stable, the numbers are the exact eigenvalues of a nearby matrix, and nearby is doing an enormous amount of work in that sentence.