The vector moves with its eigenvalue

The middle eigenvectors come back exact put the Sylvester–Kac matrix’s eigenvectors against the exact ones, which are Krawtchouk polynomials evaluated at the integers and can be computed in rational arithmetic to any order. It found that a feared loss did not happen. The matrix’s middle eigenvalues lose up to eleven digits at order 80 — they are the most ill-conditioned — and their eigenvectors, computed by three steps of inverse iteration with the computed eigenvalue as the shift, from a random start, are exact to fourteen digits at every order. An earlier result that looked like a loss had been a start vector with no component along half the eigenvectors.

Its closing section named the route on which the loss should happen after all: “Back substitution in the computed quasi-triangular factor and multiplication by the accumulated orthogonal transformation is how LAPACK computes eigenvectors, and it inherits the similarity’s backward error. The prediction with a sign is that on that route the middle eigenvectors’ normwise error grows like the middle eigenvalues’ error divided by the gap — about $10^{-5}$ at order 80 — so that the earlier essay’s prediction holds there and fails here, and the difference between the two routes is the whole of the answer.”

It holds, in its shape and its size. Its constant is a different number with a different reason, and the reason is the more interesting half.

Two routes from one Schur form #

Both routes start from the same computation. The Francis double-shift iteration reduces the Sylvester–Kac matrix to real Schur form, $A = ZTZ^{\mathsf T}$ with $Z$ orthogonal and $T$ upper triangular — every eigenvalue of this matrix is real, and up to order 80 the computed form has no two-by-two blocks. The diagonal of $T$ is the computed spectrum, the numbers an answer that is known measured against the integers $-(n-1), -(n-3), \ldots, n-1$.

The first route is the earlier essay’s: take a computed eigenvalue $\hat\lambda$, factor $A - \hat\lambda I$ once, and solve with it three times from a random vector. The second is LAPACK’s. An eigenvector of the triangular $T$ for its $j$-th diagonal entry has zeros below position $j$, a one at $j$, and entries above found by back substitution, each a division by $t_{ii} - t_{jj}$; the eigenvector of $A$ is $Z$ times that. It costs no factorisation and returns every eigenvector at once, which is why it is the standard route. The form a real matrix can reach is the account of why the Schur form is the object every dense eigensolver aims at, and two shifts that are never formed of how the Francis iteration reaches it without complex arithmetic; neither route here touches the iteration itself, only what is done with its output.

The figure is the prediction’s verdict. By inverse iteration the worst middle eigenvector is off by $2 \cdot 10^{-16}$ to $5 \cdot 10^{-14}$ at every order. From the Schur form it is off by $10^{-14}$ at order 16, $2 \cdot 10^{-9}$ at 48, $8 \cdot 10^{-7}$ at 64 and $2.3 \cdot 10^{-5}$ at 80 — the prediction’s $10^{-5}$, within a factor of three — and it climbs in step with the middle eigenvalues’ error, never more than a little under it. The two routes start from the same computed eigenvalues and end nine orders of magnitude apart.

Every eigenvector, just under its eigenvalue #

The middle is where the loss is largest, but not where it happens.

At order 80 the eigenvalue errors make an arch over the spectrum, from $10^{-12}$ at the ends to $3 \cdot 10^{-5}$ in the middle, and the Schur route’s eigenvector errors make the same arch a fixed step below it, at every eigenvalue where the arch is above $10^{-11}$. Inverse iteration’s errors are flat at $10^{-15}$ across the whole spectrum, rising to $10^{-13}$ only in the middle. Turn the dial down and the arch flattens, at order 48 to a peak of $3 \cdot 10^{-9}$, and the red points follow it down; the blue ones do not move.

The two dips in the arch near $\pm 24$ are eigenvalues the Francis iteration happened to converge to more accurately; their eigenvectors from the Schur form dip with them. Nothing about the eigenvector’s own sensitivity changes there. Whatever the Schur route does to an eigenvector, it does in proportion to what happened to that eigenvalue.

A ratio that is not one over the gap #

The prediction’s constant was one over the gap: an eigenvector’s sensitivity to a perturbation is the perturbation over the distance to the neighbouring eigenvalues, and here every distance is two, so the error should be half the eigenvalue’s. That is the reasoning the gap decides the eigenvector made precise for symmetric matrices.

On 274 such eigenvectors the ratio is between 0.76 and 0.84, and nearly all of them are within a hundredth of 0.79. It does not depend on the order, on where the eigenvalue sits in the spectrum, or on how large its error is across six decades from $10^{-11}$ to $10^{-5}$. It is not a half.

