The gap decides the eigenvector

This site has measured eigenvalue sensitivity twice. Symmetry is worth more than precision establishes the good case: a symmetric matrix’s eigenvalues move by at most the size of the perturbation, however clustered they are, while an unsymmetric one’s can move by the eighth root of it. A condition number for one eigenvalue gives each eigenvalue its own sensitivity, 1/|yᵀx| — the reciprocal of the cosine between its left and right eigenvectors — and checks it against a closed form.

Both of those are about eigenvalues. The other half of an eigenpair is governed by a completely different quantity, and the two are routinely conflated:

An eigenvalue’s sensitivity is set by the angle between its own two vectors. An eigenvector’s sensitivity is set by the gap to its neighbours.

A symmetric matrix has no angle to speak of. Its left and right eigenvectors are the same vector, so every eigenvalue’s condition number is exactly one, and on that axis the problem is as good as a problem can be. Its eigenvectors can still be arbitrarily badly conditioned, and this essay measures how badly, on matrices whose eigenvalues are as well behaved as the previous essays promise.

The experiment #

A symmetric matrix with a spectrum chosen so that exactly one thing varies. Two eigenvalues at 1 and 1 + gap; the rest at 4, 5, 6, 7, 8, 9. The matrix is VΛVᵀ for a fixed random orthogonal V, so the eigenvectors are the same in every run and only the gap changes.

The perturbation is symmetric — otherwise the comparison would be between a symmetric problem and an unsymmetric one, which puts the eigenvalue condition numbers back into play and is a different essay’s subject. It is normalised to a chosen Frobenius norm and generated from a fixed seed, so the same E is added at every gap.

Then three quantities are measured on the same run: how far the eigenvalues moved, how far the two eigenvectors turned, and how far the plane they span turned. Measuring all three from one perturbation is what makes the comparison a comparison rather than three anecdotes.

The eigenvalues do not notice #

Across six decades of gap at a perturbation of 10⁻⁶, the largest eigenvalue shift is 3.867·10⁻⁷ — every time. The spread of that number over the whole sweep is a factor of 1.0000001.

That is Weyl’s inequality, which says the k-th eigenvalue of A + E differs from the k-th of A by at most ‖E‖₂, and it holds with no reference to the spectrum’s shape. Two eigenvalues a billionth apart are each perturbed by no more than the perturbation, exactly as two eigenvalues a hundred apart would be.

So a reader who asked “is this problem well conditioned?” and was handed the eigenvalue sensitivities would be told: perfectly. Every condition number is one. Nothing in the answer would mention what happens next.

The eigenvectors notice a great deal #

The same runs, the same perturbation, the two eigenvectors of the cluster:

The first five rows are a straight line of slope −1 on a logarithmic plot: the angle is proportional to ‖E‖/gap, which is the sinθ theorem of Davis and Kahan. The bound is drawn on the figure and the measurement sits an order of magnitude below it, in the way this site’s other bounds do — a bound is a bound and the slack is worth printing.

At 10⁻⁹ the angle is 0.478 radians, which is 27°. The computed eigenvector and the true one have a cosine of 0.888. Nothing has failed, no algorithm has been sloppy, and every digit of the arithmetic is intact.

Where the gap comes from #

The bound is not an empirical fit and the derivation is three lines, so it is worth having rather than citing.

Perturb A by E and expand the eigenvector of λᵢ to first order in the eigenbasis of A. Writing the perturbed vector as xᵢ + δ and requiring (A + E)(xᵢ + δ) = (λᵢ + δλ)(xᵢ + δ), the terms of first order give, for each j ≠ i,

xⱼᵀ δ  =  (xⱼᵀ E xᵢ) / (λᵢ − λⱼ)

so the perturbation of the eigenvector is a sum over the other eigenvectors, each weighted by how much E couples it to xᵢ divided by the distance between their eigenvalues. The eigenvalue’s own first-order shift, by contrast, is xᵢᵀExᵢ — no denominator at all, which is Weyl’s inequality’s first-order shadow and the reason the eigenvalues in the table above do not move.

Every term in that sum is bounded by ‖E‖ over its own gap, so the whole displacement is bounded by ‖E‖ over the smallest gap, and the smallest gap is the one to the neighbour. That is the sinθ theorem in its crudest form; Davis and Kahan’s contribution is that the same bound holds without the first-order assumption, which matters because the interesting regime here is precisely the one where the expansion stops converging.

