A tree leaks along its leaves

The worst residual belongs to the route turned a slack bound into a measured map. Every stable least-squares route factorises A before it sees b, so for a fixed matrix its error is, to first order, a linear function of the right-hand side, and with the consistent part held fixed it is a linear map from the residual — the part of b outside A’s range — to the error in the solution. On a 64 × 16 problem that map has one dominant direction for every route, so a residual aimed along it draws exactly $\sqrt{m - n}$ times what a random one does, and the direction is the route’s own: aimed at Householder, the residual leaves block modified Gram–Schmidt far below its own worst, and the other way round.

It closed on the shape a distributed code factorises in. “A tree of Householder factorisations commits its roundings leaf by leaf and then in the combining steps, so its leak has its own direction. Appending b to every leaf would inherit Householder’s accuracy at each leaf; whether the tree’s aimed residual is worse than the sequential one’s, and whether the combining steps or the leaves set its direction, is the block question asked where a distributed code asks it.”

A tree has its own direction, set by how the rows are split, and its aimed residual draws about what the sequential route’s does. The appended column, which looked like the interesting design choice, is not a choice at all.

Three routes and their trees #

The problem is the earlier essay’s: a 64 × 16 matrix with condition number $10^8$, a solution $x$ with entries from 1 to 2.5, and a right-hand side $b = Ax + sz$ with z a unit vector orthogonal to the range and s set so that the relative residual is one per cent. Two placements of the ill-conditioning are measured. In the first it is spread through the columns, the matrix a generic one with singular values falling geometrically. In the second it sits between four blocks of four columns, each block orthonormal and the blocks nearly dependent on each other — the case where the earlier essays found block Gram–Schmidt’s appended block falling behind.

The tree splits the 64 rows into leaves. With two leaves each has 32 rows; with four, 16. Each leaf is factorised by Householder. The leaves’ triangular factors are stacked in pairs and factorised again, and so on up a binary tree to a single 16 × 16 triangular factor at the root — the reduction that a tall, skinny QR on many processors performs, where each processor owns a leaf and only the small triangles travel. The right-hand side goes up the tree one of two ways. Carried, each leaf’s orthogonal factor is applied to its slice of b, and the resulting coefficients travel up beside the triangles and are transformed at every combining step. Appended, b is the seventeenth column of every leaf, factorised with it, and its coefficients and residual norm travel up inside the triangles.

For each route the error map is measured exactly as before: for each of the 48 directions orthogonal to the range in turn, the right-hand side is formed, the route solves, the exact rational least-squares solution is subtracted, and the difference is one column of the map. The map’s singular values and singular vectors are then the route’s constants and directions, scaled to the bound’s $\kappa^2 u \rho$.

A tree’s leak is one direction too #

The trees’ maps have the shape the sequential routes’ had. On the spread placement the first singular value holds 99.9 per cent of the sequential map’s square, 99.9 per cent of the two-leaf tree’s and 100.0 of the four-leaf tree’s, and the next is smaller by orders of magnitude. Between column blocks the shares are 94.7, 99.5 and 99.0 per cent. So the aimed residual of a tree draws $\sqrt{48} \approx 6.9$ times a random one’s error, as every route’s did — 6.93 for both trees on the spread placement, 6.91 and 6.89 between blocks — and the earlier essay’s arithmetic carries over: a tree has a single worst residual, and the map’s top singular value is its constant.

The mechanism is the earlier essay’s too. A route’s error on the residual’s part comes from the computed factorisation’s range being tilted, by an amount of order $\kappa u$, away from A’s exact range; the residual has a component along the tilt, and the triangular solve amplifies it by another κ. A tilt of one matrix by a backward error is almost rank one in its effect, because only the smallest singular direction gets the second κ. A tree tilts the range differently from the sequential factorisation, but it still tilts it once, and the second κ still falls on one direction.

Of Householder’s order #

The figure above is the constants. With the ill-conditioning spread through the columns, the sequential route’s aimed residual draws 0.113 of $\kappa^2 u\rho$; the two-leaf tree’s 0.090 and the four-leaf tree’s 0.130. A tree of two is slightly better and a tree of four slightly worse, both within a fifth. With the ill-conditioning between column blocks, sequential Householder draws only 0.014 — this placement is kind to it, since its reflections process the blocks’ directions in turn — and the trees draw 0.042 and 0.031, two to three times as much.