The 0.79 can be found without any eigensolver. The matrix’s own recurrence — the first equation of $(A - \lambda I)x = 0$ gives $x_1$ from $x_0$, the second $x_2$ from $x_1$ and $x_0$, and so on — can be run at any $\lambda$, not only an integer, exactly, in rational arithmetic on the double’s binary value. At an eigenvalue it gives the eigenvector; at a nearby $\lambda$ it gives a vector that satisfies every equation but the last. Differencing the unit vectors at $\lambda \pm 2^{-20}$ measures how fast the eigenvector turns as its eigenvalue moves.

It turns at 0.787 per unit in the interior of the spectrum at every order from 48 up, rising to 0.91 at the two ends. On every eigenvector of orders 48 and 80 whose error is measurable the speed matches the Schur route’s ratio to within three per cent. And measured directly, on nine eigenvectors checked at orders 48, 64 and 80, the Schur route’s vector sits exactly as far from the exact one as the recurrence’s vector at the computed eigenvalue does, to three digits — on the opposite side of it, to within a fifth of that distance.

So the Schur route’s error has a plain description. It computes the eigenvector of the eigenvalue it was given, and the eigenvalue it was given is wrong. Back substitution in $T$ for the diagonal entry $t_{jj}$ is an exact computation, to rounding, of a vector that the triangular matrix maps to $t_{jj}$ times itself; multiplying by $Z$ carries it to the matrix that $Z$ and $T$ represent, which is $A$ to within rounding. What is wrong is $t_{jj}$, and the vector moves with it at the speed the matrix’s eigenvectors move — 0.79, for this matrix, not the 0.5 a gap would suggest, because the gap governs how a perturbation of the matrix mixes in the neighbouring eigenvectors, and this error is not a mixing. It is a slide along the curve the eigenvector traces as its eigenvalue changes.

Why inverse iteration does not slide #

Inverse iteration uses the same wrong eigenvalue and returns the right vector, and the difference is what the computation does with the eigenvalue. In back substitution $\hat\lambda$ is part of the answer: the vector is defined as the one that $\hat\lambda$ makes singular. In inverse iteration it is only a shift. A solve with $A - \hat\lambda I$ multiplies the component along each eigenvector by one over its distance from the shift, so the wanted component is amplified by $1/|\lambda - \hat\lambda|$ and the others by about a half; after three solves the wanted direction dominates by a factor that is enormous whether $\hat\lambda$ is right to fifteen digits or to five. The vector it converges to is $A$’s eigenvector, not $\hat\lambda$’s, because the iteration never asks which eigenvalue it belongs to. At order 80 the worst middle eigenvalue is off by $3 \cdot 10^{-5}$, so each solve amplifies the wanted direction by about thirty thousand against a half for its nearest neighbours: a ratio of sixty thousand per solve, $2 \cdot 10^{14}$ after three, which is why three solves from a random start reach rounding at every order measured.

That is the earlier essay’s result seen from the other side. Its explanation was that a shifted solve is insensitive to the shift’s error; this one adds that the route which uses the shift as a definition is exactly as sensitive as the eigenvector’s speed allows. The practical repair follows directly. Take the Schur route’s vector as the start of inverse iteration, shifted by the same computed eigenvalue, and solve once. At orders 48 and 64 every eigenvector is then within $8 \cdot 10^{-14}$ normwise; at 80, within $1.9 \cdot 10^{-10}$, from a slide of up to $4 \cdot 10^{-5}$, and a second solve would finish the job. One solve multiplies the slide by about the eigenvalue’s error over the gap, which is the shifted solve’s insensitivity working on a start that is already close. LAPACK offers the inverse-iteration route as an alternative to back substitution, working on the Hessenberg form; on this matrix the choice between them is the difference between five digits and fourteen.

A condition number for one eigenvalue is the account of why the middle eigenvalues are wrong by so much — the non-normality the left and right eigenvectors carry — and the Schur route passes that error on to the vectors at a fixed rate while inverse iteration does not pass it on at all.

The ends, lost twice #

Normwise, the Schur route’s vectors at the ends of the spectrum are as good as inverse iteration’s. Componentwise they are not.

An eigenvector near the ends of the spectrum has entries spanning many decades, and the earlier essay found that inverse iteration’s componentwise error is the rounding unit times that span, within a few units, at every order. The Schur route’s is a hundred to two thousand times past it at orders 48, 64 and 80. The vector is formed as $Z$ times the triangular eigenvector, and each entry of a matrix-vector product carries an error relative to the largest terms in its sum, not to its own size; the small entries at the far end of an end eigenvector are sums of large terms that cancel. Inverse iteration computes each entry as part of a solve whose residual is small entry by entry. Both are backward stable, normwise; only one of them preserves the vector’s small entries.

What the prediction was right about #

