The worst residual belongs to the route

The residual the appended block cannot remove gave a least-squares right-hand side a component outside the range of $A$ and watched every stable route’s error rise with it. The perturbation theory says why: the error of a backward-stable solve has a term of order $\kappa u$, paid at any residual, and a term of order $\kappa^2 u\rho$, where ρ is the residual relative to the right-hand side, which no route avoids. Householder and block modified Gram–Schmidt with $b$ appended as one more block both climbed on that second term — at about a fiftieth of it. The bound’s constant, that essay said, “on these matrices is generous”.

Generous is not an explanation, and the essay named the one it suspected. Its residual was a fixed combination of the directions orthogonal to the range, $\sin(k+1)$ times the $k$-th of them, chosen without regard to anything. “The worst case for the $\kappa^2$ term is a residual … aligned with how those directions sit in the full space. Aiming it there and sweeping again would say whether the knee moves by the factor of fifty the bound’s slack suggests, and whether the appended block’s threshold of one per cent shrinks with it.”

The aim turns out to be computable exactly, it closes a factor of seven, and the direction it finds belongs to the algorithm rather than to the problem.

The map from residual to error #

Every route here factorises $A$ before it looks at $b$. Householder computes its reflections from $A$ alone; block modified Gram–Schmidt orthogonalises the blocks of $A$ and only then runs the appended right-hand side through the same projections; column modified Gram–Schmidt on $[A\;b]$ treats $b$ as the last column, after every column of $A$ is done. So for a fixed matrix, the factorisation’s rounding is fixed, and to first order the solution’s error is a linear function of the right-hand side.

Hold the consistent part $b_0 = Ax$ fixed and add $s\,z$, with $z$ any unit vector orthogonal to the range. The error is then $E z$ times $s$, plus the consistent problem’s own error, for a fixed matrix $E$ that maps the $m - n$ dimensions of the complement into the $n$ dimensions of the solution. It can be measured one column at a time: take each of the 48 basis vectors of the complement of a 64 × 16 matrix in turn, solve, subtract the exact solution — computed in rational arithmetic from the stored doubles, as every reference in this series is — and divide by $s$. Forty-eight exact solves per route.

The residual that $E$ amplifies most is its top right singular vector, and that is the aimed residual: not a guess at what “weak directions” means, but the direction in the complement that this factorisation of this matrix turns into the most error. A random direction draws, on average, the root-mean-square singular value of $E$. And the prediction can be checked by doing what the measurement predicts — a fresh solve with the residual along the aimed direction, from scratch, against a fresh exact reference. On both placements and all three routes the fresh solve reproduces the map’s top singular value to within five per cent; on the generic matrix Householder’s predicted 0.1127 and measured 0.1127 agree to four figures.

This is the measurement the condition number is an amplifier would recognise, turned on a solver instead of a problem: perturb the input along every direction available, look at how far the output moves, and take the largest ratio. The difference is which output. The problem’s condition number is the largest ratio over all perturbations of $A$ and $b$, and it is what the bound is built from. The aimed constant is the largest ratio over residuals for the rounding one particular algorithm actually committed, on one particular matrix — a number that exists only after the factorisation and is never larger than the bound’s.

Aimed, the constant rises sevenfold #

The figure at the top of the page is a generic 64 × 16 matrix at κ = 10⁸ and a relative residual of one per cent, with the error expressed as a multiple of $\kappa^2 u\rho$, so that the bound is the dashed line at one. Each dot is one of twenty-four random residual directions; each bar is the aimed residual; each ring is the fixed direction the earlier essay used.

For Householder the random directions draw from 8·10⁻⁴ to 0.043 of the bound, median 0.014, and the fixed direction drew 0.011 — an ordinary draw. The aimed residual draws 0.113. For block modified Gram–Schmidt with $b$ appended, the median is 0.0073 and the aimed residual 0.082; for column modified Gram–Schmidt on $[A\;b]$, 0.0031 and 0.033. On all three the aimed residual is eight to eleven times the median random one, and still a factor of nine to thirty below the bound.

The $\kappa^2$ in that bound is the same squaring the road that squares the problem found the normal equations paying on every right-hand side, and here it is paid only through the residual: with ρ = 0 the stable routes lose κu and nothing more. That is the whole case for avoiding $A^{\mathsf T}A$, and the $\kappa^2 u\rho$ term is the part of the squaring that no algorithm can avoid, because it is in the problem’s own sensitivity rather than in the method.

