The repair the drift did not need

Stable once, and three thousand times slid a window of rows along a stream, adding one row and removing one at every step, and carried a Cholesky factor of the window’s Gram matrix through three thousand updates without ever looking at the data again. No single step amplified by more than 2.72 and no downdate failed, and at the end the factor in memory was 3.9·10⁻¹⁴ from the matrix it was supposed to factor — six hundred times the growth any single step allowed.

That essay then converted the drift into the quantity a user cares about, and the conversion is where this essay starts. The factor is used to solve RᵀRβ = Aᵀy for the window’s coefficients β. A perturbation of relative size δ in RᵀR moves β by at most κ·δ, and the Gram matrix there had κ = 7.2·10⁵. So a drift of 3.9·10⁻¹⁴ was worth about 2.8·10⁻⁸ in the coefficients — eight digits, from a run in which nothing went wrong. Three conclusions followed: that the drift is amplified by the conditioning like every other error in the collection; that a stream conditioned at 10¹² would lose the answer entirely while looking exactly the same from inside; and that a periodic refresh of the factor should therefore be a default rather than a tuning option.

All three were argued and none was measured. They can be measured, because the stream is integers and the coefficients can be computed exactly, and when they are, each comes out differently from how it was stated.

Exact coefficients, which the conversion never had #

That essay measured the drift exactly, entry by entry, against a Gram matrix accumulated in big integers. It did not measure the coefficients, and the conversion from one to the other is a bound — the largest amount a perturbation of that size could move the answer — rather than a measurement of how far this perturbation did.

The coefficients can be had exactly by the same route. The window’s rows are integers and its responses are integers, so AᵀA and Aᵀy are integers, formed without rounding in big-integer arithmetic, and the system AᵀAβ = Aᵀy is solved in rationals. Where the stream’s columns carry a power-of-two scaling, the scaling is undone afterwards by dividing by powers of two, which is exact. The only rounding in the reference is the final conversion of each rational coefficient to a double. An answer that is known is the collection’s standing reason for building problems this way, and here it is what separates a bound from a number.

At each sampled step the window’s coefficients are then solved three ways from the same exact Aᵀy, and each is compared coefficient by coefficient with the rational solution, keeping the largest relative error:

• with the carried factor, the one that has drifted;
• with a factor recomputed from the rows currently in the window — the earlier essay’s repair;
• with the carried factor followed by one step of correction against the window’s own rows: form r = y − Aβ from the rows, solve RᵀRδ = Aᵀr with the same carried factor, and take β + δ.

Householder least squares on the same rows is drawn beside them as the accuracy a backward-stable method attains on each window.

Two ways to make a window ill conditioned #

A window’s Gram matrix can be badly conditioned for two reasons, and they are not the same reason in any sense that matters to the coefficients.

The earlier essay’s stream gets its conditioning from column scaling: its integer columns are multiplied by powers of two, so the Gram matrix’s condition number can be pushed as high as desired while the columns themselves stay as far from parallel as they were. The second stream gets it from collinearity: each column is an integer base vector plus integer perturbations of size spread, so the columns are nearly parallel and no diagonal scaling can separate them. Every entry is still an integer and every Gram entry an integer below 2⁵³, so the reference is still exact.

Both streams can reach any condition number. Only one of them has a condition number the solve can feel.

On the original stream the conditioning never arrives #

Take the earlier essay’s stream, keep its rows, and push the column scaling from nothing to twelve decades. The Gram matrix’s condition number goes from 4.2 to 1.2·10¹² — past the 10¹² at which the earlier essay predicted the answer would be lost.

Medians over three seeds, at the end of the run:

The coefficient error is the same number at every setting, to three figures: 6.91·10⁻¹⁴ at a condition number of 4 and at a condition number of 10¹². The drift reads the same from four decades on. The conversion, drift times κ, would have predicted 3.9·10⁻¹³ at the first row and 2.9·10⁻² at the last; at twelve decades it is pessimistic by eleven orders of magnitude.

The reason is a property of the Cholesky factorisation that the conversion does not see. Multiplying column j of A by 2ᵉ multiplies row and column j of AᵀA by 2ᵉ, and multiplies row and column j of R by 2ᵉ as well — exactly, with no rounding introduced, because a power of two only moves the binary point. The triangular solve then works on the same significands at every scaling, the same roundings happen in the same places, and a relative error in each coefficient is unchanged. A condition number produced by diagonal scaling is a condition number the factor carries exactly and the solve never meets. A condition number scaling cannot move and two condition numbers of one matrix are about precisely this gap between the number κ(AᵀA) reports and the componentwise sensitivity the arithmetic actually has.

On this stream the earlier essay’s repair does work, and it is worth saying so plainly: in the figure above, a recomputed factor gives 1.8·10⁻¹⁶ against the carried factor’s 4.5·10⁻¹⁴, and so does a single correction. The drift is the error here. It is simply an error of 10⁻¹⁴, not 10⁻⁸, and a stream conditioned by scaling at 10¹² loses nothing it did not lose at 10⁰.

