A triangle where the scalar was

One number that has to be right reduced Householder’s orthogonality to a single quantity per reflection. The reflector’s direction is free — perturb every component by a relative one per cent and ‖QᵀQ − I‖ does not move off 10⁻¹⁵ — and the scalar is not: perturb it by the same amount and the departure is six per cent, linear in the perturbation with a constant near four.

That is the unblocked algorithm, which no production code runs. Applying reflections one at a time is a sequence of matrix–vector products, and a matrix–vector product moves as many numbers as it multiplies. Every level-3 implementation therefore gathers a block of r reflectors into a single matrix and applies it with matrix–matrix products, using the compact WY form

with the reflectors as the columns of Y and T upper triangular of size r. The identity is exact: the product of r reflections is exactly of that shape, and T is determined by Y and the scalars.

Determined, and computed — by a recurrence whose inputs at column i are the columns of T already computed.

What the triangle holds #

A stable block is not a stable basis is the other blocked orthogonalisation this field measures, and reading the two together is the point of the section below. The recurrence is three lines. $T_{11} = \beta_1$; for each $i$, form $w = Y_{1:i-1}^{\mathsf T} v_i$, set the first $i-1$ entries of column $i$ to $-\beta_i T_{1:i-1,\,1:i-1}\,w$, and put $\beta_i$ on the diagonal.

So column i is a triangular solve against the block already built, scaled by the reflector’s own β. Every entry of T above the diagonal is a computed number that depends on every entry above and to the left of it, which is a different dependency structure from the unblocked form’s — where each β depends on its own sum of squares and on nothing else.

The counts are the point. A block of r reflectors stores r(r+1)/2 numbers in T:

Sixteen reflections applied one at a time have sixteen numbers that have to be right. The same sixteen applied as one block have a hundred and thirty-six.

It is worth being clear that the identity itself is not in doubt. The product $H_1H_2\cdots H_r$ of $r$ reflections, each of the form identity minus $\beta_i v_i v_i^{\mathsf T}$, expands to the identity minus a sum of terms each of which is a product of some of the $v_i v_i^{\mathsf T}$ with some of the scalars $\beta_i$, and every one of those terms has the form (a column of $Y$)(a row of $Y^{\mathsf T}$) with a scalar between — so the whole sum is $YSY^{\mathsf T}$ for some $S$, and that $S$ is triangular because term $i$ cannot involve reflector $j > i$. The recurrence is that expansion read off in order.

What the expansion does not carry across is the invariance. Each $H_i$ is orthogonal for a reason internal to itself; a product of orthogonal matrices is orthogonal because each factor is; and I − YTYᵀ is orthogonal because T happens to be the matrix the expansion produces. Given an arbitrary T there is no reason for it to be anything, and the set of T for which it is orthogonal is a surface of dimension r(r+1)/2 minus r(r+1)/2 — a point, once Y is fixed.

And what that costs #

Over six matrices of 24 rows and 16 columns at condition number 10⁸, ‖QᵀQ − I‖ for the assembled Q reads

A factor of five, over a sixteen-fold increase in the block and a hundred-and-thirty-six-fold increase in the stored numbers. That is the answer a code wants and it is the uninteresting half of the measurement: blocked Householder is stable, everybody knew it was, and nothing here is a warning.

The growth is roughly like the square root of the block — 7.8·10⁻¹⁶ times 4 is 3.1·10⁻¹⁵ against a measured 3.9·10⁻¹⁵ at r = 16 — which is what an accumulation of independent errors over a triangular solve of depth r looks like, and is far below the r(r+1)/2 that counting numbers would suggest. The entries of T are not independent; the recurrence’s triangular solve is against a matrix whose diagonal entries are the βs, all of order one, so the conditioning of the solve does not grow with the block.

The guarantee has changed shape #

The interesting half is what happens when T is wrong. Perturb every entry by a relative 10⁻⁶ and the same five block sizes give 2.27·10⁻⁶, 2.34·10⁻⁶, 3.98·10⁻⁶, 6.47·10⁻⁶ and 9.79·10⁻⁶.

Linear in the perturbation — at 10⁻¹⁰ the same column reads 2.27·10⁻¹⁰ through 9.79·10⁻¹⁰, the identical numbers four decades down — and growing with the block by the same factor of about four that the unperturbed column grows by.

