The sum can be single, the residual cannot

A correction that reads every row took a sliding least-squares window — 24 rows over six nearly parallel columns, a new row in and an old one out at every step — and asked whether its correction could read fewer rows. The window moves its answer by the update its two changed rows imply, then corrects it once: $d = G^{-1}A^{\mathsf T}(y - A\beta)$, a Newton step on the window’s least-squares problem, with $G = A^{\mathsf T}A$ held as a Cholesky factor. Read from any subset of the rows, even 23 of 24, the correction converged to the subset’s own answer, which on noisy collinear data is wrong by a relative error of order one. Every row has to be read.

The essay ended on the other way to make the correction cheaper. “A correction’s cost is reading every row and forming the residual. Computed in single precision and accumulated in double, the residual costs less to form and still sums over every row. The prediction with a sign is that a full correction with a single-precision residual ends within a factor of ten of the double-precision one at $\kappa(A) = 3 \cdot 10^4$ and more than a hundred times above it at $3 \cdot 10^6$, where the residual’s own rounding, amplified by $G^{-1}$, is the larger term.”

The second half holds, by five decades more than it asked. The first half fails by seven. The mechanism it named was wrong in a way that matters, because the rounding of the residual is not amplified by the square of the conditioning, and the half of the correction that is amplified that way turns out to be the half that can sometimes be made cheap.

Three ways to form a correction #

The streams are the earlier essays’: rows of six integer columns built from one large base value plus small perturbations of a chosen spread, so the columns are nearly parallel and the window’s rows have a condition number set by the spread. Across the windows of a stream $\kappa(A)$ ranges over a factor of two, from 1.6 to 3.4 times its median; the earlier essays named the worst stream by the top of that range, $3 \cdot 10^6$, and its median is $2.2 \cdot 10^6$. This essay uses four streams, at medians of about $2 \cdot 10^3$, $2 \cdot 10^4$, $2 \cdot 10^5$ and $2 \cdot 10^6$, with an integer response carrying noise from none to the size of the data. Every number is checked against the window’s exact least-squares answer, computed in rationals from the integer data, at every fiftieth step.

Each step’s correction reads all 24 rows and is formed four ways. In double, as before. With the residual $r_i = y_i - a_i^{\mathsf T}\beta$ formed in single — every operand rounded to binary32 and every product and difference rounded too — and $A^{\mathsf T}r$ then summed in double. With the same single residual but the answer split into its binary32 head and the double remainder, so that only the head’s products are rounded. And the other way round: the residual in double, and the sum $A^{\mathsf T}r$ formed in single.

Seven decades at both conditionings #

At the worst conditioning and noise one the double correction ends at $3.3 \cdot 10^{-9}$, the number the earlier essay measured, and the single residual at $3.7 \cdot 10^{-2}$: eleven million times worse, the hundredfold the prediction asked for with five decades to spare. At $\kappa(A)$ about $2 \cdot 10^4$ the double correction ends at $2.2 \cdot 10^{-12}$ and the single residual at $3.1 \cdot 10^{-4}$, a hundred and forty million times worse, where the prediction said ten. It is not the noise: at that conditioning the single residual ends between $6 \cdot 10^{-5}$ and $3 \cdot 10^{-4}$ at every noise level from one to none, and the double between $10^{-13}$ and $2 \cdot 10^{-12}$. Six and a half to nine decades separate them on every noisy stream, at every conditioning.

The single residual’s error falls by a factor of 9 to 12 for each decade of $\kappa(A)$: $3.7 \cdot 10^{-2}$, $3.6 \cdot 10^{-3}$, $3.1 \cdot 10^{-4}$ and $3.4 \cdot 10^{-5}$ from the worst stream to the best. That is one decade of error per decade of conditioning, not two, and it sits on a line parallel to single’s unit roundoff times $\kappa(A)$: between 0.04 and 0.5 of $2^{-24}\kappa(A)$ on every noisy stream. The double correction sits on the matching line for double at three of the four conditionings and a decade above it at the worst.

Why the residual’s rounding is amplified once #

The prediction reasoned that the residual’s rounding is “amplified by $G^{-1}$”, and $G^{-1}$ carries the square of the condition number. But the rounding does not reach $G^{-1}$ on its own. It enters as a perturbation $\delta r$ of the residual, and the correction applies $A^{\mathsf T}$ to it before it applies $G^{-1}$: the step’s error is $G^{-1}A^{\mathsf T}\delta r = A^{+}\delta r$, the pseudo-inverse applied to the perturbation. The pseudo-inverse’s norm is $1/\sigma_{\min}$, which is $\kappa(A)/\|A\|$, one power of the conditioning and not two. Half of $G^{-1}$'s amplification is undone by the $A^{\mathsf T}$ in front of it, and that is the same reason the seminormal equations, which solve through $R^{\mathsf T}R = A^{\mathsf T}A$, are not as bad as the road that squares the problem makes the normal equations look: what squares the conditioning is forming the Gram matrix in floating point, and the correction here never forms anything new from it.