A factor of three between blocks is not a factor any user would notice on one solve; at a residual of one per cent on a matrix with $\kappa = 10^8$ the bound $\kappa^2 u\rho$ is about $10^{-2}$, and the aimed errors are $1.6 \cdot 10^{-4}$ for the sequential route and $4.7 \cdot 10^{-4}$ for the two-leaf tree — about four correct digits against three and a half. It matters to the extent the direction does: for most right-hand sides neither route is near its worst. Neither number is large beside the bound. Every constant here is below a seventh of $\kappa^2 u\rho$, so the earlier essays’ finding that the bound is loose by a factor of seven to ninety on these matrices survives the change of shape. What changes between the placements is how much the sequential order’s particular sequence of reflections happens to suit the matrix. Between blocks, the sequential route meets each block’s four orthonormal columns together and keeps them together; a tree splits the rows, so every leaf sees every block’s columns as sixteen generic vectors and has no such luck.

Each factorisation leaks its own way #

The trees’ aimed directions are not the sequential route’s. With the ill-conditioning spread, the two-leaf tree’s lies at a cosine of 0.58 from it and the four-leaf tree’s at 0.40; between blocks, 0.18 and 0.13 — nearly orthogonal in a 48-dimensional space where a random pair of directions would sit at a cosine of about 0.14. And the two trees disagree with each other: 0.75 on the spread placement, 0.13 between blocks. Change the number of leaves and the worst residual moves to a different direction.

That answers half of the earlier essay’s question. Backward stability is a statement about the size of the perturbation a route commits, not about its direction, and every route measured here — Householder, the block methods of the earlier essays, two trees — is backward stable with its own perturbation. A residual that is worst for one of them is a random residual for another.

The leaves decide where it leaks #

The other half — whether the leaves or the combining steps set the direction — shows in where the aimed residual lives. With the ill-conditioning spread, sequential Householder’s aimed residual has three quarters of its square on the first sixteen rows: the first reflections, which act on every row but are chosen from the first columns’ full length, commit the rounding the map reads most. The two-leaf tree’s aimed residual is split 55 to 45 between its two leaves, the halves of the rows. The four-leaf tree’s puts 38 and 39 per cent on two of its four leaves and 10 and 13 on the others. Between blocks the four-leaf tree’s spreads almost evenly, 22 to 27 per cent a quarter.

The four-leaf tree is the sharper test of what sets this, because its leaves are square. A 16 × 16 leaf has a full orthogonal factor and leaves nothing of its slice of b outside its range; every rounding that reaches the residual’s part of the error is committed either in that orthogonal factor or in the combining steps above it, which work on stacked triangles with no rows of the original matrix in them. And still its aimed residual sits mostly on two of its four leaves, and still it lies far from the two-leaf tree’s direction. The combining steps’ rounding reaches the residual only through the leaves’ orthogonal factors, which map it back onto the leaves’ rows, so whatever commits it, the split of the rows shapes where it lands. A different split is a different direction; that much the measurement separates. Whether the leaves’ own factorisations or the combining steps commit the larger share of the tilt, it does not.

Two Gram–Schmidts found one argument changed in an inner product moving an orthogonalisation’s accuracy by orders of magnitude. The tree is the opposite case: the order of the arithmetic is reorganised entirely, and the size of the error barely moves while its direction moves completely.

Appending b at the leaves is the same computation #

The design choice the earlier essay expected to matter does not. With the ill-conditioning between blocks, the two-leaf tree carrying Qᵀb draws 0.04239 of the bound and with b appended 0.04239; the four-leaf trees 0.03064 and 0.03064; on the spread placement 0.09020 against 0.09020 and 0.12988 against 0.12988. The aimed directions agree to a cosine within $2 \cdot 10^{-10}$ of one on every pair.

The reason is that for Householder the two are one computation in exact arithmetic and nearly one in floating point. Appending b as a column applies the leaf’s reflectors to b as they are formed; carrying it applies the same reflectors, accumulated into an explicit orthogonal factor, to b afterwards. The difference between the two is a rounding of order u in the coefficients — the consistent part’s error, of order $\kappa u$ in the solution — and the residual’s part of the error, which the map measures, comes from the tilt of the triangular factor, which both share. What the appended block inherits found the opposite for block Gram–Schmidt, where appending b rescues a route whose Q is not orthogonal; Householder’s Q is orthogonal to working accuracy at every leaf, and there is nothing to rescue.

