The factor nobody forms

Eight blocks and sixty-four reflections asked whether a blocked Householder factorisation multiplies the departure from orthogonality that one block carries, and found that it does not: blocks add in quadrature, and fewer, larger blocks end closer to orthogonal than sixty-four single reflections — 1.42·10⁻¹⁴ in blocks of eight against 2.91·10⁻¹⁴ one at a time. Its conclusion was that blocking, adopted for speed, “also buys a factor of two in orthogonality, at every width a machine is likely to choose.”

Its own last section named the gap in that conclusion. Every number in it was a property of $Q$ formed — the blocks, each stored as a pair $(Y, T)$ with $I - YTY^{\mathsf T}$ its factor, multiplied out into an $m \times m$ matrix and then asked how far $Q^{\mathsf T}Q$ is from the identity. “Most uses of a QR factorisation apply $Q$ or $Q^{\mathsf T}$ to vectors through the stored $(Y, T)$ pairs and never form $Q$,” it said, and whether the applied operator’s departure follows the same law, “it is the same product, applied in a different association,” was left open.

It does not follow the same law, and the difference reverses the conclusion.

Two ways to ask whether Q is orthogonal #

The factorisation is the one the earlier essay measured: textbook blocked Householder QR of 96 × 64 matrices with condition number $10^8$, five seeds, at block sizes 1, 2, 4, 8, 16, 32 and 64. Each panel of $b$ columns is factored by $b$ reflectors, the reflectors are assembled into the compact form $I - YTY^{\mathsf T}$, the trailing columns are updated, and the pair $(Y, T)$ is kept.

Formed, the factor is the product of the blocks’ $m \times m$ matrices, and its departure is the largest singular value of $Q^{\mathsf T}Q - I$. The earlier essay used the Frobenius norm of the same matrix, which is larger by roughly the square root of the dimension and follows the same trends; here both forms are measured in the 2-norm so that they can be compared as operators.

Applied, the factor is never an array. $Qx$ is computed by applying the blocks right to left, each as $x - Y\,(T\,(Y^{\mathsf T}x))$, and $Q^{\mathsf T}y$ by applying them left to right with $T^{\mathsf T}$ in place of $T$. Its departure is the largest round-trip error $\|Q^{\mathsf T}(Qx) - x\|$ over unit vectors $x$, found by iterating the map $x \mapsto Q^{\mathsf T}(Qx) - x$ from a random start — for an operator this close to the identity the iteration settles on its worst direction in a few dozen steps. The same stored numbers go into both; they are rounded in a different order.

The same numbers, a different verdict #

Formed, the factor behaves as the earlier essay found. One reflector at a time it departs by 7.5·10⁻¹⁵; in blocks of two, 5.5·10⁻¹⁵; in blocks of four, 4.1·10⁻¹⁵; from eight to sixty-four it sits between 3.6 and 4.2·10⁻¹⁵. Blocking halves it.

Applied, the block size does not matter. The medians are 3.4, 3.1, 3.5, 3.4, 3.2, 4.0 and 3.7·10⁻¹⁵ from one reflector at a time to one block of sixty-four: a spread of a factor of 1.27 with no trend in it, and a seed-to-seed scatter at each size — 2.3 to 4.0·10⁻¹⁵ one reflector at a time, 2.6 to 5.5·10⁻¹⁵ in blocks of sixteen — larger than any difference between the medians. The applied factor departs about thirty units of roundoff from orthogonality whichever way its reflectors were grouped.

The two agree at the large block sizes and part company at the small ones. One reflector at a time, the formed product departs by 67 units of roundoff and the applied operator by 30: forming the product doubles the departure of a factor that, applied, is as good as any other. So the factor of two that blocking bought is real, and it was bought back from forming. A caller who never forms Q never pays it and has nothing to buy back.

Where the formed product’s extra departure comes from #

Followed reflector by reflector, the two forms separate early and keep separating. One reflector at a time, both start near 9·10⁻¹⁶ after the first reflector — the same computation, a single reflection applied twice. The formed product then grows steadily: 3.5·10⁻¹⁵ after eight reflectors, 5.4·10⁻¹⁵ after twenty-four, 7.5·10⁻¹⁵ by the last. The applied operator grows for a few reflectors — 1.3·10⁻¹⁵ after two, 2·10⁻¹⁵ after three — and then stops, varying between 3 and 4.5·10⁻¹⁵ from twenty reflectors to sixty-four with no further trend. The dashed quadrature line, the first reflector’s departure times the square root of the number applied, reaches 7·10⁻¹⁵ at the end, and the formed product follows it closely; the applied operator ends at half of it.