So the blocked form has exactly the sensitivity the unblocked form has, per stored number, and it has more stored numbers. Q = I − YTYᵀ is orthogonal when $T$ satisfies $T + T^{\mathsf T} = T^{\mathsf T}Y^{\mathsf T}YT$, which is $r(r+1)/2$ equations in $r(r+1)/2$ unknowns, and the recurrence is a way of solving them rather than a structure that makes them hold. Nothing about the form is a reflection any more: I − YTYᵀ with a perturbed T is a perfectly ordinary matrix with no invariance to fall back on.

That is the difference this essay is about, and it does not show up in the flat line. Both forms produce orthogonality at rounding, for two different reasons: the unblocked one because one number is computed accurately, the blocked one because a hundred and thirty-six are.

One reading of the sensitivity deserves care, because it is easy to take the wrong one. The blocked form is not more sensitive than the unblocked one in any per-number sense: perturb one β by 10⁻⁶ in the unblocked algorithm and the departure is about 4·10⁻⁶, and perturb all 136 entries of T by 10⁻⁶ and it is 9.8·10⁻⁶. So 136 numbers perturbed together cost about twice what one number perturbed alone costs, which is the accumulation of independent errors and not an amplification.

What has changed is not sensitivity. It is that the unblocked form has one quantity to get right and a structural argument that everything else is free, and the blocked form has a hundred and thirty-six quantities and an accuracy claim about a recurrence. Those are different things to have to believe, and the measurement’s job is to say what believing the second is worth: 3.9·10⁻¹⁵ rather than 7.8·10⁻¹⁶.

Why the exchange is worth making anyway #

A block of r reflectors applied to the trailing n − k columns costs the same arithmetic either way — about 4mr(n−k) operations — and moves a different amount of memory. Unblocked, each of the r reflections reads and writes the whole trailing block, so the data crosses the cache boundary r times. Blocked, $YTY^{\mathsf T}$ is formed once and applied with two matrix–matrix products, and the trailing block is read and written once.

On any machine where a multiply costs less than a memory reference — which is every machine built since about 1985 — that is the whole of the difference between a QR that runs at a few per cent of peak and one that runs at sixty.

What the measurement adds is the price in the other currency, stated exactly: a factor of five in ‖QᵀQ − I‖, from 7.8·10⁻¹⁶ to 3.9·10⁻¹⁵, both of which are rounding. Nobody would decline the exchange. The reason to have the number is that it is the kind of quantity that stops being negligible somewhere, and knowing it is a factor of five at r = 16 in double precision says roughly where to start worrying — at a block whose square root of size times the working precision’s rounding reaches whatever the application needs.

The same shape one level up #

A stable block is not a stable basis found the other blocked orthogonalisation behaving very differently: block Gram–Schmidt orthogonalises twice over, between blocks and inside them, and the second pass cannot undo what the first lost. The contrast is worth drawing because both are “the same algorithm, blocked”.

Block Gram–Schmidt blocks an algorithm whose orthogonality was never structural, so blocking it changes where the loss happens without changing that there is one. Blocked Householder blocks an algorithm whose orthogonality was structural, and the blocking replaces the structure with arithmetic — which is a real change in kind, and which costs a factor of five because the arithmetic happens to be accurate.

The general lesson is the one this fleet keeps arriving at from different directions: an algebraic identity between two expressions is not a numerical identity between two computations, and the place to look for the difference is wherever one of them computes something the other one knew. The road that squares the problem is the same statement about the normal equations, two Gram–Schmidts about the order two subtractions are done in, and orthogonal is a number is where the quantity all three are measured in was established.

That last figure is the check that the factor is about the recurrence and not about double precision. If the block’s extra departure came from anything that scaled differently with the working precision — a cancellation, a sum whose length matters — the ratio would move when the precision does. It does not, which is what a straightforward accumulation of independent rounding errors looks like and is the least interesting explanation available, which here is the right one.

What a code could check #

The measurement suggests a cheap self-test that no library ships, and it is worth writing down because the arithmetic is favourable.

$T$ is orthogonality-producing exactly when $T + T^{\mathsf T} = T^{\mathsf T}(Y^{\mathsf T}Y)T$. Both sides are $r$-by-$r$, $Y^{\mathsf T}Y$ costs $r^2m$ operations, and the two triangular products cost $r^3$ — so checking the identity costs $O(r^2m + r^3)$ against the block’s own $O(mr(n-k))$ application. At r = 16 and n − k in the hundreds that is a few per cent, and it is a direct check on the quantity the guarantee rests on rather than on the answer.

