Every trip comes back shorter

The plateau is set by the length measured how an applied Householder factor — the product of reflectors $I - \beta vv^{\mathsf T}$ kept as its vectors and never multiplied out — loses orthogonality one reflector at a time. The departure, the worst $\|Q^{\mathsf T}(Qx) - x\|$ over unit $x$, grows for the first dozen or two reflectors and then stops, and where it stops is set by the length $m$ of the vectors: $u\sqrt{m}$ times 3.1 to 3.4, at 96, 192 and 384 rows.

The square root was explained by the rounding. Each reflector’s work is an inner product of length $m$ followed by an update, and to nearest, the rounding errors in an inner product have random signs, so $m$ of them add up like a random walk, to about $\sqrt m$ units. The essay then named the case that would break the explanation: “Directed rounding — every operation rounded toward zero, as some low-precision hardware does — gives errors of one sign. The prediction with a sign is that under rounding toward zero the plateau grows like $m$ rather than $\sqrt m$, by a factor near 2 rather than 1.4 for each doubling of the length, while the count to reach it stays near a dozen.”

Both halves fail, and the reason they fail is a distinction the prediction blurred: errors of one sign relative to what.

Rounding modes that JavaScript does not have #

The arithmetic every essay here runs in is double precision rounded to nearest, and there is no switch for any other mode. So this measurement is made in single precision — the format what a float can hold lays out bit by bit — where every other mode can be emulated exactly. The product of two single-precision numbers fits in a double with no rounding at all. Their sum fits too, with the part that does not fit recovered exactly by the two-sum trick. A quotient or a square root can be bracketed by multiplying the candidate back. With the exact result in hand, or bracketed, each operation’s result is rounded to single precision either to the nearest single or to the largest single not above it in magnitude — toward zero. The emulation is checked against two thousand random operations: the toward-zero result never exceeds the exact one in magnitude, and its neighbour one step away from zero always does.

Every reflector is built in the arithmetic being measured, and applied in it. The matrices are Gaussian with 64 columns and 96, 192, 384 or 768 rows, nine of each; the departure after $j$ reflectors is found by power iteration on the round trip. To nearest in single precision, the earlier law reappears: the departure levels off at 2.7 to 3.4 times $u\sqrt m$, with $u = 2^{-24}$ now.

Toward zero, the departure after 64 reflectors is 109 units at 96 rows, 129 at 192, 92 at 384 and 178 at 768. The length grows eightfold and the departure by 1.6. The prediction’s line would end at 870 units; even the square root’s would end at 310. On these matrices the departure toward zero grows more slowly with the length than the departure to nearest does — 27 to 97 units, a factor of 3.6 — and it starts four times higher.

One sign relative to the partial sum #

The prediction’s mechanism is real, and it is easy to find where it holds. Take an inner product of 96 to 768 terms, every term positive, and round each operation toward zero. Every partial sum is positive, every rounding cuts it, and the cuts add up: the result falls short of the exact value by 0.34 units of roundoff per term, relative to the sum of the terms’ magnitudes, at every length. That is linear in $m$ with no ambiguity — 35 units at 96 terms and 262 at 768.

Give the terms random signs and the picture changes completely. Toward zero, a rounding cuts the partial sum towards zero, and the partial sum of a mixed-sign inner product changes sign as it goes — it wanders, positive for a stretch, negative for a stretch. Each cut is towards zero relative to the current partial sum, which is a different direction each time the sum crosses zero. The median relative error is under one unit at every length, which is to say there is no bias, and the typical size grows by a factor of 1.44 to 1.53 for each doubling, the square-root law again. It is larger than to nearest’s — 4 units at 96 against a quarter of one — because each error is a whole unit rather than half of one on average, and because runs of one sign are long enough to add up within each run.

The direction the error leans measured this distinction from the other side: rounding toward infinity accumulated with an exponent of 1.01 in the number of operations, where to nearest gave 0.47, on a sum of positive numbers. “Errors of one sign” means one sign relative to the running value, and whether that is also one sign in absolute terms depends on whether the running value keeps its sign. In a Householder reflector it does not. The inner product $v^{\mathsf T}x$ that sizes each update is a sum of terms of both signs, and so is every entry of the update. Three walks and one bound found three unrelated accumulations — a sum, a chain of rotations and a residual recurrence — all growing with an exponent of a half under rounding to nearest; the mixed-sign inner product rounded toward zero joins them, with a larger constant and the same exponent, and the order they are added in is a reminder that which partial sums appear, and so which signs they carry, is a choice the summation makes.

