Two shifts that are never formed

The algorithm the libraries actually run establishes that the QR algorithm is unusable without a shift and cubic with one, and measures both: 248 iterations against 6 on a matrix with a spectrum of 8, 4, 2, 1.8, and outright failure on a matrix with eigenvalues ±1.

The shift it uses is Wilkinson’s, and Wilkinson’s shift is a real number. On a matrix whose trailing block has a conjugate pair, that real number is trying to approximate an eigenvalue which is not on the real line, and it cannot.

The measurement, before the algorithm #

The failure is worth running rather than describing, and the site runs it.

Take a 4×4 real matrix built from the spectrum {1 + 2i, 1 − 2i, 4, −1}. Give it to the single-shift algorithm with Wilkinson’s shift, 200 iterations, a tolerance of 10⁻¹². It does not converge. The subdiagonal entry it is chasing moves by 5.8·10⁻¹⁵ between the hundredth iteration and the two hundredth, which is to say it does not move at all.

That is the same signature as the ±1 failure in the symmetric field, and the cause is the same shape: the shift is stuck at a distance from the eigenvalues it is trying to reach, so the ratio governing convergence is 1. The difference is that on ±1 a shift fixes it in one step, and here no real shift exists that would.

The double shift settles the same matrix in seven double steps.

Shifting by both at once #

The idea is one sentence. If μ is a complex shift then so is μ̄, because the matrix is real and its complex eigenvalues come in pairs. Do both, one after the other:

$$(A - \mu I)(A - \bar\mu I) = A^2 - sA + tI, \qquad s = \mu + \bar\mu, \quad t = \mu\bar\mu$$

and both s and t are real — a sum of conjugates and a product of conjugates. So the operator being factorised is a real matrix, the QR factorisation of it is real, and the resulting similarity is real. The complex shift is used and never formed.

In practice s and t come straight from the trailing 2×2 block: s is its trace and t is its determinant, so the two shifts are that block’s eigenvalues and neither is ever computed. When the block’s eigenvalues happen to be real, the same step performs two real shifts at once and works perfectly well; nothing has to detect which case it is in.

That is the whole of the mathematical content, and it would be an unusable algorithm as stated. Forming A² costs O(n³), destroys the Hessenberg structure that the form that makes it affordable established is the only reason the iteration is cheap, and squares the condition number of everything in sight.

The implicit version #

What is done instead uses a theorem to avoid the product entirely.

The implicit Q theorem says that if two orthogonal similarities both reduce A to unreduced upper Hessenberg form and their first columns agree, then they agree up to the signs of their columns. So a step is determined by its first column alone.

The first column of A² − sA + tI is available for three multiplications, because A is Hessenberg: it has exactly three nonzero entries, and they are

$$x = a_{11}^2 + a_{12}a_{21} - s\,a_{11} + t, \quad y = a_{21}(a_{11} + a_{22} - s), \quad z = a_{32}a_{21}.$$

Build the 3×3 Householder reflector that maps [x, y, z] to a multiple of e₁ and apply it as a similarity. That is the first column right, and by the theorem the rest of the step is determined — so all that remains is to restore Hessenberg form without disturbing the first column, which fixes everything.

Restoring it is the bulge chase, and it is what the figure draws.

The bulge, measured #

Applying that first reflector to rows 1–3 and columns 1–3 puts entries where the Hessenberg form has none. The site locates them rather than describing them, and the shape is exact:

A three-entry triangle, moving down and right by exactly one each time, shrinking to a single entry and then gone. On a 6×6 that is five reflectors, which is n − 1, and each one costs O(n) work — so a whole double step is O(n²), the same as a single shifted step, and half as many of them are needed.

One of the five creates the bulge and the other four chase it, which is where the customary “n − 2 reflectors” comes from: it counts the chase and not the creation. Both counts are right about different things, and the figure shows the whole of it.

The positions are the part worth reading rather than the count. At reflector two the bulge is at (4,2), (5,2), (5,3); at reflector three it is (5,3), (6,3), (6,4). Every index has gone up by one — the same triangle translated one step down the diagonal — and each bulge overlaps the one before it in exactly one entry, the corner that the next reflector is built to annihilate.

The site asserts the leftmost column of the bulge is the reflector’s own index, which is the version of the claim with no off-by-one available to it. The first version of that check looked for “two entries below the subdiagonal” and measured the lowest row, which starts at 3 rather than 2 — the right-hand application of a reflector reaches one row further than the left-hand one, so the bulge is a triangle rather than a pair. The assertion was wrong about the shape of the thing it was asserting about, and the figure it produced would have been wrong in the same way.

