A reduction that changes the order

The first essay in this field reorders an elimination and gets the identical factorisation, bit for bit. That is a comfortable result and it is not the general case.

This essay reorders a reduction — a QR factorisation of a matrix with far more rows than columns — and gets a genuinely different sequence of roundings, a genuinely different R, and the same backward error. Which is the more interesting outcome, and the reason it holds is a property this site has been measuring since its foundation phase.

The shape of the problem #

A tall-skinny matrix — 512 rows, 12 columns — is the shape that comes out of a least-squares fit with many observations, out of a block Krylov method, and out of the panel factorisation the previous essay left as the part that does not block.

A Householder sweep on it is a loop over the twelve columns. Each step computes a reflector from the column below the diagonal, applies it to the remaining columns, moves on. The arithmetic is 2mn², which at these sizes is 147,456 operations — a small number. The traffic is not small: each of the twelve steps reads the whole trailing part of the matrix, so the matrix crosses the memory twelve times.

That is a level-2 operation in the vocabulary of the previous essay, and a matrix with 512 rows is exactly the case where the reads dominate. What is wanted is a factorisation that touches the matrix once.

The tree #

Split the rows into p blocks. Factor each block independently — every leaf touches only its own rows and nothing else, so the whole matrix is read exactly once and the p factorisations do not communicate at all. That leaves p small R factors, each 12×12.

Stack them in pairs and factor those. Then stack the results in pairs and factor again. After log₂p levels there is one R left, and it is the R of the original matrix.

The traffic is m words for the leaves plus a negligible amount for the tree — the internal nodes work on 24×12 matrices — against 12m for the sweep. And the leaves are independent, which is why this is the standard algorithm on any machine with more than one processor: the leaves need no communication whatever, and the tree needs log₂p rounds of it.

Why it is the same R, and why that is not obvious #

Every step of this is an orthogonal transformation, and the composition of orthogonal transformations is orthogonal. So the whole reduction is Qᵀ A = R for some orthogonal Q, and the R that comes out satisfies RᵀR = AᵀA — which pins it down to a choice of signs.

That is a statement about the algebra and it is exactly true. The question this site asks is what happens in the arithmetic, and the answer is not a corollary. The tree computes twelve reflectors per leaf, then more at each internal node, in an order that has nothing to do with the sweep’s. Different numbers are subtracted from different numbers at every step. There is no reason for the answers to agree bit for bit and they do not.

What is measured instead is ‖AᵀA − RᵀR‖/‖AᵀA‖, which is the strongest statement available without forming Q — and Q is deliberately not formed, because accumulating it is where the method’s real implementation difficulty lives and eliding that would be dishonest.

The tree returns 1.65·10⁻¹⁵. The sweep returns 9.95·10⁻¹⁵. The tree is, on this matrix, slightly better.

And the error does not grow with the depth #

This is the claim that makes the method a method rather than a curiosity, and it is the one a reader should want checked.

A deeper tree is more levels of rounding. If the backward error grew with the depth then splitting the matrix further — which is what one does to move it less, and what one must do to use more processors — would be paid for in accuracy, and the method would have a size beyond which it stopped being usable.

Measured on the same matrix at four depths:

It does not grow. It wanders inside a factor of three and if anything falls, and the assertion requires the spread across the four depths to be under twenty — which is loose enough not to be measuring noise and tight enough to fail on any real growth.

The reason is the same one as before, applied to the tree rather than to the algorithm: an orthogonal transformation does not amplify. ‖Qx‖ = ‖x‖ exactly, for every x, so whatever error is present when a level begins is neither magnified nor attenuated by that level — it is carried. Errors that are carried rather than amplified accumulate like the number of levels at worst, and log₂16 = 4 levels is not a number that shows up against 10⁻¹⁵.

What classical Gram–Schmidt has that neither of them does #

The bottom of the figure is this site’s oldest measurement, and it is here because it is the control.