The expression also says which perturbations matter. Only the component of E that couples xᵢ to xⱼ appears, so a perturbation orthogonal to that coupling moves nothing however large it is — which is why the four seeds at a fixed gap give four different answers rather than four noisy copies of one.

And it says where the denominators are small. A cluster of m eigenvalues has m(m−1) small gaps between its members and m(n−m) large ones to everything outside, so the sum splits into a part with tiny denominators — confined to the cluster’s own plane — and a part with ordinary ones. That split is the whole of the next essay, visible here in the derivation before it is measured.

Saturation, and what it means #

Below a gap of about 10⁻⁷ the angle stops climbing. Measured at 10⁻⁹, 10⁻¹⁰ and 10⁻¹¹ it is 0.478, 0.483 and 0.484 — a spread of 1.3% across three decades of gap, where the previous decade had cost a factor of ten each.

The bound has not been violated; ‖E‖/gap is 1,000 there and the angle is a radian, so the bound has simply stopped being informative. What has happened is more interesting than saturation at a limit.

Once the gap is well below the perturbation, the perturbed matrix’s behaviour inside that plane is dominated by the perturbation’s own 2×2 block, not by A’s. The eigenvectors that come back are essentially the eigenvectors of E restricted to the plane. So they have stopped being a property of A at all — and the test for that is to change E rather than the gap.

At a gap of 10⁻¹⁰, four different perturbations of the same size 10⁻⁶ return angles of 0.483, 1.466, 0.536 and 0.167 radians. Nothing about the matrix changed between those four runs. The answer is decided entirely by which perturbation was applied, which is to say by which rounding errors the arithmetic happened to make.

That is the sentence worth carrying away. Below the saturation point the computation has not degraded — it has changed subject.

A double eigenvalue has no eigenvector #

Push the gap to zero and the conditioning language stops applying, because there is no quantity left to be badly conditioned.

If λ is a double eigenvalue of a symmetric A, then every vector in the corresponding plane satisfies Ax = λx. There is no correct eigenvector, so there is nothing for an algorithm to fail to return. Three perturbations of size 10⁻¹⁴ — fourteen orders of magnitude below anything a tolerance would notice — return first eigenvectors at 0.034, 1.056 and 0.318 radians from each other, and each of them satisfies Ax = λx with a residual of 5·10⁻¹⁵.

None of them is wrong. That is the point of measuring the residual: it refuses the reading that the disagreement is error. Three algorithms, or one algorithm on three machines, or the same code compiled twice, would produce three different answers and all three would be right.

No residual test can see it #

The standard defence against a wrong answer on this site is the backward error: compute the residual and, if it is small, the answer is the exact answer to a nearby problem. That defence is exactly useless here, and the reason is the same one a small residual is not a small error is about.

Take the computed eigenvector at a gap of 10⁻⁹, the one 27° from the true one. Its residual ‖Ax̂ − λ̂x̂‖ is at rounding — around 10⁻¹⁵ — because x̂ is very nearly an eigenvector of a matrix very near A, and being 27° from the true vector costs nothing in that measurement. Every vector in the plane has a residual of about 10⁻¹⁵, and the residual therefore cannot distinguish between them.

The general statement, which this subject has in a sharper form than most: for a symmetric matrix, a small residual bounds the distance from λ̂ to some eigenvalue, and it bounds the angle to the corresponding eigenvector only after division by the gap. So the residual is a complete guarantee about the eigenvalue and a guarantee about the eigenvector that is worth exactly as much as the gap is large.

An eigensolver’s convergence test is a residual test. It is doing its job — the eigenvalue really is accurate — and it has nothing to say about the quantity a caller is likely to use. The gap is computable, cheaply, from the eigenvalues the same routine already produced, so the information is available and is not usually reported.

The instrument has a floor, and it looks like a result #

The plane’s angle is measured by principal angles, and computing them the way they are defined does not work. This is worth its own section because the defect produces a plausible figure.

The definition: the principal angles between two subspaces are the arccosines of the singular values of Q₁ᵀQ₂ for orthonormal bases Q₁ and Q₂. For a small angle those singular values are near one, where the arccosine is infinitely steep — acos(1 − δ) ≈ √(2δ) — so an error of 10⁻¹⁶ in the singular value comes out as an error of 1.5·10⁻⁸ in the angle. Half the digits, and a floor no small angle can be measured below.

