Numerical Linear Algebra
where the algebra stops being the arithmetic
Eigenvalues, singular values, rank

The form that makes it affordable

One Householder reduction, done once, turns every subsequent iteration of the eigenvalue algorithm from cubic to quadratic cost. It changes no answer at all, which is why it is easy to describe as an optimisation and wrong to.
16 min read 7 figures Where the algebra stops being the arithmeticMeasured, not assertedTwo routes to a numberIdentical algebra, different arithmetic

The algorithm the libraries actually run is about the ingredient that makes the QR algorithm converge. This one is about the two that make it affordable, and they are a different kind of thing entirely.

A shift changes the answer’s arrival time and can change whether it arrives at all. The two steps here change nothing about the answer whatsoever. Run the algorithm with them and without them, and the eigenvalues that come out agree to 10⁻¹⁰ and the iteration counts are 17 and 18.

Which is exactly why they are worth an essay. A step that changes no answer and no iteration count, and without which the method is unusable, is a step whose value is invisible in every measurement the method reports.

Householder reduction to Hessenberg form, on a symmetric matrixTwo matrices side by side. The first is full; the second has zeros everywhere below the first subdiagonal, and being symmetric is tridiagonal.-0.16-0.63-0.493.43-0.025-0.632.3-0.69-0.95-1.70.076-0.49-0.694.3-1.6-1.41.93.4-0.95-1.65.40.0141.63-1.7-1.40.0144.60.055-0.0250.0761.91.60.0553.1A, symmetric→-0.164.600004.65.92.500002.51-1.20000-1.24.90.3800000.385.80.2400000.242H = QᵀAQ, tridiagonal‖A − QHQᵀ‖/‖A‖1.1·10⁻¹⁵below the subdiagonal0worst eigenvalue movement7.1·10⁻¹⁵a similarity, so the spectrum is untouched — and every later step is O(n²) rather than O(n³)one reduction, then every iteration is cheapthe eigenvalues did not move
Fig. 1 A symmetric matrix and its Hessenberg form, which for a symmetric matrix is tridiagonal. The badge carries the reconstruction residual and the largest distance any eigenvalue moved. Nothing about the spectrum changed; everything about the cost of the next step did.

What the reduction does

Apply Householder reflections as similarities — P A Pᵀ rather than P A — chosen to annihilate everything below the first subdiagonal, one column at a time.

The result is upper Hessenberg: one nonzero subdiagonal, structural zeros everywhere below it. On a 6×6 matrix that is 21 structural zeros, asserted as an exact count rather than as a norm, because the claim is about structure and a structure claim should be checked as one.

And a symmetric matrix comes out tridiagonal. The reduction is a similarity, so symmetry survives it, and a symmetric Hessenberg matrix has zeros above the first superdiagonal as well. The site checks ‖H − Hᵀ‖ below 10⁻¹² rather than assuming it.

The cost is O(n³), paid once — about 10n³/3 operations for a general matrix, or 4n³/3 for a symmetric one, which is comparable to a single LU factorisation and is therefore not the expensive part of anything.

Why the algorithm does not undo it

The property that makes the whole arrangement work, and it is not obvious.

A QR step preserves Hessenberg form. Factorise a Hessenberg matrix and Q is Hessenberg; multiply R by it and the result is Hessenberg again. So the reduction is done once and every subsequent iteration starts from the form it produced.

That is what converts the cost. A QR factorisation of a dense n×n matrix is O(n³). A QR factorisation of a Hessenberg matrix is O(n²) — there are only n−1 subdiagonal entries to annihilate, one per column, each a Givens rotation acting on two rows.

For a run of a few dozen iterations that is the difference between O(kn³) and O(n³ + kn²), and since k is comparable to n in practice it is the difference between n⁴ and n³. On a 1,000×1,000 matrix that is a factor of a thousand.

The reduction does not make the algorithm faster. It makes it possible.

