Eight blocks and sixty-four reflections

A triangle where the scalar was took apart the representation every level-3 QR factorisation uses. A block of $r$ Householder reflectors is applied not one at a time but as a single matrix, $Q = I - Y T Y^{\mathsf T}$, with the reflectors as the columns of $Y$ and an upper triangular $T$ computed by a recurrence whose inputs are its own earlier columns. One reflection’s orthogonality rests on one computed scalar; a block of sixteen rests on 136 computed numbers. Measured on a 24 × 16 matrix, the block departed from orthogonality by 3.9·10⁻¹⁵ against the single reflector’s 7.8·10⁻¹⁶: a factor of five.

That measurement was of one block. A factorisation of a 512-column matrix in blocks of 64 is eight blocks in sequence, each built from a panel the earlier blocks have already transformed, and its orthogonal factor is their product. The essay closed on the question that decides whether the factor of five matters: do the blocks’ departures add, add in quadrature, or cancel? If they add, a factorisation of many blocks is many times worse than the measurement suggests. If they cancel, it is no worse than one.

They add in quadrature — and the consequence is that a blocked factorisation ends closer to orthogonal than an unblocked one.

A factorisation of many blocks #

The measurement is the textbook blocked Householder QR. Factor a panel of $b$ columns by ordinary reflectors, one at a time, applied within the panel. Assemble the panel’s reflectors into $(Y, T)$ by the recurrence. Apply $I - Y T^{\mathsf T} Y^{\mathsf T}$ to the columns to the right, and accumulate $Q \leftarrow Q\,(I - Y T Y^{\mathsf T})$. Record $\lVert Q^{\mathsf T}Q - I\rVert$ after every block, and each block’s own $\lVert (I - YTY^{\mathsf T})^{\mathsf T}(I - YTY^{\mathsf T}) - I\rVert$ beside it. With $b = 1$ the same procedure is unblocked Householder, one reflection per factor.

The matrices are 96 × 64 with condition number $10^8$, five of them, and every number below is the median over the five.

How eight blocks accumulate #

In blocks of eight the accumulated departure is 5.4·10⁻¹⁵ after the first block, 8.1·10⁻¹⁵ after the second, 1.11·10⁻¹⁴ after the fourth and 1.42·10⁻¹⁴ after the eighth. The blocks’ own departures run from 5.4·10⁻¹⁵ for the first down to 3.2·10⁻¹⁵ for the last, because each later panel is shorter. Summed, they reach 3.29·10⁻¹⁴; their root-sum-square is 1.18·10⁻¹⁴. The accumulated departure tracks the root-sum-square from above, 1.2 times it at the end, and stays well below the sum.

That is the behaviour of independent errors. Each block’s factor is orthogonal up to a small symmetric error of its own, and the product of several such factors is orthogonal up to, to first order, the sum of their errors after each has been rotated by the factors around it. Rotations do not shrink an error and do not align it with the others, so the errors add as random vectors do: their squared sizes add. The fifth of excess over the root-sum-square is what is left of alignment between the blocks, which share a matrix and are not entirely independent.

Taken one reflector at a time, the same factorisation accumulates 64 errors instead of 8. Each reflector’s own departure is smaller — about 2.5·10⁻¹⁵ at the start and 1.5·10⁻¹⁵ at the end — but there are eight times as many, and the product ends at 2.91·10⁻¹⁴, twice the blocked factorisation’s.

Why the square root #

The quadrature law can be written down in two lines, and writing it down says what could break it. Let each block’s computed factor be $\hat H_j = H_j (I + E_j)$, where $H_j$ is exactly orthogonal and $E_j$ is small. Then $\hat H_j^{\mathsf T}\hat H_j - I = E_j + E_j^{\mathsf T}$ to first order, which is the block’s own departure. The product of eight such factors is $H_1 \cdots H_8$ times a correction, and to first order

$$\hat Q^{\mathsf T}\hat Q - I \approx \sum_j G_j^{\mathsf T}\,(E_j + E_j^{\mathsf T})\,G_j, \qquad G_j = H_{j+1}\cdots H_8 ,$$

a sum of the blocks’ own departures, each rotated by the exact factors that come after it. A rotation changes no norm, so each term is as large as its block’s own departure. What decides the size of the sum is whether the terms point the same way. Terms that were produced by different roundings, in different panels, and then rotated by different products have no reason to align, and a sum of $N$ unaligned terms of size $\varepsilon$ has size about $\sqrt{N}\,\varepsilon$, not $N\varepsilon$.

So the law would fail for errors that are systematically aligned — a rounding mode that always biased in one direction, say, or a construction in which every block saw the same panel. Neither is the case here, and the measured factor of 1.2 above the root-sum-square is the small alignment that remains.

Fewer, larger factors #

The comparison generalises across block sizes, and the two effects it balances can be read separately.

