One number that has to be right

A reflection cannot stop being one explained why Householder QR holds ‖QᵀQ − I‖ at 10⁻¹⁵ whatever the condition number of A, where Gram–Schmidt loses orthogonality in proportion to it. The explanation was structural rather than a bound: the reflection’s matrix is I − 2vvᵀ with v a unit vector, so rounding v gives the reflection about a slightly different plane, and a reflection about any plane is orthogonal. Nothing in that sentence mentions A, which is why nothing about A can spoil it.

The sentence is right and it is not the whole of the structure. I − 2vvᵀ is orthogonal exactly when ‖v‖ = 1, and the site’s own implementation makes that true by dividing v through by its norm. Every production code does something else. LAPACK stores v unnormalised with its leading entry fixed at one, and carries a separate scalar β = 2/vᵀv beside it, because that saves a division per component and lets v live in the storage the zeros it creates have just freed. The two are the same reflection in exact arithmetic. They are not the same object numerically: one has a constant where the other has a computed number.

So the structural guarantee has a dependency, and it has exactly one. This essay finds it by breaking it.

Two degrees of freedom, not one #

Cancellation takes the answer is the other place a single arithmetic quantity decides a construction’s worth in this field, and the sign choice the earlier measurement measured is the same shape one step earlier. Write the reflection as H = I − βvvᵀ with no assumption on either. Then

by three lines of algebra, and the expression vanishes exactly when β = 2/‖v‖². That is one equation. The direction of v does not appear in it at all: v enters only through its norm, and the whole of v’s angular freedom — n − 1 parameters — leaves HᵀH exactly equal to I.

Write β = (2/‖v‖²)(1 + η) and the coefficient becomes 4η(1 + η)/‖v‖² times ‖v‖², which is 4η to first order. Write instead H = I − 2vvᵀ with ‖v‖² = 1 + δ and the coefficient is 4δ. It is the same expression under two names, because it is the same one degree of freedom: the product β‖v‖², which has to be 2.

The algebra says the departure is linear in the error with a constant of four. Computing HᵀH for a reflection built from an eight-component vector and reading the coefficient off the entry where |vᵢvⱼ| is largest gives 3.99966·10⁻¹², 3.99999917·10⁻⁹, 4.000004·10⁻⁶ and 4.004·10⁻³ at η = 10⁻¹², 10⁻⁹, 10⁻⁶ and 10⁻³, against the expansion’s 4.000·10⁻¹², 4.000·10⁻⁹, 4.000004·10⁻⁶ and 4.004·10⁻³.

The coefficient rather than the largest entry, because the largest entry of $vv^{\mathsf T}$ is $\max_i |v_i|^2$, which is below one for every vector that is not a coordinate axis. Reading the departure off an entry reports 2.63·10⁻⁹ where the algebra says 4.00·10⁻⁹, and the discrepancy is the vector rather than the arithmetic.

It is worth pausing on how few constraints that is. A reflector in eight dimensions is eight numbers, and it is stored with a scalar beside it, so nine numbers go into each reflection. Eight of the nine are free: any values at all produce a matrix that is exactly orthogonal, provided the ninth satisfies one equation. That is the whole of the structure the earlier measurement called structural, and it is a great deal of freedom — which is why perturbing a reflector’s direction, which is what rounding does to it, costs nothing.

Gram–Schmidt has no such arrangement. Its Q is assembled column by column from differences of vectors that are nearly parallel when A is ill-conditioned, and there is no equation that, if satisfied, makes the result orthogonal; orthogonality is what the subtractions are attempting rather than what the representation guarantees. Two Gram–Schmidts measured what changing one line of that attempt is worth and orthogonal is a number measured the attempt’s failure directly. The contrast this essay sharpens is between a method that computes orthogonality and a method that computes a direction and gets orthogonality for nothing.

What the direction is free to do #

Take the QR of an 8×8 matrix with condition number 10⁸, and inject a relative perturbation of ε into every component of every reflector — a different perturbation per component, so the direction genuinely moves rather than being rescaled. Then normalise, as the implementation does, and factorise.

Across ε from 10⁻¹⁴ to 10⁻², ‖QᵀQ − I‖ reads 1.99, 1.59, 1.57, 1.51, 1.34, 1.84 and 1.49 — all times 10⁻¹⁵, against 1.70·10⁻¹⁵ unperturbed. Twelve decades of perturbation and the quantity does not move.

That is not a small effect being drowned; the perturbation is doing exactly what it should. ‖A − QR‖/‖A‖ over the same range goes 5.15·10⁻¹⁵, 5.05·10⁻¹³, 5.05·10⁻¹¹, 5.05·10⁻⁹, 5.05·10⁻⁷, 5.05·10⁻⁵, 5.01·10⁻³ — linear in ε with a constant of a half, which is what a perturbed reflection reflecting the wrong column onto the wrong axis looks like.