The same statement holds for the storage. A dense QR step needs workspace proportional to n²; a Hessenberg one needs the n−1 rotations, which is O(n). On the matrices this is run at, that is the difference between an algorithm whose memory grows with the square of the problem and one whose memory is the matrix.

The two routes, again

The claim that a similarity preserves eigenvalues is a theorem, and this site checks theorems on the objects it actually computes with.

Six eigenvalues, computed from the matrix before reduction and from the matrix after, by the foundation’s Jacobi solver. The largest movement is below 10⁻¹².

That check is doing real work. The reduction is thirty lines of index arithmetic applying reflections to rows and to columns, and getting one of the two applications wrong produces a matrix that is Hessenberg, looks entirely reasonable, and is not similar to the original. Every eigenvalue would be wrong and nothing in the shape of the output would say so.

The reconstruction residual is the second guard: ‖A − QHQᵀ‖/‖A‖ = 4.5·10⁻¹⁶, on the badge. That one catches a Q that has drifted from orthogonality, which is a reflection cannot stop being one’s subject — Householder reflections are orthogonal by construction rather than by normalisation, which is why the accumulated Q survives dozens of applications.

Why it stops at one subdiagonal and not at zero

The obvious question, and its answer is the reason an eigenvalue algorithm needs an iteration at all.

If Householder similarities can annihilate everything below the subdiagonal in O(n³), why not annihilate the subdiagonal too and read the eigenvalues off a triangular matrix directly?

Because they cannot. A Householder reflection acting on the left annihilates entries in a column; as a similarity it must also act on the right, and that action fills the entries back in. Reducing to Hessenberg works precisely because the reflection is chosen to leave the first subdiagonal alone — the right-multiplication then cannot reach the zeros just created. Aim one column further and the right-hand action undoes the left-hand one.

There is a deeper reason it cannot work, and it is worth stating because it explains why the whole subject is shaped as it is. A finite sequence of similarity transformations computed by a formula would be an algebraic procedure producing the eigenvalues of a general matrix — the roots of its characteristic polynomial. Abel’s theorem says no such procedure exists for degree five and above.

So the eigenvalue problem is necessarily iterative, not because nobody has found the direct method but because there is none to find. Hessenberg reduction goes exactly as far as a finite process can, and the iteration takes over at the boundary. That is why the reduction and the iteration are separate phases with different costs, and why the reduction is not simply the first few steps of the algorithm.

Deflation, the second free win

When a subdiagonal entry becomes negligible, the matrix splits. The eigenvalue in the corner is converged, and the remaining problem is the leading block — strictly smaller, and the algorithm continues on it alone.

Without this every iteration would cost the full O(n²) however much of the matrix had already finished. With it the active block shrinks by one each time an eigenvalue converges, so the total cost is a sum over shrinking blocks rather than n times the largest.

Deflation is also what makes the convergence rate measurable at all, and this is where it caused trouble. The recorded history of subdiagonal magnitudes runs on past each deflation into the next subdiagonal, which starts large again — so a rate fitted to the concatenated history measures whichever convergence contributed the most points. At an eigenvalue ratio of 0.3 the fit came back 0.5, which is the ratio of a different pair.

The library now returns the history split per deflation. A flat array was hiding the fact that it described several separate convergences laid end to end, and the fix was a data-structure change rather than a numerical one.

Driving the subdiagonal to zero, with λ₄/λ₃ = 0.90A semi-logarithmic plot of the magnitude of the subdiagonal entry against iteration count for three shift strategies. The unshifted curve is a straight line; the two shifted curves plunge to the bottom of the plot within a few steps.0357010514017521024528010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration|subdiagonal entry|no shiftRayleighWilkinsontwo routes to one raterate, from the spectrum0.9rate, measured0.9iterations, none / Wilkinson45symmetric 4×4, spectrum 8, 4, 2, 1.8the dashed line is the prediction
Fig. 2 What the reduction and the deflation are in service of. The subdiagonal driven to zero, one at a time — each of those curves is one deflation, and the algorithm moves up to the next subdiagonal when the one below it is negligible.