The reason is in what each rounding touches. Forming the product multiplies an $m \times m$ matrix by the next factor’s $m \times m$ matrix, and every entry of the result is an inner product of length 96 that is rounded; the rounding error of each multiplication is spread over every direction of the space, and the next multiplication cannot undo it. Sixty-four multiplications add sixty-four such errors in quadrature, which is the law the earlier essay found, and the formed factor carries all of them in every direction.

Applying a reflector to a vector rounds only in proportion to how much the reflector moves the vector. A Householder reflector moves $x$ by twice its component along $v$, and the rounding of that update is proportional to the size of the component. A unit vector cannot have a large component along many of sixty-four nearly orthogonal directions at once, so the worst vector for the applied operator is moved substantially by only a few reflectors and sees only their roundings. That is the saturation: once there are more reflectors than a unit vector can be substantially moved by, adding more adds nothing to the worst case.

The argument can be checked on the worst vectors themselves. Taking the direction the iteration settles on for each of the five one-reflector-at-a-time factorisations, and recording its component along each reflector’s direction as it arrives there, gives the same shape every time: one component between 0.69 and 0.93, a handful between 0.2 and 0.5, and the rest small — only one to six of the sixty-four reflectors see a component above 0.3. A random unit vector is spread more evenly, with none or two above 0.3 and eighteen to twenty-two above 0.1. The worst direction is, in effect, the direction of one reflector, chosen so that its own rounding is as large as a rounding can be and the others barely touch it. Adding more reflectors adds more candidates for that role; it does not add their errors together.

Turn the dial to blocks of eight or thirty-two and the same picture appears compressed. Both forms are already near their final size after the first block — eight or thirty-two reflectors in — and the applied operator ends where it began, 3.4 against 3.1·10⁻¹⁵ over eight blocks of eight. The formed product grows a little over its eight blocks, 3.1 to 3.9·10⁻¹⁵, because there are only eight multiplications to add instead of sixty-four.

The earlier derivation, with one term added #

The earlier essay’s two-line account of the quadrature law can be extended to say exactly where the forms differ. It wrote each computed block as an exact orthogonal factor times $I + E_j$ and found the product’s departure to be a sum of the $E_j + E_j^{\mathsf T}$, each rotated by the exact factors after it. Forming the product adds a second kind of term the derivation left implicit: the rounding of each matrix multiplication, a perturbation $F_j$ of the running product that is as large as $u$ times the product’s norm times a factor for the inner products’ length, and that points in no particular direction. There are as many $F_j$ as multiplications, and they add in quadrature with the $E_j$. One reflector at a time there are sixty-four of them; in blocks of sixteen, four.

Applying the factor adds the rounding of each application instead, and that rounding is not a matrix perturbation of fixed size. It is proportional, block by block, to how much the block moves the particular vector in hand — zero for a vector the block leaves alone. The $E_j$ are the same stored numbers in both forms and are small in both; what separates the forms is whether the arithmetic that combines them is charged once per factor to every direction, or once per factor to the directions that factor acts on. Counted that way, blocking reduces the formed product’s charge by reducing the number of factors and cannot reduce the applied operator’s, which was never charged per factor to begin with.

Every run #

Run by run the same division holds. The five one-reflector-at-a-time runs lie well below the diagonal: formed between 6.7 and 8.7·10⁻¹⁵, applied between 2.3 and 4.0·10⁻¹⁵. The thirty blocked runs cluster near the diagonal, and seven of the thirty-five lie above it — blocked runs in which the applied operator’s worst direction happened to depart a little further than the formed matrix’s. The largest applied departure in the whole set, 5.6·10⁻¹⁵, is a blocked run; the largest formed departure, 8.7·10⁻¹⁵, is a single-reflector one. No blocked arrangement is systematically better or worse when applied; forming is systematically worse one reflector at a time.

The typical vector #

A caller applies Q to particular vectors, not to the worst one, and for a random unit vector the round trip is about two and a half times more accurate than the worst case. Formed, the typical error again falls with blocking, from 3.0·10⁻¹⁵ one reflector at a time to 1.3·10⁻¹⁵ from sixteen on. Applied, it is 1.4·10⁻¹⁵ one reflector at a time and rises gently with the block size, to 1.6·10⁻¹⁵ at sixty-four: a large block does a little more arithmetic per application than a reflector does, because $T$ is a dense triangle and $Y^{\mathsf T}x$ and $T\,(\cdot)$ each round. The two curves cross at blocks of sixteen.

