The form a real matrix can reach

The algorithm the libraries actually run builds the QR algorithm and runs it on symmetric matrices. The form that makes it affordable reduces to Hessenberg first and measures what that buys. Both essays end by saying that the general case is harder, and neither says in what way — which is the kind of sentence this site is supposed to notice.

Three things are harder and they are separable. This essay is the first of them, and it is not about an algorithm at all. It is about what the answer looks like, and the answer looks different for a reason that has nothing to do with arithmetic.

The argument, which is one line #

Suppose a real orthogonal matrix Z existed with ZᵀAZ = T upper triangular and real. The eigenvalues of a triangular matrix are its diagonal entries, and a similarity does not change the spectrum. So the eigenvalues of A would all be real.

A real matrix can perfectly well have complex eigenvalues — the rotation by ninety degrees has i and −i — so no such Z exists. The obstruction is the field, not the algorithm. No amount of iteration, no better shift and no more careful arithmetic produces a real triangular form for a matrix whose spectrum is not real, and a routine that appeared to would be returning something false.

What is reachable is the real Schur form: A = ZTZᵀ with Z real orthogonal and T quasi-triangular — upper triangular except for 2×2 blocks on the diagonal, one for each conjugate pair. Each block has the pair as its own eigenvalues, and no real rotation removes it.

Over the complex numbers none of this arises: every square matrix has a complex Schur form that is genuinely triangular. The real form exists because real arithmetic is faster, because a real matrix’s eigenvectors and invariant subspaces are more naturally described over the reals, and because a routine that returned complex output for real input would force every caller into complex storage for a result that is mostly real.

Reading the blocks #

A 2×2 block

$$\begin{pmatrix} a & b \\ c & d\end{pmatrix}$$

has eigenvalues (a + d)/2 ± √(((a − d)/2)² + bc), and it stays in the form precisely when that discriminant is negative. When it is negative the pair is a ± iβ with a = (a + d)/2 and β = √(−disc), and both are read off the block directly.

When the discriminant is not negative the block holds two real eigenvalues and it should not be in the form at all — a rotation finishes it, and the site applies one. That is not tidiness. If 2×2 blocks with real eigenvalues were left in place, the number of blocks would depend on when the iteration happened to stop, and “the number of conjugate pairs” would become an artefact of a tolerance rather than a property of the matrix.

So the site standardises every such block, and then asserts the count.

The count, checked at four spectra #

The measurement that makes the argument above into evidence: build a matrix with a known spectrum and check that the number of blocks is the number of pairs.

Building it is the exact-ground-truth habit applied where the subject usually offers nothing. A general matrix’s eigenvalues are exactly the thing there is no closed form for — so the matrix is built from its spectrum instead. Each pair a ± iβ becomes the block with a on the diagonal and ±β off it, which is a rotation scaled by the modulus; each real eigenvalue is a diagonal entry. An orthogonal similarity then hides the block structure completely without moving a single eigenvalue.

Four spectra, and the count is exact in all four:

The third row is the control the argument needs, and it is the one that would catch a routine leaving blocks it could have split. With a real spectrum the form comes out genuinely triangular, every subdiagonal entry below 10⁻¹² of the norm, and the site asserts that separately.

The entry that does not go away #

The blocks are not nearly-converged. On the matrix drawn above the surviving subdiagonal entry is 2.0, against a Frobenius norm of 5.8.

That number is worth pausing on because every other convergence claim in this field is about a subdiagonal entry going to zero. The algorithm the libraries actually run measures the unshifted algorithm failing on ±1, and the failure there looks like this: the entry does not move. The entry in a 2×2 block does not move either, and the two situations are completely different. The first is an algorithm that cannot make progress; the second is an algorithm that has finished.

Telling them apart is not a matter of watching the entry. It is a matter of asking whether the block it belongs to has a real discriminant, which is a question with an answer.