What the reduction costs in accuracy

Nothing measurable, and the reason is worth having because it is not automatic.

A reduction to Hessenberg form is a sequence of about n Householder similarities, each applied twice — once from the left and once from the right. That is roughly 2n orthogonal transformations accumulated onto the matrix, and a reader who has been through two Gram–Schmidts should expect that to cost something.

It does not, and the reason is the same one that essay ends with. A Householder reflection is orthogonal by construction: it is I − 2uuᵀ for a unit u, and it is orthogonal for any u whatsoever, including a u that rounding has perturbed. There is no invariant being maintained approximately, so there is nothing for rounding to erode. The measured ‖QᵀQ − I‖ for an accumulated Householder Q on this site is 1.4·10⁻¹⁵ across eleven decades of condition number — flat, and a reflection cannot stop being one is the essay for why.

So the reduction is backward stable: the computed H is exactly the Hessenberg form of a matrix within rounding of A. Which is the strongest thing available and is what licenses treating the reduced problem as the original one — and it is why the eigenvalue-movement check comes back at 10⁻¹² rather than at something that would need explaining.

Had the reduction been built from Gram–Schmidt orthogonalisations instead, that number would degrade with the conditioning of the matrix, and the essay would be about a trade rather than about a free win.

Why a symmetric matrix is so much cheaper

The tridiagonal form is where the symmetric eigenvalue problem separates from the general one, and the separation is large.

A tridiagonal matrix has 3n − 2 nonzero entries. A QR step on it costs O(n) rather than O(n²), because there are n−1 rotations each acting on a constant number of nonzeros. So the symmetric algorithm is O(n³) for the reduction plus O(n) per iteration, against the general case’s O(n³) plus O(n²).

That is the cost side of the trade symmetry is worth more than precision describes from the accuracy side. There, symmetry buys eigenvalues that move by at most the size of a perturbation, with no constant, where a defective matrix’s move by a fractional power. Here it buys an order of magnitude in the exponent.

Symmetry is not a convenience that makes the analysis tidier. It changes the accuracy guarantee and it changes the complexity, and a problem that can be posed symmetrically should be.

What is asserted here

The reduction produces Hessenberg form, with the norm below the first subdiagonal under 10⁻¹³.

A symmetric matrix comes out tridiagonal, ‖H − Hᵀ‖ below 10⁻¹².

QHQᵀ is the matrix it reduced, to 4.5·10⁻¹⁶.

Every eigenvalue survived, largest movement below 10⁻¹² — the claim that a similarity is a similarity, checked on the computed object.

The reduced matrix is structurally zero below the subdiagonal, as an exact count of 21 positions on a 6×6 — not a norm, a count.

And both routes get the same eigenvalues, reduced and unreduced, to 10⁻¹⁰, in 17 and 18 iterations.

The last measurement is the interesting one

17 and 18. The reduction saves one iteration out of eighteen.

That is the essay’s real content and it is why the check was written to report iteration counts at all. Anybody measuring the value of Hessenberg reduction by timing iterations, or counting them, or comparing accuracy, finds almost nothing — the two routes agree on every quantity the algorithm reports.

The saving is entirely in what an iteration costs, and the algorithm does not report that. The site’s assertion is therefore about the structure rather than about a timing: 21 structural zeros below the subdiagonal, which is the thing that makes each Givens rotation touch a constant number of entries instead of a row.

That is an unusual shape for an assertion on this site, where nearly everything is checked by measuring a number the computation produced. Here the number the computation produces is the same either way, and the only honest check is on the structure that governs a cost the computation never mentions.

What “negligible” means, and who decides

Deflation needs a threshold, and the threshold is a judgement rather than a constant — which is the same shape rank is a decision is about, arriving in the iteration rather than in the spectrum.