One block’s own factor departs by a median of 2.0·10⁻¹⁵ at a block of one, 3.1·10⁻¹⁵ at four, 4.3·10⁻¹⁵ at eight, 5.4·10⁻¹⁵ at sixteen, 9.6·10⁻¹⁵ at thirty-two and 1.4·10⁻¹⁴ at sixty-four, where a single block is the whole factorisation. Quadrupling the block multiplies its departure by one and a half to two and a half. The whole factorisation is the product of $64/b$ of them and departs by 2.91·10⁻¹⁴ unblocked, 1.66·10⁻¹⁴ at four, 1.42·10⁻¹⁴ at eight, 1.24·10⁻¹⁴ at sixteen, 1.27·10⁻¹⁴ at thirty-two and 1.39·10⁻¹⁴ at sixty-four.

If each factor contributes independently, the total goes as the square root of the number of factors times the departure of one. The number of factors falls as $1/b$ and the departure of one rises, very roughly, as $b^{1/2}$, so the product should go as $\sqrt{64/b}\,\cdot\, b^{1/2}$, which is flat — and the measured total is nearly flat above a block of eight, with a shallow minimum at sixteen. Below eight the total rises because the departure of one factor stops falling: a single reflector is not free, and 64 of them in quadrature cost more than eight blocks.

So the factor of five the earlier essay measured is real and is paid per block, and a factorisation pays it fewer times than the unblocked factorisation pays its own smaller price. Blocking was adopted for speed — the essay on the triangle priced it at matrix–matrix products instead of matrix–vector ones, and a block size is a property of the machine showed how the best width is set by the fast memory rather than by the matrix — and on this measurement it also buys a factor of two in orthogonality, at every width a machine is likely to choose.

Signs that do not agree #

There is a smaller version of the same fact one level down, in the quantity being measured. Orthogonal is a number made $\lVert Q^{\mathsf T}Q - I\rVert$ the site’s measure of orthogonality: the entries of $Q^{\mathsf T}Q$ are inner products of computed columns, and the departure is the size of how far they are from zero and one. An inner product with no fixed sign looked at those entries across ten ways of partitioning one computation and found 125 of 1,128 off-diagonal entries taking both signs, while the total moved by only 2.4 per cent: the entries are rounding, and rounding has no preferred sign.

The quadrature law between blocks is that observation at the scale of whole factors. Each block contributes a matrix of small errors whose entries have no preferred sign, and eight such matrices, rotated by different products, overlap in their signs about as often as not. Where two blocks’ errors reinforce in one entry they cancel in another, and the total grows as the square root of the count. The unblocked factorisation is sixty-four such contributions. Counted alone, at the square root, eight times as many terms would make it 2.8 times the blocked total if each were as large as a block’s; each is in fact about half as large, which predicts 1.3. It ends at 2.05 times, because its sixty-four errors are more aligned than eight blocks’ are — its total is 1.7 times their root-sum-square where the blocks’ is 1.2. Why single reflections align more is not isolated here; the likeliest reason is that successive reflections share more of their arithmetic than successive blocks do, and what they share, they round alike.

How many numbers the guarantee rests on #

The earlier essay counted the numbers that have to be right: one scalar per reflection, $b(b+1)/2$ per triangle. Counted over the whole factorisation, the blocked form rests on more of them, not fewer. Sixty-four reflections carry 64 scalars. Eight blocks of eight carry eight triangles of 36, 288 numbers. Four blocks of sixteen carry 544, and one block of sixty-four carries a triangle of 2,080.

And the whole factorisation’s departure falls as those counts rise, from 2.91·10⁻¹⁴ on 64 numbers to 1.24·10⁻¹⁴ on 544. The count of computed numbers was a natural thing to worry about and it is the wrong measure. What matters is how many independent factors the errors are committed in, because each factor’s error is one term in the quadrature sum, and a triangle’s 36 numbers produce one term, not 36.

The same arithmetic runs through the same arithmetic at a different price, which found a blocked and an unblocked elimination performing the same 72,568 operations associated differently. Association is also what decides the rounding: the blocked factorisation groups its operations into eight factors and the unblocked into sixty-four, and the rounding follows the grouping.

Independent of the conditioning #

The totals do not depend on the matrix. Blocks of eight end between 1.37·10⁻¹⁴ and 1.42·10⁻¹⁴ at condition numbers from $10^2$ to $10^{14}$; one reflector at a time, between 2.78·10⁻¹⁴ and 2.98·10⁻¹⁴. That is the property a reflection cannot stop being one established for a single Householder QR — its orthogonality is a property of how the factor is built, not of the matrix it is built from — and it survives blocking whole. A modified Gram–Schmidt factorisation of the same matrices would lose orthogonality in proportion to the condition number, and none of the blocked forms here do.

A triangle with entries of 10²⁶ #

The earlier essay’s other open question was about the triangle itself. The recurrence solves against its own earlier columns, whose diagonal holds the reflectors’ scalars $\beta = 2/v^{\mathsf T}v$, and nothing makes the βs comparable to each other. A reflector whose $v$ is small has a large β. A block containing a nearly dependent column — one that is almost a combination of the columns before it — should produce such a reflector, since what is left of the column after the earlier reflectors have acted is small. Whether that makes $T$ badly conditioned, and the block’s departure worse, was the measurement it asked for.

