A reflection cannot stop being one

Look again at the flat line.

Every other quantity here degrades as the problem gets harder. Forward errors rise with κ, loss of orthogonality in Gram–Schmidt rises with κ or κ², approximation errors rise, rank decisions get harder. The Householder line is flat, and the flatness is not a matter of the method being better by a large factor. It is a different kind of thing, and the reason is worth understanding because it generalises far beyond this algorithm.

The construction #

Each step of Householder QR takes a column x and reflects it onto a multiple of the first coordinate axis. The reflector is

H = I − 2vvᵀ, where v = (x − αe₁) / ‖x − αe₁‖ and α = ±‖x‖.

v is a unit vector by construction — it is divided by its own norm — and H is the matrix that reflects space in the hyperplane perpendicular to it.

The structural argument #

Now the point.

Suppose the arithmetic rounds while computing v. What comes out is v̂, a slightly different vector, and — after the normalisation, which is also rounded — a vector of length 1 ± u.

The matrix I − 2v̂v̂ᵀ is a reflection in the hyperplane perpendicular to v̂ rather than the one perpendicular to v — a slightly different plane, and still orthogonal, still norm-preserving, still an involution.

It is a reflection exactly when ‖v̂‖ = 1 exactly, and the normalisation is rounded like everything else. Writing ‖v̂‖² = 1 + ε,

HᵀH  =  I − 4v̂v̂ᵀ + 4v̂(v̂ᵀv̂)v̂ᵀ  =  I + 4ε·v̂v̂ᵀ

so ‖HᵀH − I‖ is exactly 4|ε| — four times the squared-norm error, which is eight times the norm error, and independent of the size of the matrix. Measured by perturbing the normalisation on purpose: a relative δ of 10⁻¹², 10⁻⁸ and 10⁻⁴ gives 8.000δ, 8.000δ and 8.002δ, and the same 8δ at four times the dimension.

That is a sharper statement of the structural argument rather than a hole in it. The departure from orthogonality is exactly the error in one normalisation and nothing else — nothing about the matrix, the column, or the condition number enters it. The invariant is maintained to the accuracy of the representation, which is what makes the line flat; it is not maintained exactly, which is what makes it a line at 10⁻¹⁵ rather than a line at zero.

The claim that it is the accuracy of the representation and nothing else is testable by moving the representation, and the figure has a slider for exactly that.

Across 20, 23, 26, 29, 32, 35, 38, 41, 44, 47, 50 and 53 bits Householder’s loss at κ = 10¹¹ reads 1.6·10⁻⁵, 1.8·10⁻⁶, 2.3·10⁻⁷, 3.6·10⁻⁸, 4.4·10⁻⁹, 4.9·10⁻¹⁰, 6.0·10⁻¹¹, 6.6·10⁻¹², 8.3·10⁻¹³, 1.3·10⁻¹³, 1.4·10⁻¹⁴ and 1.8·10⁻¹⁵. Divided by that format’s own unit roundoff, those twelve numbers are 16.8, 15.1, 15.4, 19.3, 18.9, 16.8, 16.5, 14.5, 14.6, 18.3, 15.8 and 16.2 — sixteen, give or take two, over a factor of 8.6·10⁹ in u.

So the structural argument is measured rather than argued: the departure from orthogonality is a fixed small multiple of the unit roundoff, at every precision and at every conditioning. Neither axis of the problem is in it. The 8δ of the algebra above is one reflector’s worth; sixteen is what n of them accumulate to on this matrix, and it stays sixteen when the arithmetic changes by nine orders.

The other line is the one that says what that is worth. Classical Gram–Schmidt at κ = 10¹¹ loses 5.5, 4.6, 4.5, 4.5, 4.4, 4.1, 3.5, 3.5, 3.2, 2.6, 2.5 and 2.4 over the same twelve precisions. Thirty- three extra bits of arithmetic — a factor of 8.6 billion in the unit roundoff — buy it a factor of 2.3, and it never gets within nine orders of being orthogonal at any of them.