It would catch exactly the failure this essay injects. The measured departure of ‖QᵀQ − I‖ is about four times the relative error in T, so a residual of 10⁻⁶ in the identity predicts an orthogonality of 4·10⁻⁶, and a library that computed it could report the number instead of resting on the theorem.

Nothing here argues that it should be on by default — the recurrence is accurate, the flat line says so, and a check that always passes is a check nobody runs. What it argues is that the check exists and is affordable, which is not obviously true of most numerical guarantees and is worth knowing about this one.

What must fail for any of this to be wrong #

Four claims. That the triangle at r = 16 holds 136 numbers, which is arithmetic and is required so that the count in the caption is the count the recurrence actually produces. That every block size is orthogonal to rounding. That the level does not grow with the block by more than a factor of forty — it grows by five. And that a relative perturbation of the triangle costs orthogonality, with a refusal fed the claim that the blocked form’s orthogonality is structural in the way the reflection’s is.

The compact WY assembly and the unblocked reflection are built from one reflector routine, so the r = 1 column of every table is the unblocked algorithm rather than a re-implementation of it. That is what makes “a factor of five” a comparison rather than a coincidence between two programs.

The trade, priced in both currencies #

Putting the two costs side by side is worth doing once explicitly, because the arithmetic is favourable enough that it is easy to state and easy to get backwards.

The blocked form does the same number of floating-point operations as the unblocked one, to leading order, plus the $O(r^2m)$ of forming $T$. What it changes is how those operations are arranged: two matrix–matrix products where the unblocked form does r matrix–vector ones, and the trailing block read and written once instead of r times.

On a machine where a multiply is free relative to a memory reference — which has been every machine since about 1985 — that rearrangement is the difference between a few per cent of peak and sixty. The order the entries are touched in is where this field sets that argument out for a factorisation whose arithmetic is even simpler.

Against sixty per cent of peak, ‖QᵀQ − I‖ rising from 7.8·10⁻¹⁶ to 3.9·10⁻¹⁵ is not a consideration. Both are rounding, the difference between them is a factor of five, and the quantity any downstream computation actually sees is bounded by the larger of the two times whatever that computation’s own condition number is.

So the exchange is not close and nothing here suggests otherwise. What the measurement supplies is the number — so that a future representation trading more of the guarantee for more of the speed has something to be compared against, and so that a reader who was told the blocked form is “the same algorithm” knows in what sense it is not.

What this does not settle #

One shape — 24 rows, 16 columns — one condition number, and six matrices at each block size. The r = 16 block is the whole factorisation rather than one block of several, so nothing here measures how the departures of successive blocks accumulate, which is the quantity a real factorisation of a large matrix would be governed by.

The perturbation is injected into T after the recurrence has run, uniformly and independently across entries. A rounding error in the recurrence is neither uniform nor independent — it is correlated down each column by the triangular solve — so the sensitivity measured is to a perturbation of the wrong shape, chosen because it is the one that isolates T from everything else.

The forward accumulation is never compared against the storage form’s own application. Real codes never form Q; they apply it, and I − YTYᵀ applied to a matrix is a different sequence of operations from I − YTYᵀ formed and multiplied. The orthogonality measured here is of the formed matrix.

And no timing appears anywhere in it. The memory-traffic argument above is standard and is not measured, which is a deliberate limitation — a timed test fails on a busy machine, and operations are the honest currency here rather than seconds.

Still open: accumulation across blocks, and the triangle’s own conditioning #

Several blocks. The r = 16 measurement factorises sixteen columns as one block. A 512-column matrix in blocks of 64 has eight blocks whose departures compose, and whether they add, root-sum-square, or cancel is what decides whether the factor of five is a constant or the start of a trend.

The triangular solve’s conditioning. The recurrence solves against T[1:i−1, 1:i−1], whose diagonal is the βs. Nothing makes those comparable to each other — a reflector whose vᵀv is small has a large β — so a block containing one nearly-null column could give a badly conditioned triangle. Whether such a block arises in practice, and what it does to the factor of five, is the measurement that would say whether the flat line has a worst case.

The storage form’s own application. Applying I − YTYᵀ without forming it changes the rounding, and it is what a code does. Whether the applied form’s departure is the formed form’s, larger or smaller, is one experiment and is the one a library would want.

And the rotation’s block. Givens rotations have no compact form, which is part of why they lost to reflections for dense factorisations and part of why they are still used where the matrix has structure. Whether the argument above — structure replaced by arithmetic when the operation is blocked — has an analogue for a sequence of rotations is the question that would generalise this essay past one algorithm.