Stable once, and three thousand times

The essays on a rank-one correction measure one update. What it costs, when the shortcut is safe, and why removing an observation whose leverage is near one cannot be done at all — every claim there is about a single correction, which is how the literature states them and how a textbook proves them.

A recursive least squares over a sliding window is not a single correction. It adds one row and removes one row at every step, it does that for as long as the data keeps arriving, and it never looks at the rows again. After three thousand steps the triangular factor in memory has been through six thousand rotations and has not seen the matrix since the first one.

Nothing about a per-step bound says what that factor is a factor of.

An exact reference, which is what makes this a distance #

The measurement is only worth having if the thing it is measured against is right, so the stream is built for that: every entry is an integer times a power of two.

That is not a stylistic choice. An integer times a power of two is exact in a binary format, the product of two of them is exact, and the sum of a few hundred of them is exact — so AᵀA computed in double is AᵀA, and ‖RᵀR − AᵀA‖ is a distance from the answer rather than a difference between two float computations. The claim is checked rather than argued: the Gram matrix is formed a second time in BigInt from the integer parts and the two agree bit for bit at every entry.

The conditioning comes from the exponents, which move κ without moving anything else. Column j is scaled by a power of two, spread so that the Gram matrix has a condition number of 7.2·10⁵ — five orders of ill-conditioning with no inexact number anywhere in the matrix.

Every step is safe #

The two halves of a sliding-window step are not the same kind of operation, and the collection already has an essay about why.

Adding a row is a Givens rotation: orthogonal, c² + s² = 1, no entry can grow by more than the vector being folded in, and there is no case in which it fails.

Removing a row cannot be done that way, because an orthogonal transformation cannot make a matrix smaller. It is a hyperbolic rotation, c² − s² = 1, and its amplification is 1/√(1 − h) where h is the leverage of the row being removed — a number with no upper bound.

Over the three thousand steps drawn here, the worst leverage met is 0.6944 and the worst amplification is 2.72. Every downdate succeeds. By any per-step standard this run is uneventful — there is no near-failure anywhere in it, and a code watching for one would report nothing for the whole run.

And the chain is not #

The drift climbs anyway:

3.5·10⁻¹⁵ at 25 steps → 1.1·10⁻¹⁴ at 625 → 2.8·10⁻¹⁴ at 1,225 → 3.9·10⁻¹⁴ at 2,975

The factor in memory ends 350 times the unit roundoff away from the matrix it is supposed to be a factor of, and it is still climbing when the run stops.

Nothing about that is visible from inside. Every step succeeded. Every rotation was correctly computed. The factor is upper triangular, its diagonal is positive, and a solve with it returns a vector. What has happened is that three thousand roundings have been folded into an object that is never compared against anything, and there is no residual anywhere in the algorithm that would notice.

The bound is right, and it is loose by a square root #

The standard statement about accumulation of this kind is a bound linear in the number of steps: the drift after k steps is at most a constant times k·u. At k = 2,975 that is 3.3·10⁻¹³, and the measured drift of 3.9·10⁻¹⁴ is comfortably inside it.

But the bound’s shape is wrong, and the shape is what somebody reasoning about a long run will use. Fitted on log axes over two decades of steps, the drift grows with a slope of 0.554. Not 1.

The reason is the one this collection already measured about a sum: the roundings have signs. A bound has to assume they all point the same way, because it is a bound; what they actually do is a random walk, and a walk of k steps is √k long. So the bound is correct, it is loose by a square root, and a reader who takes it for the behaviour is wrong by a factor of √k — which at three thousand steps is fifty-five, and at a million is a thousand.

The repair, which is the field’s repair #

Rebuild the factor from the rows every m steps.

At every twenty-five steps the drift is held at about sixty times the unit roundoff for a quarter more work, and the climb becomes a sawtooth whose height is set by the period rather than by the length of the run.

Those four numbers are measured at every step. An earlier version of this table sampled the drift every twenty-five steps, which is the last row’s own period, and the next section is about what that did. That is the difference that matters for a stream: a policy of never rebuilding has a drift that grows without limit in the run length, and a policy of rebuilding periodically has one that does not.

The cost is the exchange rate again. A window refactorisation is O(wp²) and a step is O(p²), so rebuilding at every step would cost a factor of w — the whole point of the recursive form. Rebuilding every m steps costs a factor of w/m, and the table is that ratio priced at three values of m.

Swept continuously the trade is linear, and one number in it does not move at all.

The drift is proportional to the period, so the exchange rate is a constant. Over the same 2,976 steps, refreshing every 1000, 200, 100 and 50 gives 3, 15, 30 and 60 refactorisations and worst drifts of 4.21·10⁻¹⁴, 1.13·10⁻¹⁴, 5.68·10⁻¹⁵ and 4.79·10⁻¹⁵. The product of drift and refactorisation count is 1.3·10⁻¹³, 1.7·10⁻¹³, 1.7·10⁻¹³ and 2.9·10⁻¹³ — flat across the middle, so halving the period halves the drift and doubles the work, with no knee to find.