That is the whole essay in one comparison. Precision buys Householder everything and buys the other method nothing, because one of them is spending it on a single normalisation and the other is spending it on subtractions that have already cancelled. A property enforced by a representation improves at the rate the representation improves; a property achieved by computation improves at the rate the computation’s worst step improves, which here is not at all.

The rounding perturbed which orthogonal transformation was applied, and perturbed the orthogonality only through the one quantity the form depends on. Orthogonality is not something the arithmetic computes entry by entry — it is a property of the form I − 2vvᵀ together with a single scalar constraint, and that form survives having its v rounded at a price of 8δ.

And the flat line is flat in κ, not in n #

The other thing the figure holds fixed is the size. ‖QᵀQ − I‖ for Householder QR at κ = 10⁸:

A factor of sixteen across a factor of thirty-two in n — about n^0.8, which is n reflectors each contributing its own normalisation’s rounding, adding up not quite as a walk and not quite in a line.

So the flatness the essay is about is flatness in the hardness of the problem, and the size still enters at a rate no orthogonal structure removes. Both halves matter to a reader deciding what to expect: a κ of 10¹¹ costs nothing and an n of 10⁴ costs three orders, and only one of those is on the figure.

The two halves also answer a question the structural argument raises and does not settle: if the property is enforced by the representation, why is there any loss at all? Because there is one number the representation does not enforce. A reflector is determined by a direction and a normalisation, and only the direction is free — the length is a constraint the arithmetic has to satisfy rather than a degree of freedom it can spend. Every rounding of that one scalar is a departure of 8δ, and n of them accumulate.

Which suggests the repair, and the reason nobody makes it. Renormalising v̂ exactly — computing ‖v̂‖ in a compensated pair and dividing by it — would drive δ below u² and the loss with it. It costs a few operations per reflector against the O(mn) the reflector is applied with, so it is nearly free, and it is not done because 10⁻¹⁵ was never the problem. The measurement above says what the repair would be worth if it ever became one: everything, because there is nothing else in the error.

Compare Gram–Schmidt. Its Q is assembled by subtraction and division, entry by entry, and nothing about the assembly enforces the property. When the subtractions are inaccurate, the result is a matrix, and there is no structure to fall back on. Its orthogonality is achieved, and it can therefore fail to be achieved.

Achieved against enforced #

That distinction is the transferable lesson, and it is worth stating outside the subject.

An invariant maintained by a representation survives perturbation. An invariant maintained by computation does not.

Store a rotation as an angle and it is a rotation no matter what the arithmetic does to the angle. Store it as a 2×2 matrix and repeated multiplication drifts it away from being orthogonal, which is why graphics and robotics code re-orthonormalises rotation matrices periodically, and why quaternions are popular: a normalised quaternion is a rotation by construction.

The same idea explains why the algorithms on this site are made of orthogonal steps rather than inverses. An orthogonal step has κ = 1 and therefore amplifies nothing, and it cannot cease to be orthogonal. Two properties, both structural, and together they are most of what makes an algorithm stable.

The sign, which is the one place it can go wrong #

There are two choices of α — plus or minus ‖x‖ — and they differ only in which of the two directions along the axis x is reflected onto. Mathematically it does not matter.

Numerically one of them is a cancellation. If x already points mostly along e₁ with x₁ > 0, choosing α = +‖x‖ makes the first component of v equal to x₁ − ‖x‖, a difference of two nearly equal positive numbers. The direction that defines the reflection is then computed from cancelled digits, and every operation afterwards inherits it — cancellation takes the answer is the mechanism.

Choosing α with the opposite sign to x₁ makes it a sum. Nothing cancels. The figure at the top of this essay reports what the unsafe choice would cost on that particular vector: v would come out at a small fraction of the size of its operands, and the digits would be gone.

The fix is one sign flip and it is in every serious implementation. It is in almost no derivation, because the derivation is correct either way, which is the recurring shape of identical algebra, different arithmetic.

What Householder costs #

For an m×n matrix with m ≥ n, Householder QR costs about 2mn² − 2n³/3 operations, which for a square matrix is 4n³/3 — roughly twice Gaussian elimination and about twice Gram–Schmidt.