The site’s implementation deflates when a subdiagonal entry falls below a tolerance times the Frobenius norm of the whole matrix. That scaling matters: an absolute threshold would be meaningless on a matrix whose entries are all 10⁻⁶, and would deflate immediately on one whose entries are 10⁶.

LAPACK’s criterion is finer. It compares the subdiagonal entry against the sum of the magnitudes of its two neighbouring diagonal entries, which is a local scale rather than a global one, and it makes a real difference on a graded matrix — one whose entries span many decades — where a subdiagonal entry that is negligible relative to the whole matrix may be entirely significant relative to the 2×2 block it sits in. Deflating it there would discard a small eigenvalue and replace it with noise.

That is the same failure mode as a rank decision made with a tolerance drawn from the arithmetic rather than from the data: a threshold chosen at the wrong scale, producing an answer of the right shape. And it has the same fix, which is to set the threshold from something local to the quantity being judged rather than from a global summary.

The general shape

Worth extracting, because it recurs and it is a genuinely hard thing to notice in one’s own code.

A step that changes no output cannot be validated by comparing outputs. Hessenberg reduction, the symbolic phase of a sparse factorisation, an elimination ordering — all three change the cost by orders of magnitude and change the answer by nothing at all. Two ends of the same arrow is the same observation in the sparsity field: two factorisations, one dense and one with no fill, both exact to 10⁻¹⁶, and the only difference is a cost nothing in the answer records.

The consequence for testing is direct. A test suite that compares answers is blind to this entire class of decision, and it will stay green while the code becomes unaffordable. What has to be asserted instead is the structure — the count of zeros, the fill total, the bandwidth — which requires knowing in advance which structure the cost depends on.

Where the same form appears elsewhere on this site

Hessenberg is not a curiosity of the eigenvalue problem. It is the shape a Krylov method produces whether or not anybody asked for it, and noticing that the two are the same object is worth the paragraph.

Arnoldi — the orthogonalisation inside GMRES — builds an orthonormal basis V of the Krylov subspace and, alongside it, a small matrix H satisfying AV = VH. That H is upper Hessenberg, for the same structural reason: each new basis vector is orthogonalised against all the previous ones and against nothing after them, so the coefficient matrix has one subdiagonal and nothing below.

The consequence in that setting is the same as here. A least-squares problem with a Hessenberg coefficient matrix is solved by n Givens rotations, one per column, and adding a step costs one more rotation rather than a fresh factorisation — which is how the spectrum that predicts nothing can plot a residual at every GMRES step without ever forming the solution.

And in the symmetric case Arnoldi’s Hessenberg matrix is tridiagonal, exactly as the reduction here produces a tridiagonal form, which is what turns Arnoldi into Lanczos and the long recurrence into a three-term one. One subdiagonal is what makes an O(n³) step an O(n²) one, in both places, for the same reason — and the symmetric collapse to tridiagonal is what makes it O(n) in both places too.