A stream whose conditioning is real #

The collinear stream is the one a design matrix usually is: columns that measure nearly the same thing, which no rescaling separates. Its condition number reaches the solve, and the question becomes whose fault the error is.

At κ(A) = 2,025 the carried factor’s coefficients are wrong by 3.5·10⁻¹⁰ while its drift is 9.5·10⁻¹⁵. The recomputed factor, which has no drift at all, is wrong by 1.5·10⁻¹⁰ — better by a factor of two, and nowhere near Householder’s 2.4·10⁻¹³. Whatever is costing the coefficients three orders of magnitude, a fresh factor has it too.

Across six conditionings, medians over three seeds and 1,500 steps:

From κ(A) = 202 upward the recomputed factor is within a factor of five of the carried one at every setting, and on the medians it is slightly worse at five of the six. Both columns rise like κ(A)²u, the last column, and sit within a factor of five of it. Only at the gentlest setting, κ(A) = 20.5, does the refresh improve on the carried factor by a useful factor — 5.7·10⁻¹⁴ to 8.5·10⁻¹⁵ — because only there is the problem well enough conditioned for the drift’s share to be visible above the solve’s.

The mechanism is the one the road that squares the problem measured, arriving in a different setting. A Cholesky factor of the Gram matrix, however accurately formed, is used through the seminormal equations RᵀRβ = Aᵀy, and two triangular solves with R and Rᵀ compose to the operator AᵀA whatever produced R. Its condition number is κ(A)², and a backward-stable solve through it loses κ(A)²u. The carried factor loses that; a fresh factor loses that; the drift, at 10⁻¹⁴, is a rounding error on top of it. The earlier essay’s conversion attributed to the drift an error that belongs to the solve.

That also explains the coincidence in the opening figure. The collinear stream at κ(A) = 2·10⁴ has a carried-factor error of 2.8·10⁻⁸ — the very number the earlier essay predicted for its own stream from the drift. The two share digits and nothing else: one was κ(AᵀA) times a drift of 3.9·10⁻¹⁴ on a stream whose conditioning is a scaling, and the other is κ(A)²u on a stream whose conditioning is real.

One correction from the rows #

The fourth column is the repair, and it is not a refresh. The same essay that measured the seminormal equations’ κ² measured their cure: one step of iterative refinement in which the residual is formed from A itself rather than from the factor. Form r = y − Aβ from the rows in the window, solve RᵀRδ = Aᵀr with the same carried factor, and add δ to β. These are the corrected seminormal equations.

One correction brings the carried factor to within ten times κ(A)u at every setting through κ(A) = 2·10⁵ — 1.3·10⁻¹⁵, 3.5·10⁻¹⁴, 3.9·10⁻¹³, 8.0·10⁻¹² and 3.4·10⁻¹¹ — and in the single-seed figures it lands on Householder’s own accuracy: 7.1·10⁻¹³ against 1.1·10⁻¹² at κ(A) = 2·10⁴, 1.6·10⁻¹¹ against 1.7·10⁻¹¹ at 2·10⁵. The carried factor, drift and all, delivers a backward-stable answer after one pass over the window’s rows.

The reason the correction works where the refresh does not is that its residual reads the data. The error in β from the seminormal solve is of relative size κ(A)²u; the correction computes the residual that error leaves in the rows, where it is of ordinary size, and solves for it with the same factor, which reduces that error by another factor of about κ(A)²u. From a start at κ(A)²u, one step lands at about κ(A)⁴u², and that is below κ(A)u exactly when κ(A)³u is below one.

At κ(A) = 2·10⁶ the product κ(A)³u is about 900, and one correction is not enough: the median lands at 2.0·10⁻⁸, ninety times κ(A)u. A second correction takes it to the level of Householder least squares.

When a second correction helps, and when it cannot #

The κ(A)³u boundary makes a prediction about the second correction as well as the first, and the four figures above test it. Where one correction has already reached Householder’s accuracy there is nothing left for a second to remove, and it should move the answer only within rounding. Where one correction has not, the second should close the gap.

At κ(A) = 2,025, where κ(A)³u is about 10⁻⁶, one correction gives 1.5·10⁻¹³ and a second 2.2·10⁻¹³, against Householder’s 2.4·10⁻¹³: the second step moves the answer by less than the difference between two correct answers. At κ(A) = 2·10⁴, where κ(A)³u is about 10⁻³, one correction gives 7.1·10⁻¹³ and a second 3.7·10⁻¹², against Householder’s 1.1·10⁻¹²: the second is three times worse than the first and still inside a factor of four of Householder, which is what rounding noise on an answer already at the backward-stable level looks like. At κ(A) = 2·10⁵, where κ(A)³u is about 0.9, one correction gives 1.6·10⁻¹¹ and a second 5·10⁻¹² — a genuine improvement by a factor of three, at the boundary where the rule says a single step starts to fall short. And at κ(A) = 2·10⁶, where κ(A)³u is about 900, the median goes from 2.0·10⁻⁸ after one to 3.2·10⁻¹⁰ after two, below Householder’s 1.3·10⁻⁹.