The prediction’s claim was that the route is the whole of the difference, and it is: the same eigenvalues, the same matrix, nine orders of magnitude between the routes in the middle at order 80. Its size was right within a factor of three. What it had wrong was the mechanism behind the constant, and the difference matters for any matrix other than this one. One over the gap is the right constant for an eigenvector perturbed by a perturbation of the matrix. The Schur route’s error on this matrix is dominated by something more specific — the triangular form’s diagonal is wrong by the eigenvalue’s error, and the eigenvector is computed from the diagonal — and for that, the constant is the speed at which the eigenvector moves along its own family as the eigenvalue changes, which a gap does not determine.

Balancing moves the eigenvalue, and the vector follows #

The reading makes a prediction that can be checked at once. Balanced is not symmetric found that balancing the Kac matrix — a diagonal similarity $B = D^{-1}AD$ chosen to even out its rows and columns — improves its computed eigenvalues by a constant factor without changing how fast they are lost with the order. If the Schur route’s eigenvector error is the eigenvalue’s error times the vector’s speed, it should improve by exactly the same factor, since the eigenvector family of $A$ is fixed whatever similarity is used to compute it.

It does. Run through the Schur form of the balanced matrix, with each eigenvector of $B$ multiplied back by $D$, the middle eigenvalues at order 48 are off by $2.7 \cdot 10^{-11}$ instead of $3.0 \cdot 10^{-9}$, and their eigenvectors by $2.2 \cdot 10^{-11}$ instead of $2.4 \cdot 10^{-9}$; at 64, $1.2 \cdot 10^{-8}$ and $9.8 \cdot 10^{-9}$; at 80, $2.2 \cdot 10^{-6}$ and $1.7 \cdot 10^{-6}$. Over every lossy eigenvector of the three orders, 136 of them, the ratio is 0.76 to 0.83 — the same speed, at errors a hundred times smaller. Balancing bought the eigenvalues a factor of 14 to 110 and the vectors exactly the same factor, and nothing else.

It also rescues the ends. The balanced route’s componentwise error is within 6 to 29 units of the dynamic range’s bound at orders 48 to 80, where the unbalanced route’s was a hundred to two thousand times past it: $D$ carries the vector’s range, so the product with $Z$ is formed among entries of comparable size and the small ones are no longer differences of large ones. Symmetry is worth more than precision found that the symmetrised Kac matrix, a similarity by a different diagonal, computes its eigenvalues exactly; on that route the eigenvalue error vanishes and with it, by this reading, the Schur route’s slide.

What one matrix does not show #

The Sylvester–Kac matrix has every gap equal, so it cannot distinguish a constant that comes from the gap from one that does not except through its value, and it has an eigenvector family with a closed form, which is what makes the speed computable. On a matrix with unequal gaps both effects would be present at once, and only a matrix with known eigenvectors could separate them. The Schur form here has no two-by-two blocks; for complex eigenvalues LAPACK’s back substitution solves small two-by-two systems, and how they carry an eigenvalue’s error is not measured. And the orders stop at 80, where the middle eigenvalues still have five digits; the Schur forms stay real to order 104, but the eigenvalues there have three digits, and an eigenvector that slides by a thousandth is no longer a small perturbation of the right one.

The refinement’s cost is not measured either. A solve with a tridiagonal matrix is cheap, and the Sylvester–Kac matrix is tridiagonal, so here one solve per eigenvector is negligible beside the Schur form; on a dense matrix it is a factorisation per eigenvector, $n^3/3$ operations each, and refining all $n$ of them costs more than the Schur form itself. That is why back substitution is the default and why a code that wants accurate vectors for a few eigenvalues and not all of them would do well to refine those few.

Still open: the left eigenvectors, and a matrix with unequal gaps #

The left eigenvectors. The left eigenvectors of the Kac matrix satisfy the transpose’s recurrence and are exact rationals too. From the Schur form they come from the transposed triangular system. The prediction with a sign is that they slide with their eigenvalue at a speed of their own — not 0.79 — and that the product of the two speeds, left and right, times the eigenvalue’s error, accounts for the error in the eigenvalue condition number computed from the two vectors to within a factor of two.

Unequal gaps. A tridiagonal matrix with the Kac matrix’s off-diagonal pattern but a varying diagonal has unequal gaps and still has eigenvectors given by a three-term recurrence, so their speed is computable. The prediction is that the Schur route’s error is the speed times the eigenvalue error on every eigenvector to within ten per cent whatever the gaps, and that one over the gap predicts it only where the speed happens to be near one over the gap — so that the slide, not the mixing, is the Schur route’s error wherever the eigenvalue is computed less accurately than the vector’s neighbours are separated.