That factor of two is the price of orthogonality that does not depend on the matrix, and for a dense factorisation it is worth paying without thinking about it. Two situations complicate the choice.

Q is often not needed explicitly. The reflectors are stored as their v vectors, and applying Q or Qᵀ to something costs O(mn) per vector without ever forming the matrix. LAPACK returns the reflectors, not Q, for exactly this reason. If Q is formed explicitly the extra cost is real and frequently unnecessary.

Column-at-a-time access. Krylov methods generate their basis vectors one at a time and cannot factorise a matrix they do not yet have. Householder does not fit that pattern, which is why Gram–Schmidt with re-orthogonalisation survives there.

For everything else — a least-squares fit, a rank decision, an explicit orthonormal basis — Householder is the answer and the choice is not interesting.

Givens rotations, briefly #

The other orthogonal building block deserves a mention because it shares the structural property and is used where reflections are clumsy.

A Givens rotation acts on two coordinates at a time, zeroing one entry with a plane rotation built from a cosine and a sine satisfying c² + s² = 1. Like a reflection, it is orthogonal by construction: round c and s, renormalise, and what remains is still exactly a rotation of a slightly different angle.

Rotations cost more than reflections for a dense factorisation — they zero one entry at a time rather than a whole column — and they win where the matrix is sparse or already nearly triangular, because they touch only two rows. Updating a QR factorisation when a row is added, which is the standard move in recursive least squares, is a sequence of Givens rotations.

Both are instances of the same idea, and the idea is the reason this field exists as a separate part of the site: orthogonal transformations are the only ones that can be applied repeatedly without watching the error.

What is asserted here #

The figure’s claims are checked as it is drawn, and each is a measurement with a tolerance bracketed rather than a round number chosen.

The reflection preserves length: ‖Hx‖ equals ‖x‖ to 10⁻¹³, checked as two routes to the same number rather than as an inequality. It lands the vector on the axis: the second component of the image is below 10⁻¹³. It is orthogonal: ‖HᵀH − I‖ below 10⁻¹⁵.

And it is an involution: HH must be the identity to 10⁻¹⁴. That is the second route, and it is the one that would catch a sign error in the construction of the reflector — a wrong sign gives a matrix that is still orthogonal and still norm-preserving and is not a reflection, so the first three checks would all pass.

The cancellation claim is a measurement too: the unsafe sign is computed, and the size of the resulting v relative to what it would be with the safe sign is required to show real cancellation. The essay does not warn about it; the build measures it.

Across the loss curve, Householder’s ‖QᵀQ − I‖ is required to stay below a bound proportional to the unit roundoff and the matrix size at every condition number, which is the flat line as an assertion. If it ever rose with κ — which is what a subtle bug in the reflector application would cause — the build would stop, and the essay’s central claim would be the thing that caught it.

The cost accounting, in full #

It is worth being precise about the arithmetic, because “Householder costs twice as much” is repeated without the qualifications that make it decision-relevant.

For an m×n matrix with m ≥ n, factorising costs about 2mn² − 2n³/3 flops. Classical or modified Gram–Schmidt costs about 2mn². For a square matrix that is 4n³/3 against 2n³, so Householder is cheaper — the “twice as expensive” folklore comes from the case where Q is formed explicitly, which costs a further 2mn² − 2n³/3 and is usually unnecessary.

The reflectors are stored as their v vectors, in the space the zeroed entries of A vacated, and applying Q or Qᵀ to a vector costs O(mn) by applying the reflectors in sequence. LAPACK returns them in that form for exactly this reason, and code that calls orgqr to materialise Q is frequently spending more than the factorisation cost on something it will only use to multiply.

So the honest summary: Householder is not more expensive than Gram–Schmidt for a dense factorisation, it is more expensive than a Gram–Schmidt that also produces Q explicitly, and it is the wrong shape for algorithms that generate their columns one at a time.

Why the same argument applies to rotations #

The structural property is not special to reflections, and the second orthogonal building block shows why the argument is about form rather than about a particular matrix.