So the rule is not only a sufficient condition that happened to hold. It separates the settings where the second step is noise from the settings where it is necessary, at the boundary the arithmetic predicts. The same structure is in the gap refinement can close, where refinement at the working precision repaired a solve through a computed inverse in one step at every conditioning measured: a correction whose residual is computed from the data contracts by the error of the solve it reuses, and the only question is whether one contraction is enough.

It also bounds what the correction cannot do. The correction reads the rows the window holds and reuses the factor the window carries, and on the downdate side of the window that factor is maintained by rotations that are not orthogonal. The observation that cannot be removed found that removing a row costs a square root of 1 − h, where h is the row’s leverage, and that the arithmetic of the removal degrades as h approaches one. A correction step repairs the answer against the error the solve introduces. It does not repair a factor that a near-breakdown downdate has damaged beyond the conditioning of the problem, and no stream here came near that case: the worst leverage on the earlier essay’s run was 0.69.

What each repair costs #

A refresh recomputes the window’s factor from its rows: forming the Gram matrix of w rows on p columns costs about wp²/2 and its Cholesky factor p³/6, so the refresh is dominated by wp² — 864 multiply-adds for a window of 24 rows on 6 columns, before the factor is used.

A correction passes over the rows once to form r = y − Aβ and Aᵀr, about 2wp, and does two triangular solves with the factor already in memory, about p². For the same window that is about 330 multiply-adds, and it costs a factor of p less than a refresh as the window grows: at 200 rows on 20 columns, roughly 8,400 against 40,000.

The earlier essay priced its refresh at “a quarter more work at the period that holds the drift at the unit roundoff”. That price bought an answer with the seminormal equations’ κ(A)²u error in it. A correction at every solve costs less than that refresh and buys κ(A)u.

It is worth being exact about what the correction needs, because it is the one thing a sliding window was designed not to keep: the rows. The problem that arrives again describes the sequence setting in which a window never looks at its data twice, and a correction looks at w rows once per solve. A window that is sliding already holds its w rows to remove them later, so in the usual implementation they are in memory; a stream processed with only its factor and running cross-products is a different design, and for it the correction is not available at any price.

What the earlier essay’s three conclusions become #

The drift is amplified by the conditioning. It is amplified by the componentwise conditioning, not by κ(AᵀA). On a stream conditioned by scaling the conversion is eleven orders pessimistic at κ(AᵀA) = 10¹², and the coefficients lose nothing. On a stream conditioned by collinearity the drift is amplified too, but it is dwarfed by the solve’s own κ(A)²u, which is there whether the factor drifted or not.

A stream conditioned at 10¹² would lose the answer. A stream conditioned at κ(AᵀA) = 10¹² by scaling keeps 6.9·10⁻¹⁴. A stream conditioned at κ(A) = 2·10⁶ — κ(AᵀA) about 4·10¹² — by collinearity loses four digits to 1.4·10⁻⁴, and loses them identically with a fresh factor. Both halves of the prediction’s number are right about something; neither is right about the drift.

The refresh should be a default. The refresh is the wrong repair on the stream where repair matters. What a rebuild is worth found that the right rebuild period for a preconditioner depends on what the rebuild costs and not on the drift, and a factorisation kept past its date found a kept factor failing by ceasing to converge rather than by losing accuracy; here the right repair is not a rebuild at all, because what needed repairing was never the factor.

None of this changes the earlier essay’s central measurement. The drift is real, its growth along the chain is a square root of the step count, and a factor carried indefinitely does move away from its data. What changes is the price put on that movement, which the measurement shows was borrowed from a bound and then attributed to the wrong error.

Where this goes from here #

Carrying a triangular factor that never forms the Gram matrix. The carried factor here is updated by rotations that act on R directly, which is the stable choice for the update; the seminormal solve is what squares κ. A window that carries a QR factorisation with Q represented implicitly — or that updates R and applies the correction at every solve as a matter of course — has a different accuracy profile on the collinear stream, and the κ(A)³u boundary for one correction is the natural thing to measure on it.

How often to correct. Every solve above is corrected. A code solving at every step could correct at every kth, carrying the corrected coefficients forward as the starting point, and whether the correction’s benefit decays along the chain as the drift grows — or holds, because it reads the rows each time — is unmeasured.

And a stream whose data drifts. Both streams here are stationary: the drift is in the factor, not in the problem. A stream whose columns change their conditioning along the run, passing through a stretch of near-collinearity and out again, would put the κ(A)³u boundary on a moving target, and the question of when a correction needs a second pass becomes a question the solver has to answer as it goes.