The one number that never moves is the amplification: 2.72 at every refresh period, including none. That is the quantity a reader might expect the repair to be improving, and it is untouched, because it belongs to the problem — it is what the window’s conditioning does to whatever error reaches it — while the refresh policy decides only how much error there is to amplify. So the two halves of the identity this site is built on are cleanly separated here by a knob: the refresh period moves the backward error and nothing else, and the conditioning sits at 2.72 whatever is done about it. A policy discussion about rebuilding is therefore a discussion about one factor of the product, and the other factor is not available for negotiation.

The far end is where that stops paying. Going from every hundred to every fifty doubles the refactorisations and buys 16% of the drift, because at 5·10⁻¹⁵ the drift is approaching the floor set by the arithmetic rather than by the policy — the same shape as the halving that stops, arriving in a different currency. The useful reading is that the linear range is where the choice lives, and it ends around a hundred steps on this run.

A measurement sampled in step with the thing it measures #

The drift is recorded every sample steps and rebuilt every refresh steps, and when those two are the same number every recording lands immediately after a rebuild. The sweep’s shortest period is 25 and the recording interval was 25.

At a refresh period of 25 the drift sampled every 25 steps is 2.15·10⁻¹⁶. Sampled every step it is 6.74·10⁻¹⁵ — thirty-one times larger. At 50 the two agree to within a factor of 1.5; at 100, 1.5; at 500, 1.1; at never, exactly. Everywhere the two intervals are incommensurate some sample lands late in a tooth and the peak is nearly found. At 25 they coincide and the tooth is never sampled at all.

The immediate correction is that the repair buys a factor of ten rather than three hundred: 6.7·10⁻¹⁵ against never’s 6.6·10⁻¹⁴. That is still worth a quarter more work on a stream that runs indefinitely, because the point of the repair was never the constant — it is that a bounded sawtooth replaces a drift that grows without limit in the run length.

And the alias contradicted this essay’s own mechanism #

It is worth being fair to the original sampling, because sampling every twenty-five steps is not a careless choice. The drift measurement costs a Gram matrix — O(wp²), the same as a refactorisation — so recording it at every step would cost more than the algorithm being measured, and every figure in this essay would take twenty-five times as long to draw. Sampling was the right decision and the period was the wrong number, and the wrong number is the sweep’s own shortest period, which is the one place it could do harm.

The general repair is smaller than sampling everything: make the two intervals coprime. Sampling every seven steps costs a third of what sampling every step does and returns 6.74·10⁻¹⁵ at a refresh of 25 — the same as sampling everything, to three digits, because seven and twenty-five share no factor and the samples walk through the tooth. That is the whole fix, it costs nothing, and it would have made this section unnecessary.

The second consequence is the one worth carrying, because it says the defect was findable from inside.

Fitted against the period with honest sampling, the sawtooth’s height is 6.74·10⁻¹⁵, 7.19·10⁻¹⁵, 8.25·10⁻¹⁵, 1.56·10⁻¹⁴, 3.20·10⁻¹⁴ and 4.22·10⁻¹⁴ at periods of 25, 50, 100, 250, 500 and 1,000 — an exponent of 0.497. The square root, to three digits, and exactly the law the section above derives for the un-refreshed run: the roundings have signs, a sum of them is a walk, and a walk of m steps is √m long. A tooth is a walk of m steps, so its height is √m·u, and it is.

The aliased table said something else and said it loudly. A drift of 2.2·10⁻¹⁶ at a period of 25 against 2.9·10⁻¹⁴ at 500 is a slope of 1.6 — faster than linear, which no accumulation of roundings can be, because even the pessimistic bound that assumes every rounding points the same way is linear. The table was asserting something the essay two sections earlier proves impossible, and nothing compared them.

That is the shape worth taking away. The alias did not produce an implausible number: 2.2·10⁻¹⁶ is twice the unit roundoff, which is exactly what a freshly rebuilt factor should carry, and it reads as a success rather than as an error. What it produced was an implausible relationship, visible only by putting the repair table beside the mechanism section and asking whether the two slopes agree. Neither half looks wrong alone.

The check costs one run at a different sampling interval, and it generalises past this figure: when a measurement has a period and the thing measured has a period, vary one of them. This site has now been caught by two versions of the same thing — a bias quantised by a λ grid, and a sawtooth sampled in phase — and in both cases every individual number was correct.

Why this is the opposite case to the rest of the field #

Every other essay in this field has an outer loop that recomputes its residual from the matrix, and that recomputation is what makes a stale inner object survivable. A Newton step solved to two digits is fine because the next step measures the residual again. A factorisation four members old is fine because the chord iteration measures the residual again.