So the perturbed factorisation is a genuinely different factorisation of a genuinely different matrix, and it is exactly orthogonal. That is the structural claim, and this is what it looks like when it is true.

And what the scalar is not free to do #

Perturb β instead, leaving v exactly as computed, and ‖QᵀQ − I‖ reads 6.61·10⁻¹⁴, 6.70·10⁻¹², 6.70·10⁻¹⁰, 6.70·10⁻⁸, 6.67·10⁻⁶, 6.29·10⁻⁴ and 6.49·10⁻² at the same seven values of ε.

Six point seven times ε, over eleven decades, on the nose. The constant is not four because the perturbations are independent per reflector and the seven of them accumulate, and it is not far off four for the same reason it should not be: the departure from orthogonality of a product of nearly-orthogonal matrices is about the sum of the individual departures, and seven reflectors with independent relative errors uniform on [−ε, ε] give a root-mean accumulation of a little under $4\varepsilon\sqrt{7}/\sqrt{3}$.

The residual moves too, at 3.08ε — smaller than the direction’s 5.05ε, and for a reason worth noticing: a wrong β scales the reflection by a factor near one rather than tilting it, so the column it was supposed to annihilate is left with a small residue rather than pointed somewhere else.

Which of the two the arithmetic actually gets wrong #

None of this is an argument against the β form, and the measurement that says so is the one a code cares about.

Both forms, at condition numbers from 10 to 10¹¹, give ‖QᵀQ − I‖ between 1.36·10⁻¹⁵ and 1.87·10⁻¹⁵, with no trend in either direction. The unit form’s one number is the norm it divides through by, computed as a square root of a sum of squares; the β form’s is 2 divided by the same sum. Both carry a relative error of a few units of roundoff, and four times a few units of roundoff is a few units of roundoff.

So in working precision the two are indistinguishable, and the flat line this field started from is flat for both. What the perturbation sweep establishes is not that one form is better. It is what kind of guarantee the flat line is: not “a reflection cannot stop being one” — it can, and 10⁻⁴ in one scalar is enough — but “the reflection has exactly one number that can be wrong, and that number is computed from a sum of squares to a relative accuracy that does not depend on A.”

That is a weaker statement and a more useful one, because it says where to look. Everything that could go wrong with a Householder QR’s orthogonality goes wrong in one place, and the place is the same for every reflector, at every step, on every matrix.

Where the sum of squares could go wrong #

Which invites the obvious question and answers part of it. vᵀv is a sum of m − k squares, computed in the working precision, so its relative error is bounded by about (m − k)u under the standard model. At m = 8 and double precision that is 10⁻¹⁵, four times which is 4·10⁻¹⁵, which is what the flat line reads.

Two things could make it worse and neither is exotic. A very long reflector — the first one on a tall matrix, whose length is m — makes the sum longer, so the guarantee degrades like m·u rather than being independent of the size. That is the one dependency the structural claim has on anything, and it is on the shape rather than on the conditioning.

And a sum of squares overflows when its terms are large, or underflows to zero when they are small, which is why every implementation of the norm scales before summing. An unscaled vᵀv on a vector whose entries are 10¹⁶⁰ returns infinity, β returns zero, and the reflection is the identity — orthogonal, and the wrong reflection. The units that overflow before the answer does is that failure in its general form, and it is the one place where the structural guarantee and the practical one come apart. The pivot that reads the units is the same hazard for a decision rather than for a constraint: a quantity computed in the matrix’s own units decides something, and the units are not part of the derivation that says the decision is right.

What this makes of “backward stable” #

The standard statement about Householder QR is that it is backward stable: the computed Q̂ and R̂ satisfy Q̂R̂ = A + ΔA with ‖ΔA‖ ≤ c(m,n)u‖A‖, and Q̂ is orthogonal to within c(m,n)u. Both halves are true and they are different kinds of true, which the perturbation sweep separates.

The residual bound is a bound: it involves a polynomial in the dimensions, it is proved by tracking rounding through the algorithm, and no structural argument gives it. The orthogonality bound is a bound on one scalar per reflection, and the polynomial in front of it is the accumulation over reflections rather than anything about the matrix. A reader of the two statements side by side would reasonably assume they have the same standing. They do not: the first is what the arithmetic does to an answer, and the second is what the arithmetic does to a constraint that the representation almost enforces by itself.

Achieved against enforced is the distinction the orthogonality essays have been drawing since the first of them, and this is where the line actually falls. Householder does not enforce orthogonality — a wrong β produces a non-orthogonal matrix and nothing complains. It reduces enforcing it to getting one sum of squares right.

What must fail for any of this to be wrong #

Four claims with refusals attached. That the algebra’s 4η is the arithmetic’s coefficient, checked by two routes at four values of η. That a one per cent perturbation of the direction leaves orthogonality at rounding — a claim that the direction carries the orthogonality fails. That the same perturbation of the scalar destroys it — a claim that a reflection survives an error in its scalar fails. And that the perturbed direction does move the factorisation — a claim that a different reflection produces the same factorisation fails, which is the check that stops the first one being vacuous.