So the slack of fifty is two things, and the residual’s direction accounts for one of them. Aiming closes about seven of it. The rest — a factor of nine for Householder — is not in the residual at all, since no residual on this problem does worse than the aimed one. It is in the rounding: the bound charges the full size of a backward error $\delta A$ of order $u\lVert A\rVert$ and lets it act in the worst direction, and Householder’s actual backward error is smaller than its bound and does not point that way. A residual cannot be aimed at a quantity the solver has not yet committed.

One direction and then rounding #

Why seven, and why the same seven for three different algorithms, is visible in $E$’s singular values. On the generic matrix Householder’s are 0.82, 0.021, 0.0019, 1.4·10⁻⁴ and so on, falling by a factor of ten or more at each step; the first carries 99.9 per cent of the map’s squared size. The other two routes have the same shape. The spacing is not an accident of rounding. This matrix’s singular values are spaced geometrically, a factor of 3.4 apart, and the $\kappa^2$ term acts through $(A^{\mathsf T}A)^{-1}$, whose eigenvalues are their squares — a factor of 11.7 apart. Each step down $E$’s spectrum is one step up $A$’s, squared. The smallest singular value of $A$ dominates the error by an order of magnitude over the next one, so the map is, for practical purposes, one direction.

A map of rank one turns a residual into error only through its component along one unit vector $w$. A random unit vector in 48 dimensions has a component along any fixed $w$ of root-mean-square size $1/\sqrt{48}$, and the aimed residual is $w$ itself, so the aimed residual draws $\sqrt{48} = 6.93$ times the root-mean-square direction. The three routes measure 6.93, 6.92 and 6.93.

That is a prediction about the number of rows, and it holds.

At 32 rows the gain is 4.00 against $\sqrt{16} = 4.00$; at 128 rows, 10.58 against $\sqrt{112} = 10.58$. The aimed constant itself barely moves — Householder’s is 0.089, 0.113 and 0.137 — while the typical one falls as more rows give a random residual more directions to spend itself in. The earlier essay’s slack was therefore partly a property of the shape of the problem. On a tall fitting problem, a thousand observations of sixteen parameters, a residual made of independent noise would sit thirty times below the residual that the same solver turns into the most error, and it would look correspondingly far below the bound.

That is the same arithmetic a first vector nobody can build against relied on from the other side: a random start vector finds the direction a condition estimator needs with a component of about one over the root of the dimension, which is small but never zero, and so it cannot be fooled. Here the random direction is the residual and the thing it fails to find is the solver’s worst case.

The dial tests the obvious way to break the rank-one picture. Make the smallest singular values of $A$ equal — two, four or eight of them at 1/κ — and there is no longer one direction for $(A^{\mathsf T}A)^{-1}$ to prefer. If the leak were spread evenly over the cluster, the map would have as many equal leading singular values as the cluster has members and the gain would fall to $\sqrt{48/k}$: 4.9, 3.5 and 2.4. It does not. With eight equal singular values Householder’s map reads 1.3, 0.60, 0.17, 0.11, and the gain is 6.2; the appended block’s is 6.5 and column Gram–Schmidt’s, the one that comes closest to the prediction, 4.4. The rounding inside the cluster is not isotropic. Some directions of the cluster leak much more than others, and which ones is a property of the factorisation — the next section’s subject.

Aimed at one route, mild for another #

If the aimed direction were a property of the problem — “the weak directions” — every route would share it, and a residual aimed at one would be the worst for all three. Feed each route’s aimed residual to the other two and that is not what happens. On the generic matrix the residual aimed at Householder draws 0.11 from Householder, 0.031 from the appended block and 0.0010 from column Gram–Schmidt, which is a thirtieth of column Gram–Schmidt’s own worst. The residual aimed at column Gram–Schmidt draws 0.033 from it and 0.0035 from Householder. With the ill-conditioning between blocks the pattern repeats, with the appended block’s aimed residual drawing 0.084 from it and 0.0068 from Householder. In every one of the six cases some other route is under a fifth of its own worst.

All three maps act through the same $(A^{\mathsf T}A)^{-1}$, so all three are dominated by the same smallest singular direction of $A$ on the solution side. What differs is the other side: which residual direction each factorisation lets leak into that singular direction. Householder’s computed range basis is tilted towards the complement in one direction; modified Gram–Schmidt’s, in another; the appended block’s projections, in a third. The tilts are each of size about $\kappa u$ and their directions are set by the order in which the algorithm commits its roundings — the same reason two Gram–Schmidts found two algorithms with identical algebra losing orthogonality by eight orders apart.