Two Gram–Schmidts establishes that the classical process loses orthogonality like κ² and the modified one like κ, while Householder stays flat at 1.8·10⁻¹⁵ across eleven decades. On this matrix at κ = 7,151 that is 4.6·10⁻¹⁰ classical, 1.5·10⁻¹² modified, 2.0·10⁻¹⁴ Householder.

Classical Gram–Schmidt’s steps are projections, not reflections. A projection is not orthogonal: it has a null space, ‖Px‖ ≤ ‖x‖ with equality only for vectors already in its range, and subtracting one from a vector is precisely the cancellation this site’s arithmetic field is about. When the vector being orthogonalised is nearly in the span of the previous ones — which is what an ill-conditioned matrix means — the subtraction destroys the part that was wanted, and what is left is the rounding error.

So the property that makes Householder safe is exactly the property that survives being reassociated, and Gram–Schmidt has neither. That is not two separate facts about two algorithms. It is one fact about orthogonality, seen from two directions, and the tree is the demonstration that it is the right fact: reassociation is safe when the pieces do not amplify, and unsafe when they do.

The site said in its foundation phase that a reflection cannot stop being one. This is what that buys, four phases later, in a field that did not exist then.

What the tree costs that the sweep does not #

Three things, and the essay would be dishonest without them.

Q is harder. The sweep produces Q as an accumulation of reflectors in a standard order and applying it is straightforward. The tree’s Q is a product of the leaf Qs and the tree’s internal Qs, in a structure that has to be stored and traversed. Every real implementation stores it implicitly and there is genuine engineering in doing so. Nothing here does it, and the measurement above is chosen to be the one that does not need it.

The leaves have to be tall enough. A block with fewer rows than columns is rank-deficient by construction, its R is not the R of anything, and the tree means nothing. The generator refuses a configuration where m < np rather than clamping the number of leaves — the same choice hessenberg-reduction makes at n = 8 and jacobi-convergence at n ≥ 8 — because a clamp would draw a four-leaf tree where the essay asked for eight and say nothing.

And the arithmetic is more. The tree factors p blocks and then log₂p levels of stacked pairs, so it performs strictly more operations than the sweep — about 2mn² for the leaves plus O(pn³) for the tree. On a tall matrix the second term is small, and the point of the method is that it is cheaper anyway because operations are not what is being paid for. Which is the field’s whole argument, appearing for the third time in three essays and in the strongest form: here the reordering that moves less data also does more arithmetic, and is still the one to use.

Where the name comes from #

The literature calls this communication-avoiding QR, and this field has already been careful about that phrase once: the previous essay measures blocking’s traffic exponent at 3.182 against an unblocked 3.185 and concludes that blocking buys a constant rather than a class.

Here the phrase is better earned. The sweep reads the matrix Θ(n) times and the tree reads it once, so the ratio grows with the number of columns rather than sitting at 2.1 — and in the parallel setting, which is where the algorithm was developed, the count being reduced is the number of messages, from Θ(n log p) to Θ(log p), which genuinely is a change of class in the parameter that matters.

That distinction is worth keeping because the two essays sit next to each other and use the same words about different-sized effects. A reordering that improves a constant and a reordering that improves an exponent are both worth doing and are not the same kind of result, and this collection’s asymptotic verdict exists because they are routinely described in the same sentence.

Two reassociations, side by side #

This field now contains both possible outcomes of reordering a factorisation, on adjacent pages, and the difference between them is worth stating plainly.

Blocked LU reorders an elimination and returns a bit-identical factorisation — 4.487946226420872 ·10⁻¹⁶ from both orderings, every digit. The reassociation is invisible.

The reduction tree reorders a QR and returns a different R, a different sequence of roundings and a backward error of the same size — 1.65·10⁻¹⁵ against 9.95·10⁻¹⁵. The reassociation is visible and harmless.