What sets the size of $\delta r$ is not the residual. Each $a_{ij}\beta_j$ is a product of the size of the rows times the answer, and the residual is what is left when those products nearly cancel against $y_i$. Rounded to single, each product carries an error of about $2^{-24}|a_{ij}||\beta_j|$, whatever the residual is. On the worst stream at noise one the single residual formed at the exact answer is wrong by about one per cent of the largest residual; at noise $10^{-6}$, where the residual is far smaller than the products that make it, by thirteen per cent. So the single residual’s error is $2^{-24}$ times the size of the rows and the answer, carried by $A^{+}$, and that is $2^{-24}\kappa(A)$ times the answer’s size, up to a constant: the floor the figure shows, independent of the noise.

Measured at the exact answer, where the double residual is the true one, the single residual’s error has two regimes. With noise the size of the data it depends on the conditioning: about one per cent of the residual at $\kappa(A)$ near $2 \cdot 10^6$, a tenth of that at $2 \cdot 10^5$, and down to $1.4 \cdot 10^{-5}$ at $2 \cdot 10^3$, because the window’s answer is larger the more nearly parallel its columns are, and so are the products that have to cancel. As the noise falls the residual shrinks and the products do not, and at $10^{-6}$ every conditioning meets at the same place: the residual formed in single is wrong by 12 to 16 per cent of itself. A residual that small is a few units left over from products of millions, and binary32 spaces numbers of a million to a few million between a sixteenth and half a unit apart.

What a float can hold is the whole of the reason the size of the products matters and the size of the residual does not. Binary32 numbers near a product of the rows’ size are spaced $2^{-24}$ of that size apart, and the residual is read off as a difference between numbers of that size. A small residual is a short distance measured on a coarse ruler: it does not make the ruler finer. That is also why the single residual’s error is so much larger, relative to the residual, at low noise than at high — thirteen per cent against one — and why it makes no difference to the floor. The correction does not care how wrong the residual is as a fraction of itself; it cares how wrong it is in absolute terms, and that is fixed by the products. A small residual is not a small error made the converse point about checking an answer by its residual: the residual’s smallness is not information about the answer’s accuracy. Here it is not even information about the residual’s own.

Buying the accuracy back found the same verdict for a linear system: a factorisation in single corrected by residuals in double recovers a double-precision answer, and with residuals in single the identical algorithm recovers nothing. The window adds two things to that. The floor has a slope, one decade per decade of conditioning, which says how much a single residual costs at any $\kappa$ rather than whether it costs. And the earlier essay’s mechanism of amplification by $G^{-1}$ — which was right for what it explained, the partial correction’s fixed point at the subset’s answer — does not carry over to rounding, because rounding arrives through $A^{\mathsf T}$ and a missing row does not.

Splitting the answer does not help #

If the trouble is that the products carry the answer’s size, the obvious repair is to keep the answer’s low part out of the rounding: round only its binary32 head, and add the tail’s contribution in double. It changes nothing. On every noisy stream the split’s error is within a factor of two of the plain single residual’s: $3.2 \cdot 10^{-4}$ against $3.1 \cdot 10^{-4}$ at noise one, $7.2 \cdot 10^{-5}$ against $6.1 \cdot 10^{-5}$ at $10^{-4}$. The head’s products are still rounded to single, and the head carries almost all of the answer’s size, so almost all of the rounding is still there.

The one place the split reads zero is the stream with no noise at the two better conditionings, and it is a coincidence of the test rather than a property of the split. With no noise the exact answer is a vector of ones and minus ones, which binary32 represents exactly; once the window has converged its head is the answer, its tail is zero, and every product it rounds is an integer times one. On the worst stream the same coincidence does not save it, and on any real stream the answer is not a vector of small integers.

The other half can be single, where the residual is small #

The fourth way of forming the correction is the reverse: the residual in double, so its cancellation happens at full precision, and the sum $A^{\mathsf T}r$ in single. On the stream with no noise it is free. At every conditioning it ends within a factor of two of the double correction — $1.2 \cdot 10^{-11}$ against $1.4 \cdot 10^{-11}$ at the worst, $1.7 \cdot 10^{-13}$ against $1.6 \cdot 10^{-13}$ at $2 \cdot 10^4$. On a noisy stream it is not, and the way it fails is the opposite of the residual’s.

At noise $10^{-6}$ the single sum’s error is $4.9 \cdot 10^{-3}$, $3.9 \cdot 10^{-5}$, $5.4 \cdot 10^{-7}$ and $5.9 \cdot 10^{-9}$ from the worst conditioning to the best: a factor of 73 to 124 for each decade of $\kappa(A)$, two decades a decade. And at every conditioning it grows with the noise — at $2 \cdot 10^3$, from $5.9 \cdot 10^{-9}$ at noise $10^{-6}$ to $2.9 \cdot 10^{-7}$ at $10^{-4}$ — until at noise of the size of the data it meets the single residual’s floor.