One ranking survives every aim. Column modified Gram–Schmidt on $[A\;b]$ has the smallest aimed constant of the three on both placements — 0.033 against Householder’s 0.113 on the generic matrix, 0.011 against 0.014 between blocks — although its Q is the furthest from orthogonal. Orthogonal is a number measured that loss on the Hilbert matrix and found it invisible in the reconstruction; here it is invisible in the residual term as well, because what leaks the residual into the solution is not the angle between the computed columns but the angle between the computed range and the true one, and modified Gram–Schmidt, reading $b$ as its last column, never multiplies $b$ by its non-orthogonal Q at all.

That makes “the worst residual” a phrase with a hidden argument. A test that compares two least-squares solvers on one residual is comparing them on a direction that may be the worst for one and mild for the other, and the comparison can come out either way by choice of residual alone.

Where the appended block falls behind #

The earlier essay concluded that the appended block never falls behind Householder by more than a factor of four, and most of its measurements had it ahead. With the ill-conditioning between blocks of four that conclusion depended on the residual. Householder’s random draws there have a median constant of 0.0016; the appended block’s, 0.0082 — five times larger, on the typical direction. The earlier fixed direction happened to be kind to Householder on this matrix, 5.8·10⁻⁴, a third of its median, while the appended block drew an ordinary 0.015.

The earlier essay’s numbers were medians over three matrices of the same construction, taken separately for each route, and on this placement the three matrices disagree by a factor of three hundred for one residual recipe: Householder’s error at a residual of one per cent is 6.4·10⁻⁶ on the first, 2.5·10⁻⁴ on the second and 2.1·10⁻³ on the third. A median of each route over those three is a median over three different alignments of one fixed residual with each route’s own leak, and the ratio of two such medians compares different matrices. Read one matrix at a time, through the map, the comparison is cleaner: on a typical residual the appended block is five times behind Householder here, and on the residual aimed at it, twelve — 9.3·10⁻⁴ against 7.5·10⁻⁵, on the same right-hand side.

The appended block’s real advantage was never over Householder. It was over its own $Q^{\mathsf T}b$, which loses $\kappa^2 u$ at any residual because its Q is not orthogonal — the rescue what the appended block inherits measured and the right-hand side as one more column first found in columns. That advantage is the one with a threshold.

The one-per-cent rule, aimed #

The earlier essay’s practical result was that appending $b$ is worth its extra block — a factor of ten or more over $Q^{\mathsf T}b$ — up to a relative residual of about one per cent. On this matrix with the fixed direction the advantage is 34 at one per cent and falls to ten at about three per cent, comfortably inside that rule.

Aimed at the appended block, the curve has the same shape and sits lower. Both lines fall as one over the residual once the $\kappa^2$ term takes over, because $Q^{\mathsf T}b$’s error does not move and the appended block’s rises linearly; aiming multiplies the appended block’s constant by 5.6, from 0.015 to 0.084, and moves the line down by the same factor. The advantage at one per cent is 5.3, and it falls to ten at a residual of about half a per cent. The rule shrinks by the factor the aim gains, which is what the earlier essay predicted it would do and is the half of its question that comes out as asked.

Whether a real residual is aimed is a separate matter, and this essay’s measurement says it usually is not. A residual made of measurement noise is a random direction in the complement, and a random direction draws a median of a tenth of the aimed one here; the fixed direction behaved like a random one on the appended block’s map. The aimed residual is the worst case a solver can meet, constructed after the factorisation by someone who has seen its rounding, and a code choosing whether to append $b$ can treat the one-per-cent rule as typical and half a per cent as safe.

What these problems do not show #

Two placements of the ill-conditioning, one matrix of each, one condition number of 10⁸ and one residual size for the map; the map is linear in the residual only to first order, which holds while the $\kappa^2$ term dominates the $\kappa u$ term — from ρ of about 1/κ up — and the aimed direction found at one per cent was used unchanged across the residual sweep. The rank-one structure is a property of matrices whose smallest singular value is separated from the next; the cluster test shows it survives equal singular values, without explaining why. And everything here is sequential: a reduction that changes the order measured a tall-skinny QR done as a tree, whose rounding is committed in a different order, and whose leak would point somewhere else again.

Still open: what tilts the range, and the tree #

The size of the leak. Aiming accounted for a factor of seven of the slack and left a factor of nine for Householder. The remaining factor is the size of the tilt between the computed range and the exact one, along the one direction the map reads. It can be measured directly — the principal angle between the span of Householder’s first sixteen columns of Q and the exact range of $A$, projected on the smallest singular direction — and the prediction with a sign is that this angle times κ, divided by $\kappa^2 u$, reproduces the aimed constant of 0.11 to within the factor the next singular direction contributes, and so reduces the bound’s slack to a statement about how far a backward-stable QR’s computed range actually moves.

The tree. The question both earlier essays left stands, with a sharper form. 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.