Those are different results and the reason for the difference is not the algorithms’ subject matter. Blocked LU defers when an accumulation runs without changing which products go into it or in what order, so every floating-point operation performed is identical to the one the unblocked version performs — it merely happens later. The tree does not defer anything; it computes an entirely different set of reflectors from entirely different data, and the only thing shared with the sweep is what the result must satisfy.

So one is a scheduling change and the other is a different algorithm with the same specification, and the safety of the second rests on a property rather than on an equivalence. Which is the more interesting kind of safety, because it survives changes the first would not: any reduction tree, at any depth, over any partition of the rows, on any number of processors.

And what would have happened with a projection #

The cleanest way to see that the property is doing the work is to imagine the tree built out of Gram–Schmidt instead.

Nothing about the tree’s structure would change. The leaves would still be independent, the matrix would still be read once, the internal nodes would still stack pairs of R factors. Every claim about data movement would go through unaltered.

What would change is that the composition of the levels would no longer be an isometry, so the error present when a level began would be amplified by that level rather than carried — and amplified by a factor that grows with how nearly dependent the columns are. A four-level tree would then be four compounding amplifications instead of four carried errors, and the depth would appear in the answer.

The measurement in the figure is the evidence that this does not happen: 3.38, 4.26, 1.65 and 1.48·10⁻¹⁵ across four depths. It is a null result, and it is the whole reason the method is usable.

The panel, closed #

This field’s three essays close a loop that the second one opened.

A blocked LU factors a panel of b columns and then applies the update. The panel factorisation is unblocked by construction — it is a tall-skinny QR-shaped problem in the middle of an elimination — and the previous essay names it as the part that does not block, and as the reason the measured saving is 2.12 rather than the factor of six the count suggests.

The tree is what to do about it. A panel is n rows by b columns with n ≫ b, which is exactly the shape this essay is about, and factoring it as a reduction tree rather than as a sweep removes the bottleneck the blocking left. That is why “communication-avoiding LU” and “communication-avoiding QR” are usually described together: the second is a component of the first.

Nothing here implements that combination, and it would be the natural next measurement — three counters again, on an elimination whose panel is a tree, against one whose panel is a sweep. What is worth stating without it is that the two essays’ subjects are not two examples of one idea. They are one algorithm with the same idea applied at two levels, and the level the second one reaches is the one the first could not.

The one number the sweep has that the tree does not #

Worth stating so the ladder in the figure is not read as a ranking with a winner.

The sweep produces R and a representation of Q that is standard, compact and easy to apply: the reflectors, stored in the space the matrix vacated, in a fixed order. Applying Qᵀ to a vector is a loop over them. That representation is what LAPACK returns and what every downstream routine expects.

The tree produces R and a representation of Q that is a structure — a list of leaf reflectors plus the internal nodes’ — and applying it means traversing that structure in the right order. It is perfectly doable and it is genuine engineering, and it is why the two methods are not interchangeable at an interface even though they compute the same factorisation.

So the figure ranks four algorithms on one axis, backward error, and the two at the top of it differ on an axis the figure does not have. A reader taking away “use the tree” should take away “use the tree when the matrix is tall enough that reading it twelve times is the cost, and be ready to store a Q that is not a list of reflectors”.

What is left #

The parallel measurement, which is the one the algorithm exists for and which this site cannot make. Everything here is a count on a model of one memory; the interesting quantity is messages between processors, and there are none of those here.

Q, and the applications that need it. A least-squares solve needs Qᵀb rather than Q, which the tree can supply without forming Q — but a block Krylov method needs the orthonormal basis itself, and that is where the implicit representation has to be traversed rather than avoided.

And the other reassociation. The maturity phase found that the implicit double shift and the explicitly formed product agree to 2.2·10⁻¹⁵ entry by entry in absolute value while a norm of their difference is about 8, because “essentially unique” in the implicit Q theorem means up to the signs of the columns. The same caveat applies to the R here and is not measured: the tree’s R and the sweep’s may differ by a diagonal of signs, which is why the comparison is made through RᵀR rather than through R.