So the bulge does not shrink gradually. It is three entries at reflectors one, two and three, one entry at four, and none at five — full size until it runs out of matrix, then gone in two steps. A reader expecting it to taper has the wrong picture of what a reflector does: each one annihilates exactly the entries the previous one created and creates exactly as many again, until there is nowhere left to create them.

Checked against the thing it avoids #

The implicit Q theorem is the licence for all of this, and citing a theorem is not this site’s habit. So both routes are computed — on a matrix small enough that forming the product is affordable — and required to agree.

The explicit route forms M = A² − sA + tI, factorises it by Householder QR, and applies the factor: QᵀAQ. The implicit route computes M’s first column, builds one reflector from three numbers, and chases a bulge.

They agree to 2.2·10⁻¹⁵ on a matrix of norm 5.8, which is rounding.

The comparison is made on absolute values entry by entry rather than as a norm of a difference, and that is not a fudge. “Essentially unique” in the implicit Q theorem means up to the signs of the columns, and the two routes make different sign choices in their reflectors — a Householder vector’s sign is chosen for numerical reasons, and the two constructions choose differently. A norm of the difference would be about 8, and the claim being checked is a claim about the matrices being the same up to those signs.

Being precise about that is worth more than taking a norm and finding it small. A check that would have passed for the wrong reason is not evidence, and a check that fails because the claim was stated imprecisely sends the next reader to break something that works.

Deflation, and the shift nobody advertises #

Two details that a description of the algorithm usually omits and a working one cannot.

Deflation. After each step the subdiagonal is scanned for entries small relative to their neighbouring diagonal entries; those are set to zero and the problem splits. Without it the iteration would keep working on a full matrix long after most of its eigenvalues had converged, and the cost would be multiplied by n. The site records where each deflation happened, which is also how the block count is read off at the end.

The exceptional shift. A double step whose shifts are the trailing block’s eigenvalues can stall, on a matrix with a symmetry the shift happens to share — the same failure the unshifted algorithm has on ±1, one level up. Every production implementation carries an ad-hoc shift for it: after some number of iterations without a deflation, replace s and t with values derived from the size of the subdiagonal entries rather than from the block, break the symmetry, and carry on.

That is genuinely a hack and it is in every library. It is worth naming for the same reason the swap that is not optional is worth writing: an algorithm presented without the part that stops it looping is an algorithm presented as more elegant than it is. And an implementation without one does not fail — it runs forever, which is worse.

Why the product would have been worse than slow #

The usual objection to forming A² − sA + tI is the cost, and the cost is the least of it. Three separate things go wrong and only one of them is arithmetic.

The structure. A is Hessenberg, so it has about n²/2 nonzeros. A² is full. Forming it therefore throws away the reduction that the form that makes it affordable measures as the difference between an O(n²) step and an O(n³) one, and it throws it away at every iteration rather than once.

The conditioning. Squaring a matrix squares its condition number, which is the same failure the road that squares the problem records about the normal equations.

The tempting way to say what that costs is in eigenvalue separation, and it is worth not saying it, because the arithmetic goes the other way. λ₁² − λ₂² = (λ₁ − λ₂)(λ₁ + λ₂), so squaring multiplies a relative gap by (λ₁ + λ₂)/λ₁, which for a close pair is two. Measured at gaps of 10⁻², 10⁻⁴, 10⁻⁶ and 10⁻⁸ the ratio is 2.000 every time: a gap of 10⁻⁶ becomes 2·10⁻⁶, not 10⁻¹². Squaring separates close eigenvalues slightly rather than merging them, and the damage the squaring does is to the condition number and not to the spectrum’s resolution.

And the shifts themselves are near-cancelling. When the double shift is doing its job, s and t are close to the trailing block’s own trace and determinant, so A² − sA + tI is close to singular by construction — that is the entire point, since a nearly singular operator is what makes the QR step converge.

How close is worth measuring, because the answer changes what the objection is. Tracking the smallest singular value of A² − sA + tI through a converging run on this essay’s own 4×4:

step   σₘᵢₙ                  step   σₘᵢₙ
  0    3.01                    4    1.51·10⁻⁶
  1    2.51                    5    1.09·10⁻²³
  2    1.14                    6    9.24·10⁻³⁶
  3    4.55·10⁻²               7    5.76·10⁻⁵¹

Singular to twenty-three digits by the fifth step and to fifty-one by the eighth. That is not an ill-conditioned matrix; it is a matrix whose smallest singular value has fallen below anything a double-precision entry can carry information about. Assembling it would produce a matrix all of whose entries are the size of the largest one, and whose smallest singular value is a number that the assembly’s own rounding has entirely replaced.