The triangle does blow up, exactly as predicted. With the fourth column of a well-conditioned matrix set to the third plus ε times a random vector, the largest entry of the first block’s $T$ is 86 at ε = 10⁻², $8.6 \times 10^9$ at 10⁻⁶, $8.6 \times 10^{17}$ at 10⁻¹⁰ and $8.6 \times 10^{25}$ at 10⁻¹⁴ — growing as $1/\varepsilon^2$, the square of the small residual’s reciprocal. The orthogonality does not move. The first block’s own departure is between 4.0·10⁻¹⁵ and 5.4·10⁻¹⁵ at every ε, and the whole factorisation’s between 1.33·10⁻¹⁴ and 1.45·10⁻¹⁴.

The reason is a scale invariance. A Householder reflection $I - \beta v v^{\mathsf T}$ is unchanged if $v$ is multiplied by any $s$ and β divided by $s^2$: the reflection depends only on the direction of $v$. So is the compact form: multiplying a column of $Y$ by $s$ and the corresponding row and column of $T$ by $1/s$ leaves $Y T Y^{\mathsf T}$ exactly the same. A tiny $v$ with a huge β is a unit vector with a moderate β, written in different units. The triangle’s entries are large because its reflector is small, and they are large in exactly the proportion that cancels in the product. What must be right is not the size of $T$’s entries but their relative accuracy, which is what one number that has to be right found for the single scalar: a relative error of $10^{-2}$ in β costs the orthogonality $10^{-2}$, and a huge β computed to working precision costs nothing.

This does depend on the implementation computing β consistently with $v$, as $2/v^{\mathsf T}v$ from the stored $v$. A code that normalised $v$ to have a unit first component, as LAPACK does, would store a different pair — and the scale invariance says it would make no difference to the product, which is the point.

What changes in the account of blocking #

The earlier essay’s summary was that blocking moves the guarantee from a structure to a recurrence: a block of sixteen rests on 136 computed numbers, and they cost a factor of five. Both halves stand. What the factorisation adds is that the cost is a per-factor cost, that factors accumulate in quadrature, and that a blocked factorisation has fewer factors. The trade was presented as speed against a small loss of orthogonality; measured on whole factorisations, it is speed and a small gain.

Block Gram–Schmidt makes an instructive contrast. A stable block is not a stable basis found that orthogonalising inside each block with Householder does not save a block Gram–Schmidt factorisation, because the blocks are orthogonalised against each other by projections whose error grows with the conditioning. The Householder factorisation’s blocks are joined by multiplication of near-orthogonal factors, which is why they accumulate in quadrature rather than in proportion to κ. The block is not what makes a factorisation stable; the way blocks are combined is.

It also sharpens what a code could check, which the earlier essay proposed: after assembling each block, test $\lVert (I - YTY^{\mathsf T})^{\mathsf T}(I - YTY^{\mathsf T}) - I\rVert$ against a few units of roundoff times $b$. A block that fails has a triangle wrong in relative terms, which is the only way it can fail; a block whose triangle is merely huge passes, correctly.

What this does not settle #

One shape of matrix, 96 × 64, and Gaussian-built matrices of prescribed conditioning. The block’s own departure grows with the length of the panel as well as its width, and a tall matrix — 10,000 × 64 — would change the per-factor cost without changing the number of factors. The quadrature law should survive; the location of the shallow minimum over block sizes need not. Nor is the alignment that puts the unblocked factorisation at 1.7 times its root-sum-square, against the blocked one’s 1.2, explained by anything measured: it is the one number in the comparison that a different matrix family could move, and the comparison’s factor of two rests partly on it.

The departures are of the accumulated $Q$ formed explicitly. A code that never forms $Q$, and applies the stored blocks to vectors instead, commits the errors in a different order, and the orthogonality of that implicit $Q$ is not the same number. The earlier essay listed the applied form’s departure as the measurement a library would want, and it is still not made.

Still open: the applied form, and blocks of blocks #

The factor never formed. Most uses of a QR factorisation apply $Q$ or $Q^{\mathsf T}$ to vectors through the stored $(Y, T)$ pairs and never form $Q$. The departure that matters there is how far $Q^{\mathsf T}(Q x)$ is from $x$ for the applied operator. Whether that operator’s departure follows the same quadrature law as the formed matrix’s — it is the same product, applied in a different association — is the measurement that would carry this result into what libraries do.

Recursive blocking. High-performance codes block recursively: a panel is itself factored in blocks, and the triangles are merged. Each level adds a recurrence and multiplies the numbers that must be right, and each level also reduces the count of factors at the level above. Whether the quadrature law makes recursion a further gain in orthogonality, or whether the merged triangles cost more than they save, is a question the same instrument can answer at a few levels.