One Householder reflection, and the sign that keeps it stableA vector, the mirror line through the origin, and its reflection landing on the negative first axis at the same distance from the origin.the mirrorx, length 4.000Hx = (-4.000, 0)v = x − αe₁safe sign: α = −‖x‖, so v is formed from a sum and nothing cancelsunsafe sign: α = +‖x‖ gives ‖v‖ only 35.1% of ‖x‖ + ‖x‖ — the digits go‖HᵀH − I‖5·10⁻¹⁶‖Hx‖ − ‖x‖8.9·10⁻¹⁶second component2.2·10⁻¹⁶built from a unit vectororthogonality is structural
Fig. 3 The operation the reduction is built from. A reflection is orthogonal by construction, not by normalisation, which is why a chain of them can be applied as similarities without the accumulated transformation drifting.
How far a perturbation of size ε moves an eigenvalue, 8×8A log–log plot of eigenvalue movement against perturbation size. The symmetric case lies on a line of slope one; the non-symmetric case lies on a line of slope one eighth, and at a perturbation of ten to the minus sixteen it has already moved by a hundredth.10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/8)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision
Fig. 4 And what none of this can help with. The reduction preserves the eigenvalues of the matrix as stored, exactly; whether the matrix as stored has the eigenvalues of the matrix intended is a question about conditioning, and for a defective matrix the answer is that it does not.
GMRES on the Laplacian and on the cyclic shift, both 12×12A semi-logarithmic plot of relative residual against iteration. One curve falls steadily; the other is flat at one for every step until the last, where it drops to zero.02468101210⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1iteration‖r‖ / ‖b‖Laplaciancyclic shiftno progress at allevery eigenvalue of the shift is on the unit circleand it predicts nothing
Fig. 5 The same Hessenberg form, in the iterative field. Arnoldi produces one whether or not anybody asked, and it is what lets this figure plot a residual at every step without forming the solution.
Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 6 Why the reduction costs no accuracy. Two n-fold sequences of Householder similarities are applied and the orthogonality does not degrade, because a reflection is orthogonal by construction rather than by an invariant rounding can erode.
Jacobi sweeps on the 6×6 Hilbert matrixA plot of the off-diagonal norm against sweep number, falling from about one to ten to the minus seventeen in five sweeps, with the matrix shown at three stages beneath it.01234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweep‖off-diagonal‖eigenvalues, largest first1.6190.24240.016326.157·10⁻⁴1.257·10⁻⁵1.083·10⁻⁷κ = 1.5·10⁷‖VᵀV − I‖ = 1.9·10⁻¹⁵1.000.500.330.250.200.170.500.330.250.200.170.140.330.250.200.170.140.130.250.200.170.140.130.110.200.170.140.130.110.100.170.140.130.110.100.09as given1.620.050.00-0.040.00-0.000.050.240.000.010.00-0.000.000.000.000.00-0.00-0.00-0.040.010.000.020.000.000.000.00-0.000.000.00·-0.00-0.00-0.000.00·0.00after one sweep1.62······0.24······0.00······0.02······0.00······0.00after 8 sweepseach sweep is n(n−1)/2 rotationsthe trace is conserved
Fig. 7 The method that needs neither the reduction nor the shift. Jacobi’s off-diagonal norm falls without either, at a cost per sweep that no library is willing to pay — which is what the two structural steps in this essay are buying.

Where this goes next

0357010514017521024528010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration|subdiagonal entry|no shiftRayleighWilkinsontwo routes to one raterate, from the spectrum0.9rate, measured0.9iterations, none / Wilkinson45symmetric 4×4, spectrum 8, 4, 2, 1.8the dashed line is the prediction The algorithm the libraries actually run Factorise, multiply the factors back in the other order, repeat. That description is complete and correct and produces something nobody would use — on a matrix with eigenvalues +1 and −1 it does not converge at all, and the subdiagonal entry does not move by so much as a rounding error. the mirrorx, length 4.000Hx = (-4.000, 0)v = x − αe₁safe sign: α = −‖x‖, so v is formed from a sum and nothing cancelsunsafe sign: α = +‖x‖ gives ‖v‖ only 35.1% of ‖x‖ + ‖x‖ — the digits go‖HᵀH − I‖5·10⁻¹⁶‖Hx‖ − ‖x‖8.9·10⁻¹⁶second component2.2·10⁻¹⁶built from a unit vectororthogonality is structural A reflection cannot stop being one Householder QR holds orthogonality at 10⁻¹⁵ whatever the condition number of the matrix, and Gram–Schmidt does not. The reason is not that it is more careful. It is that its Q is built from unit vectors, and rounding a unit vector gives a different reflection rather than a broken one. 10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/8)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision Symmetry is worth more than precision A symmetric matrix gives up its eigenvalues to full accuracy however ill-conditioned it is. An unsymmetric one can move them by the eighth root of a perturbation, so the rounding involved in merely storing the matrix shifts the spectrum by a hundredth.