A Givens rotation is built from a cosine and a sine with c² + s² = 1. Round c and s, renormalise, and what remains is exactly a rotation through a slightly different angle. Same guarantee, same reason: the property is enforced by the parameterisation rather than achieved by the arithmetic.

Rotations zero one entry at a time rather than a whole column, which makes them more expensive for a dense factorisation and cheaper for a sparse or nearly-triangular one, since they touch only two rows. That is why the QR algorithm for eigenvalues uses rotations on a Hessenberg matrix, and why updating a factorisation after a new row arrives — the standard move in recursive least squares — is a sequence of Givens rotations rather than a refactorisation.

The Jacobi methods used throughout this site are the same object again: each sweep is a sequence of plane rotations, so the accumulated eigenvector matrix cannot drift away from orthogonality. Symmetry is worth more than precision reports ‖VᵀV − I‖ = 10⁻¹⁵ for the Hilbert matrix’s eigenvectors, and that number is flat for the same structural reason the Householder line is flat.

What “structural” buys elsewhere #

The general lesson is worth taking out of the subject, because it applies wherever a computation must preserve a property.

Represent a rotation as three angles and it is a rotation no matter what happens to the angles. Represent it as a 3×3 matrix and repeated composition drifts, which is why graphics and robotics code re-orthonormalises periodically and why unit quaternions are popular — a normalised quaternion is a rotation by construction, and normalising is one division.

Represent a probability distribution as unnormalised weights plus a normalisation step and it sums to one exactly, every time. Represent it as numbers that are supposed to sum to one and it drifts.

Represent a symmetric matrix by its upper triangle and it is symmetric. Store all n² entries and compute with them, and the asymmetry that rounding introduces is enough to make a symmetric eigensolver refuse — which is exactly the precondition this site’s Jacobi routine enforces, and exactly why (A + Aᵀ)/2 appears in careful code.

In each case the pattern is the same: choose a representation in which the invariant cannot be violated, and the arithmetic can be as sloppy as it likes.

Where this leaves the field #

Three essays on orthogonality, and a single practical conclusion.

Orthogonal is a number established that the property is a measurement, that it can be badly violated by a factorisation that reconstructs its matrix perfectly, and that unit columns are not perpendicular columns. Two Gram–Schmidts established that one word of difference produces eight orders of magnitude, and that the R factors agree while the Q factors do not. This essay establishes why one method escapes the problem entirely.

The conclusion: use Householder for dense factorisations, use modified Gram–Schmidt with re-orthogonalisation where the access pattern requires it, and measure ‖QᵀQ − I‖ either way. The first two are free. The third costs a matrix product and is the only thing that would reveal whether either had gone wrong.

Which is the honest summary of the trade. Elimination is a sequence of choices is faster and relies on an observation; orthogonal is a number measures what the alternative buys; and the projection and the right angle is where the guarantee is actually spent.

What this does not protect against #

Being precise about the guarantee, because “orthogonality is structural” can be over-read.

Householder QR is backward stable: the computed Q and R satisfy QR = A + E with ‖E‖ at the level of rounding, and the computed Q is within rounding of an exactly orthogonal matrix. Both of those are theorems.

Neither says the factorisation is accurate in the sense of being close to the exact QR of A. For an ill-conditioned matrix the exact R and the computed R can differ substantially, because R inherits the conditioning of A — and it must, since R and A have the same singular values. What is guaranteed is that the computed pair factorises something very close to A, which is the strongest statement available and is not the same statement.

The distinction is exactly the one the exact answer to a nearby problem sets up, applied to a factorisation rather than a solution. The structural argument in this essay is about one factor’s orthogonality; the conditioning of the problem still governs how much the factors themselves can be trusted, and orthogonal steps guarantee only that the algorithm adds nothing to it.

That is worth saying because the flat line in the first figure invites the opposite conclusion. It says the method does not degrade. It does not say the answer improves.

The reflection’s sign, which nothing checks #

A Householder step can reflect onto plus or minus the axis, and both are reflections. The choice is made for numerical reasons, it propagates into Q, and it is invisible in every quantity this site checks.