So the third objection is not “forming a nearly singular matrix and factorising it is bad practice”, which is a matter of degree. It is that after five steps the product cannot be represented, and the fast convergence the algorithm gets is exactly the thing that makes it so.

That last clause is worth sitting with, because it inverts the usual relationship between an algorithm’s speed and its numerical difficulty. Normally a method that converges faster spends less time near a hard problem. Here the convergence is the approach to a singular operator: each step makes the shift a better approximation to an eigenvalue, and a better approximation is a more nearly singular A − μI. The better the algorithm works, the less representable the object it is implicitly factorising becomes, and an implementation that formed that object would degrade in exact proportion to how well it was doing.

The implicit route is immune to all of it for one reason and it is not an accuracy argument. It never forms the product, so there is nothing whose smallest singular value can fall below the format. What it forms is three numbers — x, y and z, the first column — and one 3 × 3 reflector, and the reflector is a reflection that cannot stop being one whatever those three numbers do. The structure carries the accuracy, and the near-singular operator the whole step is about is never an object in memory at any point.

That is a strong form of a claim this collection makes repeatedly — that the inverse is never formed, that a determinant is a sign rather than a value, that a projection is a solve. Usually the argument is about cost, with accuracy as a bonus. Here it is the reverse: the cost saving is a factor of n, and the accuracy saving is the difference between an algorithm and one that stops working after five steps.

The implicit route sidesteps all three at once, and it does so without any accuracy argument of its own: it never forms the product, so there is nothing for the product’s conditioning to damage. The only thing the implicit route computes from the product is its first column, which is three entries built from six multiplications, and a reflector built from three numbers cannot inherit a condition number.

That is a pattern worth naming because it recurs. The cheapest route and the most accurate route are the same route, and they are the same route for the same reason — both are consequences of not forming an intermediate object that nobody wanted. Householder QR against the normal equations is the same shape; computing a residual rather than an inverse is the same shape; and an orthogonalisation nobody calls one makes the same observation about the Krylov basis. It is one of the few places in this subject where there is no trade to describe.

What is asserted here #

The bulge is three entries and moves one column per reflector, at every reflector but the last two, with the column index checked against the reflector’s own index.

A whole double step takes n − 1 reflectors, and the last one leaves nothing below the subdiagonal — which is the claim that the step restores Hessenberg form rather than merely reducing the damage.

The implicit step is the explicit one, entry by entry in absolute value, to 2.2·10⁻¹⁵.

A = ZTZᵀ and ZᵀZ = I, printed on every figure in this family, because realSchur returns a factorisation and this site does not draw one without them.

The refusal: that a single real shift converges on a matrix with a conjugate pair. It is run for two hundred iterations and refuses, and the drift of the entry it is chasing is separately asserted to be under 10⁻⁶ of the norm — so the refusal is about an entry that does not move rather than about one that is moving slowly.

What Francis’s paper actually changed #

It is worth being clear about which part of this was the contribution, because the double shift and the implicit step are usually described together and only one of them was new.

Shifting by a conjugate pair to stay in real arithmetic is the obvious move once the problem is stated, and it was known. What was not obvious is that the whole step can be performed without forming the shifted operator — that a reflector built from three numbers plus a mechanical restoration of Hessenberg form is provably the same similarity. That is the implicit Q theorem being used as an algorithm rather than as a fact, and it is what turned an O(n³)-per-iteration method into an O(n²)-per-iteration one.

The consequence is that the algorithm in LAPACK’s dhseqr today is recognisably the 1961 algorithm. It has acquired multishift versions, which chase several bulges at once so that the arithmetic can be done as dense matrix-matrix products and reach a useful fraction of a machine’s peak throughput; it has acquired aggressive early deflation, which detects convergence sooner by examining a trailing window rather than a single entry; and it has acquired a great deal of care about the exceptional shift. None of those change the step. They change how many are done at once and how quickly a converged eigenvalue is noticed.

That longevity is unusual in this subject and it is worth a sentence about why. The methods that get replaced are the ones whose cost model was wrong for the hardware — the normal equations, explicit inverses, Jacobi’s sweep over all pairs. The methods that survive are the ones whose inner loop is a small orthogonal transformation applied to a structured matrix, because that inner loop keeps being implementable well on whatever the machines turn out to be. Householder QR is the other example, and a reflection cannot stop being one makes the same observation about it from the accuracy side.

The bulge chase is the same object seen from the cost side: a sequence of tiny orthogonal transformations, each touching three rows and three columns, arranged so that the structure the expensive reduction created is never lost.