And the block is exactly where the site’s rule earns its place. realSchur returns a factorisation — A = ZTZᵀ with Z orthogonal — so every figure in this family prints ‖A − ZTZᵀ‖/‖A‖ and ‖ZᵀZ − I‖, and the reader can see that the thing with a 2.0 in it is a similarity of the original matrix rather than a computation that has gone wrong.

What that rule caught #

It caught a real bug, on this build, and the bug is the site’s own subject.

Splitting a 2×2 block with real eigenvalues requires a rotation whose first column is an eigenvector of the block. The obvious form of that eigenvector is [b, λ − a] with λ = (a+d)/2 + sign(p)√disc, where p = (a − d)/2. Written that way, λ − a is −p + sign(p)√disc: a difference of two same-signed quantities, and when |bc| is small compared with p² — the ordinary case for a well-separated pair of real eigenvalues — √disc equals |p| to most of its digits and the subtraction throws them away.

The rotation built from a badly determined direction fails to annihilate the subdiagonal entry. Forcing that entry to zero afterwards then puts the residual straight into the factorisation, and ‖A − ZTZᵀ‖/‖A‖ went from 3.6·10⁻¹⁵ to 9.3·10⁻¹¹ at that one line, on a matrix where every other step was exact to rounding.

Nothing else would have reported it. The eigenvalues were right to eleven digits, the form looked correct, the block count was correct, and a routine that returns the right spectrum passes every check anybody would think to write about an eigenvalue routine. The badge is what showed it, because the badge is about the factorisation rather than about the answer.

The fix is the standard one and it is cancellation takes the answer’s lesson applied: use [p + q, c] with q = sign(p)√disc, which is the same direction — (p − q)(p + q) + bc = p² − disc + bc = 0 — and adds two same-signed quantities instead of subtracting them. The residual went back to 2·10⁻¹⁵.

Why the symmetric case never met any of this #

The symmetric field is not a special case of this one with the complications removed. It is a different situation, and three things are different at once.

The spectrum is real, so there are no blocks and the Schur form is the eigendecomposition.

The Schur vectors are the eigenvectors. For a symmetric matrix T is diagonal, so the columns of Z are eigenvectors and the decomposition delivers everything. For a general matrix T is triangular and the columns of Z span nested invariant subspaces but are not eigenvectors — extracting those is a further computation, and one that can be badly conditioned even when the Schur form is not.

Every eigenvalue is perfectly conditioned. Weyl’s theorem with no constant, which symmetry is worth more than precision measures. In the general case each eigenvalue has its own condition number and they differ by decades inside one matrix, which is a condition number for one eigenvalue’s subject.

The third is the one that changes what a computed spectrum means. In the symmetric case a backward stable algorithm gives eigenvalues accurate to u‖A‖ and there is nothing more to say. Here backward stability gives the exact eigenvalues of a nearby matrix, and how far those are from the eigenvalues of this matrix is a separate question with a separate answer per eigenvalue.

What the form is good for besides eigenvalues #

Reading a spectrum off the diagonal is the obvious use and it is not the main one. Three others, and each is a reason the decomposition rather than the eigenvalues is what a library returns.

Invariant subspaces, ordered. The first k columns of Z span an invariant subspace of A, for every k that does not cut a 2×2 block in half. And the blocks can be reordered by a sequence of swaps that are themselves orthogonal similarities, so a caller who wants the invariant subspace belonging to the eigenvalues in the left half-plane — which is what a stability analysis wants — can have it as an orthonormal basis, without ever computing an eigenvector.

That matters because the eigenvector basis can be terrible. A matrix whose eigenvectors are nearly parallel has an eigenvector matrix with a huge condition number, and any computation performed in that basis inherits it. The Schur basis is orthonormal by construction, so it has condition number one, and the price is that T is triangular rather than diagonal.