That last refusal is the one worth having. A perturbation that changed nothing at all would leave ‖QᵀQ − I‖ flat for the uninteresting reason, and the residual curve is what distinguishes a free parameter from a parameter that was never applied.

Two numbers, and which one a code stores #

The choice between the two storage forms is made for reasons that have nothing to do with any of this, and it is worth setting them out because the measurement says the numerical argument is neutral.

The unit form divides v through by its norm — m − k divisions per reflector, so n(2m − n)/2 divisions over the factorisation, against the $n^3/3$ multiply-adds of the elimination itself. At m = n that is 3/(2n) of the work, which is the same fraction the pivot search costs Gaussian elimination and is negligible for the same reason. It then applies the reflection with the constant 2, which is exact.

The β form skips the divisions and stores one extra number per reflector. Its real advantage is not the arithmetic: it is that fixing v’s leading entry at one means v needs m − k − 1 numbers rather than m − k, and those fit exactly into the subdiagonal entries of the column the reflection has just zeroed. The reflectors live in the matrix, and the factorisation needs no workspace at all. That is worth more than a division count on any machine where memory traffic is the cost, which is every machine.

So a code stores β because of where v goes, not because of how β behaves, and the measurement here says the second consideration would not have decided it either way: 1.36·10⁻¹⁵ against 1.87·10⁻¹⁵ across eleven decades of conditioning is not a difference. What the sweep adds is that if a future representation put two computed numbers where this one has one — which is exactly what the blocked form does — the guarantee would be arithmetic rather than structure, and the constant in front of it would be worth measuring.

What this does not settle #

One size, 8×8, and one family — matrices built to a prescribed condition number from a random orthogonal pair. The m·u dependency argued above is not measured.

The perturbations are independent and uniform on [−ε, ε] per component, injected once per reflector. A perturbation correlated across components would move the direction less and the norm more, which is the same thing as perturbing β and is why the two experiments have to inject differently rather than differing only in a flag.

The measure is ‖QᵀQ − I‖ in the maximum-entry norm and the residual in the same, so both are within a factor of n of their spectral counterparts. The slopes reported are slopes and the constants are constants of this norm.

And the β form here computes β from the unnormalised v, which is LAPACK’s arrangement without LAPACK’s scaling. The scaled version’s error in β is smaller by whatever the scaling buys, and is not measured.

Why breaking a claim is how to test one #

A note on the method, because that measurement’s argument was of a kind this field makes often and this is the first time one of them has been tested rather than reasoned.

A structural argument says a quantity cannot go wrong. The evidence usually offered for it is that the quantity does not go wrong: a flat line across eleven decades of conditioning, which is what the earlier measurement drew. That evidence is consistent with the argument and it is consistent with several other things — with the quantity being insensitive for a different reason, with the experiment not varying whatever it is sensitive to, or with the argument being right about a mechanism that is not the operative one.

Perturbing the parts distinguishes them. If the structure is “v may be anything and the reflection is still orthogonal”, then perturbing v must leave orthogonality alone and must move something else; both halves are checked above, and the second is what rules out the perturbation simply not having been applied. If the structure has a dependency, perturbing the dependency must break orthogonality proportionally; it does, at a slope of one over eleven decades.

That is the same shape as the habit these measurements are made under: a rule is not believed until it has been fed something it has to refuse. A claim that has never been tested against something that should break it is a claim about which no measurement has been made, and a flat line is not a test — it is the absence of one.

Still open: the long reflector, the block, and the one number in a rotation #

The first reflector on a tall matrix. vᵀv is a sum of m terms at k = 0 and two at k = m − 2, so the guarantee should be worst at the first step and best at the last — the opposite end of the run from where the pivot decisions matter in elimination. Sweeping m at fixed n would say whether the flat line has a slope in the shape that the flat line in κ does not have.

The blocked form. A level-3 implementation does not apply reflections one at a time; it assembles a block into $I - YTY^{\mathsf T}$, where $T$ is triangular and every one of its entries is computed. That replaces one number per reflector with a triangle, and whether the guarantee survives is the next essay.

The rotation’s one number. A Givens rotation is a pair $(c, s)$ with $c^2 + s^2 = 1$, which is the same single constraint in two stored numbers rather than one. Whether the departure is linear in the violation of that constraint with the same constant of four, and whether renormalising c and s buys anything a code can measure, is the same experiment on the other orthogonal building block.

And what the sign does to Q. The earlier measurement noted that the reflection’s sign choice propagates into Q and is invisible in every quantity this field checks. The perturbation apparatus here would measure it directly — two factorisations differing only in a sign, both orthogonal, and the question is whether any downstream quantity can tell them apart.