This is the amplification the prediction had in mind, attached to the other half. The rounding of a sum of products $a_{ij}r_i$ is about $2^{-24}$ times the sum of their sizes, which is proportional to the residual, and it lands in $c = A^{\mathsf T}r$ directly. Nothing applies $A^{\mathsf T}$ to it first, so $G^{-1}$ amplifies it by the full square of the conditioning: an error of $2^{-24}\kappa(A)^2$ times the residual’s size over the rows’. When the data fit exactly the residual goes to zero and takes the rounding with it. When they do not, the rounding is the residual’s size times $\kappa^2$.

So the two halves of a correction cost in different currencies. The residual’s rounding is set by the size of what cancels — the rows and the answer — and is amplified once; it costs $2^{-24}\kappa$ on every stream, whatever the noise. The sum’s rounding is set by the size of what survives the cancellation — the residual — and is amplified twice; it costs $2^{-24}\kappa^2$ times the residual, which is nothing on consistent data and everything on noisy data at high conditioning.

Along the stream #

Neither cheap correction drifts. At noise $10^{-4}$ and $\kappa(A)$ about $2 \cdot 10^4$ the single residual’s error has a median of $6.1 \cdot 10^{-5}$ and a worst of $3.1 \cdot 10^{-4}$ over the sampled steps, a floor that is there from the start, and the double correction’s sits eight and a half decades below it. That is what stable once and three thousand times and the repair the drift did not need found of the window itself: every step’s error is set afresh by that step’s arithmetic, because each correction is a contraction toward the window’s answer, and nothing accumulates. Turn the dial down through the noise levels and the single sum’s line falls toward the double one, reaching it only at none, while the single residual’s line stays where it is.

What the window was asking for #

The proposal wanted a correction that costs less and still reads every row. The measurement says which half of it can be cheapened, and when. The residual cannot, at any conditioning these streams reach: a single-precision residual puts a floor of a few hundredths of $2^{-24}\kappa(A)$ under the answer, which at $\kappa = 2 \cdot 10^3$ is already $3 \cdot 10^{-5}$ — more than any double correction’s error by eight decades. The sum can, but only where the residual is small enough that $2^{-24}\kappa^2$ times it is below what the answer needs. On these streams that is the consistent one, and close to it.

The step the two rows owe found that the update the window starts from is already the right move on data that fit and three orders short on data that do not, and the answer the last window left drew the same line for carrying the previous answer forward. The precision of the sum falls on the same side of the same line: everything a sliding window can economise on is free when the data agree exactly and has a price proportional to the disagreement otherwise. The residual’s precision is the exception. It is paid on every stream, because what it rounds is not the disagreement but the size of the data.

It is worth being exact about what the cheap half would buy, since the proposal was about cost. A correction reads every row twice: once to form the residual, six products and a subtraction a row, and once to form the sum, six more products a row; then it applies the factor, two triangular solves on 21 numbers. The rows are most of the traffic, and both halves read them. Making the sum single halves the precision of one of the two passes over the rows and leaves the other, so on a machine where the cost is memory it is a quarter of the correction’s traffic at most — and it is free only where the data nearly fit. The pass that would have saved the most is the residual, which reads the rows and the response together, and it is the one that cannot be made cheaper at all.

What four streams do not show #

One family of streams, six columns, a window of 24, integer data whose exact answer can be computed. On data whose answer is small beside the rows — a well-scaled problem rather than a collinear one — the products that cancel in the residual are smaller and the single residual’s floor lower in absolute terms, though the measurement here says its slope in $\kappa$ would be the same. No fused multiply-add: a residual formed with one rounding per product-sum would round less, and how much less on a cancellation this severe is not measured. And the cost the proposal was about — memory traffic, which is what single precision saves on a real machine — is not measured at all; these are the accuracies, and whether either cheap half is worth its saving depends on a machine this essay does not have.

Still open: a single factor, and a wider window #

A factor in single, a residual in double. The half of the correction the measurement has not touched is the solve: the Cholesky factor of $G$, used in two triangular solves per step. Buying the accuracy back found a single factorisation corrected by double residuals recovers double accuracy on a linear system, as long as the single factor is good enough to contract. The prediction with a sign is that the window with its factor updated and held in single and its residual and sum in double ends within a factor of ten of the all-double window at $\kappa(A)$ up to $2 \cdot 10^5$, and fails to contract at all at $2 \cdot 10^6$, where $2^{-24}\kappa(A)^2$ passes one.

A wider window. The earlier essay’s prediction for 72 of 96 rows still stands untested. This essay adds one with a sign: with 96 rows the single residual’s floor stays within a factor of two of the 24-row window’s at every conditioning, because the floor is set by $\kappa(A)$ and the size of the answer, neither of which the width changes, while the double correction’s error falls.