Matrix functions. Computing exp(A), or A^{1/2}, or log(A) by diagonalising is the textbook route and it fails on exactly the matrices where it matters, because it inherits the eigenvector basis’s conditioning and does not exist at all when A is defective. The standard reliable algorithms compute the Schur form and then evaluate the function on the triangular factor by a recurrence along the diagonals. The 2×2 blocks are handled directly. Nothing complex is formed for a real matrix with a real-valued function of it.

And solving Sylvester and Lyapunov equations, AX + XB = C, which is the same idea: reduce both matrices to Schur form, solve the transformed equation by back-substitution one column at a time, and transform back. Control theory runs on that equation, and the Bartels–Stewart algorithm which solves it that way is essentially the reason numerical libraries expose the Schur form as a first-class result rather than as an intermediate.

The pattern across all three is one this site has met before under a different name. A decomposition whose factors are orthogonal is a decomposition whose factors do not amplify anything, which is the whole argument for Householder QR over the two Gram–Schmidts in a reflection cannot stop being one, and the whole argument for the SVD over the normal equations in the road that squares the problem. The Schur form is the same trade at the eigenvalue problem: give up the diagonal, keep the orthogonality, and every later computation is better conditioned for it.

What is asserted here #

A = ZTZᵀ to 2·10⁻¹⁵, and the orthogonality of Z to 3.7·10⁻¹⁵, on the figure, at every position of the slider. The reconstruction is what makes the rest of the claims about this matrix rather than about some matrix.

One 2×2 block per conjugate pair, and no others, at four spectra including the all-real control.

Every eigenvalue is the one built in, worst error 1.2·10⁻¹⁵ across the whole spectrum — checked against the construction rather than against a second run of the same routine.

And T is quasi-triangular: nothing survives below the subdiagonal, which is a separate check from the block count and would catch an iteration that stopped early somewhere else.

The refusal: that the real Schur form of a matrix with a conjugate pair is triangular. It is fed the block’s subdiagonal entry and must reject, and the entry is separately asserted to be of the size of the matrix rather than of the size of a tolerance — 2.0 against a norm of 5.8, so a reader cannot read the refusal as a statement about a nearly-converged quantity.

Where the count could still be wrong, and how it would show #

The block count is asserted, and it is worth being explicit about the one situation in which it is a judgement rather than a fact — because a reader who takes it as a fact will eventually be surprised.

A pair of complex eigenvalues with a very small imaginary part and a pair of very close real eigenvalues look nearly identical to the arithmetic. The discriminant of the block is a small number computed by subtracting two larger ones, so its sign — which is the entire question — is decided by a quantity that has lost most of its significant digits. Nudge the matrix by rounding error and a pair that should be complex comes out real, or the reverse.

That is not a defect in the algorithm and it cannot be repaired by computing the discriminant more carefully. It is the eigenvalues themselves that are close together, and “are these two eigenvalues a conjugate pair or two real ones” is a question with no stable answer when they are within rounding of coinciding. The condition number of the eigenvalues at a near-defective point goes as the reciprocal of their separation, which is symmetry is worth more than precision’s eighth-root measurement in a simpler form.

The site’s slider walks straight through it. complex-pair-plane runs the imaginary part from 0 to 3, and at exactly 0 the pair is two real eigenvalues while at 0.25 it is a genuine pair — so the frames either side of the bottom of the range differ in their block count, correctly, and would differ unpredictably if the range started at 10⁻⁸ instead of at 0. The figure asserts the single-shift failure only above an imaginary part of 0.2 for that reason: below it, whether the matrix has a complex pair at all is not a question the arithmetic can settle.

Which is the honest summary of the whole essay. The number of 2×2 blocks is a property of the matrix wherever the spectrum is not nearly degenerate, and a property of the rounding where it is — and the interesting thing about that sentence is that it is the same sentence as rank is a decision’s, about a different integer. Rank is a decision about a gap in the singular values; the block count is a decision about a gap between two eigenvalues. Both are integers a floating-point matrix does not quite have, both are computed by comparing something against a tolerance, and in both cases the useful output is the integer and the quantity the decision was made on.