So for the typical vector the verdict is not merely “blocking is neutral” but, if anything, the reverse of the formed one: the unblocked application is the most accurate, by a margin of a sixth that is smaller than the seed-to-seed scatter and should not be relied on. What can be relied on is that no block size makes the applied factor meaningfully less orthogonal, which is the property a caller needs.

The ratio of the two worst cases says the same thing as a single curve: 2.2 one reflector at a time, 1.8 at two, between 0.9 and 1.2 from four upward. Forming costs a factor of two when the product has the most factors and nothing measurable when it has few.

What the earlier conclusion should have said #

The earlier measurement is not wrong. The formed factor is exactly what it measured, and anyone who does form Q — to write it out, to take its columns as an explicit orthonormal basis, to compare two factorisations entry by entry — gets that departure, and blocking does halve it for them. What does not follow is the reading that blocking improves the factorisation’s orthogonality. The factorisation is the set of stored reflectors; its orthogonality, as a caller experiences it, is the applied operator’s, and that is the same at every block size.

The practical rule is therefore about forming rather than about blocking. If Q must be formed, form it from large blocks, because the formed product accumulates one rounding error per factor multiplied in, and fewer factors is fewer errors. If Q is only applied, the block size is a question of speed alone — the one a block size is a property of the machine answered from the fast memory — and has no accuracy term. A triangle where the scalar was measured what assembling one block’s triangle costs, a factor of five on one block formed; applied, even that block departs no further than a single reflection’s round trip does after twenty reflections.

This is the same distinction numerical linear algebra keeps teaching in other clothes. Orthogonal is a number turned orthogonality into ‖QᵀQ − I‖; the number is only meaningful for the object it is computed on, and “Q” names two different computations here. A reflection cannot stop being one found that a single computed reflector is orthogonal to working accuracy by construction; the applied factor inherits that, reflector by reflector, and the formed one dilutes it in sixty-four matrix products.

The use that matters most: Qᵀb #

The commonest application of the factor is the least-squares solve, which computes $Q^{\mathsf T}b$ through the stored reflectors, reads off the first $n$ entries and back-substitutes with $R$. That is the applied operator on one vector, and its accuracy is the applied departure, not the formed one. The residual the appended block cannot remove measured the analogous question for Gram–Schmidt, where $Q$ is formed by construction and $Q^{\mathsf T}b$ inherits every bit of its loss of orthogonality; the Householder factor applied to $b$ has no formed $Q$ to inherit from, and on these matrices its round trip is within thirty units of roundoff whatever the block size. Two Gram–Schmidts is the reminder of how large the difference between methods can be when the factor is formed from the data; here the difference is only in the order of the arithmetic, and it is a factor of two at most.

It also changes what one number that has to be right was protecting. That essay found the scalar β in each reflector deciding its orthogonality, because an error in β is an error in every application of that reflector. Applied, a bad β is felt in the one direction the reflector moves, and the saturation above says a caller’s worst vector will usually find the one reflector with the worst β. Formed, a bad β is folded into a dense matrix and felt everywhere. The scalar matters in both forms; which vectors pay for it differs.

What five seeds and one shape do not show #

One matrix shape, 96 × 64, and one condition number, which the earlier essay found makes no difference. The worst direction of the applied operator is found by iteration on a map that is linear only up to rounding, so it is an estimate from below; it agrees with the formed matrix’s 2-norm where both are computed on the same operator and settles within a few per cent over the last twenty iterations, but it is not a certificate. Seven of thirty-five runs have the applied departure above the formed one, which is within the scatter and means the statement is “no worse”, not “better”. And a library applies $Q$ to blocks of vectors at once, with its own ordering of the inner products; the saturation argument does not depend on the ordering, but the constants do.

Still open: blocks of blocks, and a matrix instead of a vector #

Recursive blocking. Libraries factor a panel in blocks and merge the triangles, so a factor is a block of blocks. Formed, that adds a level of products and, on this evidence, a level of accumulation; applied, the prediction with a sign is that it adds nothing, because the applied operator’s worst direction is still moved by only a few reflectors however they are grouped.

Applying to many vectors. $Q^{\mathsf T}B$ for a matrix $B$ with many columns is how a least-squares solve with many right-hand sides uses the factor. Each column is a vector and should see the saturated departure; the whole block’s departure, measured as the largest over the columns, should then grow only as the chance that some column is close to a worst direction — slowly in the number of columns.

Where the saturation begins. Twenty reflectors on 96-long vectors is where the growth stopped here. The argument above predicts the saturation count depends on the vector length and not on the number of reflectors, so that on 384-long vectors the growth would continue longer before stopping. That is one change of dimension in the same experiment.