The one-signed sums are there, and cancel #

There are sums of one sign inside a Householder reflector, and they behave exactly as the prediction expected. Each reflector is $I - \beta vv^{\mathsf T}$ with $\beta = 2/v^{\mathsf T}v$, and $v^{\mathsf T}v$ is a sum of squares — every term positive. Toward zero it comes out short, and $\beta$ comes out long: on the first reflector of nine matrices at each length, $\beta$ is too large by 54 units of roundoff at 96 rows, 106 at 192, 214 at 384 and 423 at 768. That is 0.55 units per row, and it is linear in the length with no ambiguity. A reflector built with that $\beta$ is not a reflection; it overshoots along $v$, and a round trip through it would move a vector along $v$ by four times that error — 1,700 units at 768 rows.

It does not, because the vector it overshoots on is $v$ itself, and to reflect $v$ the factor computes $v^{\mathsf T}v$ again, as the inner product that sizes the update — the same sum of squares, rounded the same way, short by the same amount. The coefficient that decides the reflection is $\beta$ times that inner product, a ratio of two equally biased sums, and it comes out as $2 - 2u$ on every one of the 36 reflectors measured: one unit from exact, at every length. The bias that the prediction located correctly is divided out by the arithmetic that created it. One number that has to be right found $\beta$ the one quantity in a reflector whose error cannot be shrugged off; toward zero it is wrong by hundreds of units and the reflection is right anyway, for the vector where the error would show.

A trip that comes back shorter #

If the inner products do not carry a bias, something else must account for the departure toward zero being four times nearest’s on short vectors, and for the second half of the prediction failing too — that the count to reach the plateau stays near a dozen. It does not stay anywhere.

On 96 rows, to nearest the departure is 26 units after sixteen reflectors and 27 after sixty-four: a plateau, reached in about a dozen as before. Toward zero it is 55 after sixteen and 109 after sixty-four, and every sample in between is higher than the one before. At every length the toward-zero departure after 64 reflectors is at least 1.45 times its value after sixteen; to nearest it is at most 1.35 times, and at 96 and 192 rows it is level or lower. Turn the dial up and the two curves come closer, because the nearest part grows with $\sqrt m$ and the other part does not.

The other part is a shrink. Send a random unit vector through the first $j$ reflectors and back, and measure only its length.

To nearest the length does not change, within two units either way at any $j$. Toward zero, it falls steadily: 20.3 units shorter after sixteen reflectors on 96 rows, 21.3 on 768 rows, and 58 and 82 units after sixty-four. Every application of a reflector updates every entry as $y_i = x_i - d\,v_i$, and the last operation in that is an addition rounded toward zero, which shortens $y_i$ by up to a unit of its own size. Every entry, every application, the same direction relative to the entry’s own sign — and an entry’s own sign, unlike a partial sum’s, does not wander. Across a vector of length $m$ the shortening is a fraction of a unit of the vector’s length, whatever $m$ is, and it repeats twice per reflector for the round trip. Over the first sixteen reflectors it comes to 1.3 units a reflector at every length measured.

The shrink is a floor #

A vector that comes back shorter by $s$ units has moved by at least $s$ units, so the shrink is a lower bound on the departure, direction by direction. The worst direction can only depart further.

That is what the measurements show at every length. Toward zero, the departure after 64 reflectors is 109 against a shrink of 58 at 96 rows, 129 against 71 at 192, 92 against 79 at 384 and 178 against 82 at 768. On short vectors the departure is the shrink plus about what nearest would give; on long ones the shrink is a smaller share, because nearest’s random part grows with the length and the shrink barely does. And the shrink grows with the number of reflectors without limit, so the toward-zero factor has no plateau: each further reflector costs every vector another unit or so of its length. The factor nobody forms found the applied factor’s departure saturating because its worst direction is touched by only a few reflectors; a shrink touches every direction, and saturation needs a direction to run out of.