Orthogonality was never the question #

It is tempting to rank the routes by how orthogonal their computed factors are, and for least squares on the residual’s part that ranking says nothing. Orthogonal is a number made orthogonality a measured quantity, $\|Q^{\mathsf T}Q - I\|$, and a reflection cannot stop being one explained why Householder holds it at the rounding level whatever the matrix: each reflector is built from a unit vector, and a rounded unit vector names a slightly different reflection, still exactly orthogonal. Every leaf here and every combining step is a Householder factorisation, so every orthogonal factor in the tree is orthogonal to working accuracy, and eight blocks and sixty-four reflections found that composing many such factors does not multiply their departure. The tree’s orthogonality is as good as the sequential route’s.

Its errors on the residual differ anyway — by a factor of three between blocks, in direction completely — because what the map measures is not how orthogonal the computed Q is but which slightly different matrix the computed Q and R exactly factorise. Two backward-stable routes factorise two different nearby matrices, and the residual’s error reads the difference between their ranges. A stable block is not a stable basis found the converse failure in block Gram–Schmidt, where each block is orthogonal and the basis is not; here every block and every basis is orthogonal, and the routes still disagree about which right-hand side is worst.

What a distributed code can take from this #

A tree is the shape a parallel least-squares solve takes, and the measurement says three things about it. Its accuracy on the residual’s part is of the sequential factorisation’s order — a little better or a little worse, by a factor of up to three on these matrices, depending on whether the sequential order happens to suit the matrix. Its worst right-hand side is a different one from the sequential factorisation’s, so a test suite that checks a parallel solver against a sequential one on a single residual direction may see agreement on a direction neither is weak on, or a disagreement that says nothing about either’s worst case. And the way b enters the tree is free to be chosen for communication rather than accuracy: appending it to every leaf costs one more column in every triangle that travels, carrying it costs a separate vector beside the triangle, and the error is the same to five figures either way.

The residual the appended block cannot remove found every stable route climbing on the bound’s second term at about a fiftieth of it; the earlier essay put the aimed constant at a ninth for Householder. A tree does not move that number by more than its factor of three. What it moves is the direction — and with it, which right-hand sides a given machine’s answer is worst on.

What two placements do not show #

One 64 × 16 matrix shape, one condition number, one relative residual, two placements of the ill-conditioning, and trees of two and four leaves. A tree with many more leaves — a hundred processors with a few rows each — has leaves wider than they are tall, which pass nearly all of the arithmetic to the combining steps; the measurement here cannot say whether the direction then still follows the split. Householder is the only leaf factorisation measured; a tree with Cholesky QR or Gram–Schmidt leaves would have leaves whose Q is not orthogonal, and there appending b might matter as it did for the block methods. The relative residual is fixed at one per cent, but the map is linear in the residual, so its constants hold at any residual large enough that the residual’s part of the error outweighs the consistent part’s, which is of order $\kappa u$ and the same for every route; below about $\rho = 1/\kappa$ the routes become indistinguishable again, whatever their directions. The trees are binary and balanced; an unbalanced reduction, as on processors that finish at different times, would weight the leaves differently.

Still open: the size of the tilt, and a tree of many leaves #

The size of the leak. The earlier essay’s first question still stands, and the trees give it a second instrument. The remaining slack is the size of the tilt between the computed range and the exact one along the direction the map reads. The prediction with a sign is that for every route here, the aimed constant equals the computed range’s principal angle with the exact range along the smallest singular direction, times κ over $\kappa^2 u$, to within a factor of two — so that the trees’ factor of three between blocks is a factor of three in tilt, visible without solving anything.

Many leaves. With sixteen leaves of four rows each, every leaf is a 4 × 16 block, wider than it is tall, whose factorisation commits almost nothing and passes the whole of the residual’s rounding to the combining steps. The prediction with a sign is that the aimed direction then spreads over the leaves — its share on each leaf’s rows close to a sixteenth, unlike the four-leaf tree’s 38, 10, 39 and 13 per cent — and that the aimed constant stays within a factor of three of sequential Householder’s, as here.