Here nothing measures anything again. The factor is the state, the state is updated by a recurrence, and the data is thrown away as it goes past. The whole algorithm exists because keeping the data would cost O(w) storage and O(wp²) work per step, and the price of not keeping it is that there is nothing left to check against.

So the two situations look identical from a single step — one update, backward stable, nothing to report — and are opposite over a run. The distinction is not how accurate a step is. It is whether anything downstream is going to look at the data again, and it is worth asking of any recurrence before trusting its per-step bound.

What the leverage essay already said, and what it did not #

The essay on the observation that cannot be removed establishes that a downdate’s error is 1/(1 − h) in the leverage of the removed row, and that near h = 1 the operation is not merely inaccurate but impossible: the square root has no real value and the algorithm reports that the matrix it was asked to produce is not positive definite.

That is a per-step statement and it is the one that governs whether a run fails. This run never gets near it. The window is 24 rows on 6 columns, so the average leverage is a quarter and the worst met is 0.69, and the amplification never exceeds 2.72.

What the leverage essay does not say — and could not, because it is about one step — is that a run made entirely of safe steps accumulates anyway. The failure mode and the drift are independent. Narrowing the window towards w = p makes the failure mode arrive; lengthening the run makes the drift arrive; and a run can have as much of one as it likes with none of the other.

What the drift actually costs a user #

A relative drift of 3.9·10⁻¹⁴ in the Gram matrix is not, on its own, an alarming number. It is worth converting it into the quantity somebody cares about, which is the coefficients the window is being solved for.

The factor is used to solve RᵀRβ = Aᵀy. A perturbation of relative size δ in RᵀR moves β by at most κ·δ in the relative sense, and the Gram matrix here has κ = 7.2·10⁵. So a drift of 3.9·10⁻¹⁴ is worth about 2.8·10⁻⁸ in the coefficients — eight digits, from a run in which nothing went wrong.

Two things follow. The first is that the drift is amplified by the conditioning like every other error in this collection, so the essay’s number is a floor on what it costs and the real cost depends on the design matrix. The second is more uncomfortable: the same run on a stream conditioned at 10¹² rather than 10⁵ would lose the answer entirely, and the run would look exactly the same from inside — same amplifications, same successful downdates, same triangular factor with a positive diagonal.

That is the argument for the refresh being a default rather than a tuning option. It costs a quarter more work at the period that holds the drift at the unit roundoff, and what it buys is that the answer’s accuracy stops depending on how long the program has been running.

The habit this generalises #

There is a shape here worth naming, because this collection has now met it three times in objects that share no arithmetic.

An algorithm that keeps a summary and discards the data has no way to check itself. A sliding window keeps a triangular factor and discards the rows. A conjugate gradient iteration keeps a residual vector and never recomputes it. A left-to-right sum keeps a running total and cannot revisit the terms. In all three the summary is exactly the thing that made the algorithm cheap, and in all three the price is that nothing downstream can notice it drifting.

The repair is the same in all three and it is always the same shape: go back to the data, at some period, and pay for it. Residual replacement in a Krylov method. Refactorisation in a window. Compensated or blocked summation in a sum. Each costs a small, bounded fraction of the work, and each converts an error that grows with the run length into one that does not.

The measurement in this essay is what that fraction buys on one of the three.

What a per-step bound is for #

Nothing above is an argument against per-step analysis, and it is worth saying what such a bound does establish, since the essay spends its length on what it does not.

A per-step backward error bound says that one step of the algorithm is the exact step of a slightly different problem. That is the property that makes an algorithm correct: it rules out a step that computes something else entirely, it is what distinguishes a hyperbolic rotation from a subtraction that happens to work, and it is checkable in isolation, which is why it is the form the literature states.

What it cannot do is compose. Two backward-stable steps applied in sequence are backward stable for the composition only if the perturbation the first one licenses is one the second is still stable against, and for a chain of three thousand there is no such statement to make — the object being perturbed is the state, and the state is the thing being carried.

So the bound and the measurement answer different questions, and the essay is not a correction to the bound. It is the observation that the question a long run asks is the second one, and that nothing in the literature’s usual form answers it.

What is asserted #

That every step is safe — worst amplification 2.72, no downdate failing, worst leverage 0.6944 — so the drift cannot be blamed on a near-failure.

That the chain drifts anyway, to 350 times the unit roundoff, having started an order of magnitude closer.

That the growth is a square root — slope 0.554 against the bound’s 1 — with the drift above √k·u and below k·u at the end of the run, so it is bracketed by the two shapes rather than merely fitted.

That the reference is exact, entry by entry against a BigInt computation, with the largest integer sum in the window checked to be inside the exactly representable range.

And the refusal: the same claim, fed a run that rebuilds the factor at every step. There the drift never leaves the unit roundoff, because there is nothing being carried, and the assertion fails as it must.