The obvious test of whether the shrink is the whole of the growth is to take it away. Rescale every reflection’s output to the length of its input — one norm computed per application, one multiplication per entry — and measure again. After 64 reflectors the departure falls from 109 to 74 units at 96 rows, from 129 to 74 at 192, from 92 to 87 at 384 and from 178 to 112 at 768: a third of it at three lengths and a twentieth at the fourth. And it does not level off. On 96 rows the rescaled departure is 20 units after sixteen reflectors and 74 after sixty-four, still rising at every sample. The shrink is the part of the toward-zero error that lies along the vector itself; rounding toward zero also cuts each entry by a different fraction of its own size, depending on where that entry’s value falls between two singles, and that uneven cut points off the vector, in a direction that changes with the entry pattern of every reflector. To nearest the same unevenness has no preferred sign and averages away into the plateau; toward zero it is always a cut, and it accumulates with the count.

So the prediction’s two numbers fail for one reason. It expected one-signed errors to pile up in the inner products, where they would have grown with the length and saturated with the count. They pile up in the update instead, entry by entry, where they grow with the count and hardly at all with the length.

What a shrink costs a factorisation #

A uniform shrink is the most benign error an orthogonal factor can have. $Q^{\mathsf T}Q = (1 - \epsilon)I$ is not orthogonal, but it is a scalar multiple of an orthogonal matrix: angles are preserved, the factor’s columns are still exactly perpendicular in the directions the shrink acts on, and a least-squares solution computed with it is scaled rather than rotated. The measured shrink is not exactly uniform — it is a median over random vectors, and the worst direction departs further — but the part of the toward-zero departure that grows without limit is this one. After 64 reflectors it is 58 to 82 units, about $5 \cdot 10^{-6}$ in single precision; a factor of a thousand reflectors would lose 1,300 units or so, near $10^{-4}$.

A reflection cannot stop being one argued that a computed reflector is exactly a reflection for slightly perturbed data, which is what makes Householder’s method backward stable. Under rounding toward zero each computed reflector is exactly a reflection followed by a slight contraction, and backward stability survives that too, with a constant that grows with the number of reflectors. The standard analysis bounds the departure by a constant times $mnu$ under any rounding mode; the measurement here says that under rounding toward zero the $n$ is the part that shows, and the $m$ is not.

A coin flip that fixes the average found stochastic rounding restoring the right answer to a sum that rounding to nearest froze. Rounding toward zero is its opposite case: the average is wrong in a direction that never changes for any one entry, and the cure for that kind of error is not more precision in the inner products but a correction of each vector’s length — and the rescaled measurement above says that even that removes only part of it, because the cut is uneven across the entries as well as one-signed along the vector.

What nine matrices of four lengths do not show #

The matrices are Gaussian, so their columns point in independent directions and the inner products that size each update have mixed signs. A matrix whose columns share a dominant direction would give inner products with a consistent sign for every reflector after the first, and there the prediction’s mechanism — a bias in the inner products growing like the length — could appear after all. The departure is a power-iteration estimate over nine seeds and is noisy, which is why the doublings jump about, and why the comparison rests on the eightfold change in length rather than on any one doubling. Single precision stands in for double: the structure of the errors is the same, the size of $u$ is not, and nothing here was measured with a hardware rounding mode rather than an emulated one.

Still open: columns that share a sign, and rounding that renormalises #

Columns with a common direction. Make every column a shared unit vector plus a small random part, so that $v^{\mathsf T}x$ has the same sign for most of its terms after the first reflector. The prediction with a sign is that toward zero the inner products’ bias then appears, the departure after 64 reflectors grows by at least a factor of three from 96 rows to 768 — more than to nearest — and that to nearest the plateau keeps its square root, so that the length matters under directed rounding only when the data give the partial sums a sign.

Which half builds the rest. A factor rounded toward zero has its reflectors built toward zero and applied toward zero, and the measurement here does not separate the two. The prediction with a sign is that reflectors built to nearest and applied toward zero give a departure after 64 reflectors within a fifth of the all-toward-zero figures at every length, and that reflectors built toward zero and applied to nearest give one within a fifth of the all-nearest figures — so that the growth without a plateau belongs to the application, entry by entry, and the reflectors’ own one-signed sums contribute nothing that survives.