That is the same catastrophic cancellation the arithmetic field draws, arriving inside a measuring instrument. A curve of subspace angle against gap, computed the obvious way, falls and then goes flat at 1.5·10⁻⁸ — and the flat part reads as the subspace angle having stopped shrinking.

The repair is standard and is not a repair everywhere: compute the part of Q₂ that Q₁ cannot reach, (I − Q₁Q₁ᵀ)Q₂, and take the arcsine of its singular values. Near zero the arcsine is flat and well conditioned, and the route follows the true angle down to 10⁻¹⁵. It is then wrong at the other end, because a right angle has sine one and the arcsine is as steep there as the arccosine was at zero: asked for the angle between two planes sharing nothing, it returns 1.5707962934 against a true 1.5707963268 — the same eight digits, lost at the opposite end of the range.

So neither formulation is right everywhere and the code uses both, taking each angle from whichever inverse function is flat for it. The floor is measured against √(2u) — 1.4901161193847656·10⁻⁸, which is the value the arccosine route bottoms out at, to every digit.

What a gap is measured against #

One loose end, and it is the sort that quietly makes a claim meaningless: a gap of 10⁻⁹ is small compared with what?

The bound ‖E‖/gap is a ratio of two quantities carrying the same units, so it is scale-free in the right way — double the matrix and both double. But the interesting comparison is between the gap and the perturbation, and the perturbation in a real computation is not chosen. It is the rounding of the arithmetic, which for a symmetric eigensolver is about u‖A‖ — the unit roundoff times the norm of the matrix.

So the practical threshold is a gap of order u‖A‖, below which the computed eigenvectors are decided by the arithmetic’s own error rather than by the matrix. On the matrices in this essay ‖A‖ is 9 and u is 1.1·10⁻¹⁶, so that threshold sits near 10⁻¹⁵ — far below the gaps drawn here, which is why every perturbation in this essay had to be added on purpose.

That is worth being clear about. A double-precision eigensolver on a well-scaled matrix has a great deal of room, and the gaps at which this essay’s effect bites in practice are the ones a problem supplies: a symmetry that makes two eigenvalues exactly equal, a near-degeneracy from a physical near-symmetry, or a cluster that tightens as some parameter is varied. Those are common, and none of them is a rounding phenomenon.

The other route to the same regime is the one the mixed-precision essays opened. At bfloat16’s eight significand bits u is 2⁻⁸, so u‖A‖ is about 0.035 on the same matrix, and a gap of a hundredth is then below the threshold. Every eigenvector of a cluster tighter than that is decided by the arithmetic. The gap that matters is not a property of the matrix alone; it is a property of the matrix and the precision together, which is the pairing this whole site is organised around.

What is asserted, and what refuses #

No eigenvalue moves by more than the perturbation, at every gap, which is Weyl checked rather than quoted.

Every eigenvector angle is inside its Davis–Kahan bound, ‖E‖/gap, at every gap on the sweep.

The eigenvalue shift does not move with the gap and the eigenvector angle moves by more than three orders of magnitude — asserted as a pair, so that the contrast is one check.

Saturation is measured in both halves: the angle stops depending on the gap below the perturbation, to within 5% across three decades, and starts depending on the perturbation instead, by more than a factor of two across four seeds.

Every returned eigenvector satisfies Ax = λx to 5·10⁻¹⁵ at the double eigenvalue, which is what refuses the reading that any of them is an error.

And the claim symmetry protects the eigenvectors is refused. “A symmetric matrix has well-conditioned eigenvectors” is fed to rejects() at a gap of 10⁻⁸ and reports an angle of a radian and a half.

What is left #

The plane, which is the next essay and is the practical answer: the invariant subspace the two vectors span is governed by the gap to the rest of the spectrum, which never got small, and it moves by 7.6·10⁻⁸ where the vectors move by half a radian.

The unsymmetric case, where the eigenvector’s sensitivity involves both the gap and the eigenvalue’s own condition number, and the two effects compound. Everything measured here is symmetric, deliberately, so that only one quantity is varying.

Inverse iteration and the algorithms that actually compute these vectors, none of which appears here. What is measured is the conditioning of the problem, which is the ceiling every algorithm runs into, and the algorithms deserve their own measurements rather than an inference from this one.