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 $\varepsilon^{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:

$$\begin{aligned} \Lambda_\varepsilon(A) &= \{\, z : \sigma_{\min}(zI - A) \le \varepsilon \,\} \\ &= \{\, z : z \text{ is an eigenvalue of } A + E \text{ for some } \|E\|_2 \le \varepsilon \,\} \end{aligned}$$

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

$$\sigma_{\min}(zI - A) = \operatorname{dist}\bigl(z, \Lambda(A)\bigr)$$

exactly, with no constant in front, because the singular values of zI − A are the moduli |z − λᵢ|. 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 — the case where an eigenvalue’s condition number is one — 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₁₀ σₘᵢₙ(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.

One more unit of coupling, and the reason for drawing all five rather than the ends is that the first step and the later ones are not the same size:

The spectral radius is 0.8000 at all five stops and the pseudospectral boundary runs from 0 to 1.687. Five matrices with identical spectra — the same eigenvalues, the same radius, to four decimal places — and a resolvent that reaches half again outside the unit disc at one end of the slider and not at all at the other. Whatever the eigenvalues are measuring, it is not this.

And the crossing happens at the first step. From no superdiagonal to one, the 10⁻³ level goes from 0.000 to 1.116 — straight past 1, which is the boundary that decides whether a power iteration grows. The remaining three units of superdiagonal add 0.215, 0.210 and 0.146: the effect is not gradual in the non-normality, it is present as soon as there is any, and then it accumulates slowly.

So the practical reading is harsher than “the spectrum is not the whole story”. A matrix with a spectral radius of 0.8 is one a stability argument would pass without hesitation, and a single unit of coupling above the diagonal is enough to put its 10⁻³ pseudospectrum outside the unit circle — at which case a perturbation of relative size 10⁻³, which is nothing, produces a matrix that genuinely does grow. The spectral radius has not become approximate here. It has become the answer to a question nobody asked.

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 σₘᵢₙ(λ̃ 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 $\varepsilon^{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 $\rho + (\delta m^{n-1})^{1/n}$ times the n-th roots of unity: a ring of n points at radius proportional to $\delta^{1/n}$, which is what the figure draws.

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

Two experiments, and the fifteen per cent between them #

Those two numbers are close and they are not the same measurement, which is worth separating because the reason they are so close is itself the point.

The formula is exact for the perturbation it describes. Perturbing the corner alone by δ, over three sizes, three superdiagonals and two perturbation sizes, the ratio of the measured ring radius to (δ·m^(n−1))^(1/n) is 1.0000 in all eighteen cases. There is nothing approximate about the derivation; it is a factorisation of a characteristic polynomial and it holds to four digits.

The figure applies a different perturbation. A random full matrix of the same 2-norm, which is what a rounding looks like, moves the ring by less — because only part of it lands in the corner. As a fraction of the corner formula’s radius:

            worst of 24 draws     median of 24
   n = 4         0.880                0.659
   n = 6         0.857                0.700
   n = 8         0.916                0.814

Constant to three digits across every perturbation size and every superdiagonal, and a weak function of n. So the 0.070 is 0.857 times the 0.0827, and the gap is this ratio rather than a rounding or a sloppy measurement.

Why it is only fifteen per cent is the interesting half. The share of a random unit-norm perturbation that reaches the one entry the formula is about is roughly 1/n — a factor of six at n = 6 — and the radius is that share raised to the power 1/n, which turns a factor of six into a factor of 1.35. The fractional power compresses the modelling error along with everything else.

That is worth carrying past this essay, and past the spectrum that predicts nothing. An ε¹ᐟⁿ law is so flat that the difference between a carefully chosen worst-case perturbation and an arbitrary one is a few tens of per cent — which is what makes the law useful as a prediction about real roundings rather than only about adversarial ones, and it is the same flatness that makes four decades of ε move the ring by one.

It also says something about how a claim of this shape should be checked. The corner perturbation is the one the algebra can handle and the random one is the one a reader has; a derivation that used the first and a figure that draws the second are answering the same question only if the constant between them is small, and nothing about the derivation says that it is. Here it is 0.86 and the essay was right to treat the two as one measurement. On a law with an integer power rather than a fractional one — an ordinary condition number, say, where the exponent is 1 — the same substitution would cost a factor of n rather than n¹ᐟⁿ, and the two experiments would be six apart at n = 6 rather than fifteen per cent.

So the licence to conflate them is not general; it is bought by the fractional power, which is the very thing the essay is about. That is a pleasant kind of consistency and it is also the reason the gap is worth printing rather than smoothing: a reader who carries “the formula predicts the measurement” to a problem with a linear sensitivity will find it does not.

A fractional power is what makes this a statement about real perturbations rather than about a limit. $\varepsilon^{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 $\varepsilon^{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 σₘᵢₙ(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 σₘᵢₙ(zI − A) = σₘᵢₙ(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ₖ₊₁ = Axₖ 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.

An equation where the same distinction decides the answer #

A non-normal matrix’s spectrum is an inadequate description of it, and this page draws that directly. The consequence for a matrix equation is that its conditioning is not the eigenvalue gap, and the gap is the number a reader is invited to consult.