Numerical Linear Algebra
where the algebra stops being the arithmetic
When the problem arrives again

Stable once, and three thousand times

A sliding window adds a row and removes one at every step and never looks at the data again. No single step of it amplifies by more than 2.72, no downdate fails, and after three thousand steps the triangular factor in memory is 3.9·10⁻¹⁴ from the matrix it is supposed to be a factor of — six hundred times growth from a per-step bound that says nothing about chains.
30 min read 24 figures Where the algebra stops being the arithmeticMeasured, not assertedExact ground truthTwo routes to a number

Worth reading first: A correction cheaper than the problem · The problem that arrives again · The exact answer to a nearby problem.

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.

How far a carried triangular factor drifts from the data, over 2976 steps of a sliding windowA window of 24 rows on 6 columns, moved one row at a time: each step folds a row in with Givens rotations and removes one with hyperbolic rotations, and the factor is never rebuilt from the rows. The rows are integers times powers of two, so AᵀA is exact in a double and the vertical axis is a distance from the answer. No single step amplifies by more than 2.72, no downdate fails, and after 2976 steps the factor is 3.87·10⁻¹⁴ from the matrix it is supposed to factor — a fitted slope of 0.554 in the step count, against a bound whose slope is 1.10²10³10⁴10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²steps taken‖RᵀR − AᵀA‖ ⁄ ‖AᵀA‖the bound, linear in the steps√k · uevery step safe, the chain notdrift after the run3.9·10⁻¹⁴the bound there3.3·10⁻¹³√k · u there6.1·10⁻¹⁵worst single amplification2.7worst leverage met0.69refactorisations1backward stable onceand three thousand times is a different claim
Fig. 1 Three thousand steps of a width-24 window on six columns, with the factor never rebuilt. The upper dashed line is the standard bound and the lower one is its square root.

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. 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.

The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 2 The collection’s other exact reference, from the field that established the habit: a rational answer against which a float one can be a distance rather than a difference.
Two condition numbers of one 12×12 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 15 to 3.3·10⁶ while the componentwise condition number stays at 13.2 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.8·10¹⁶ and 1.1·10¹⁶: a matrix whose sensitivity no scaling repairs.0123456110³10⁶10⁹10¹²10¹⁵10¹⁸spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread15κ_∞ at 6 decades3.3·10⁶cond, either end13Hilbert, equilibrated1.3·10¹⁶the solution is the same at every spreadand one of these curves knows it
Fig. 3 And the mechanism the conditioning here comes from: column scales spread apart, which moves the condition number without moving a single digit of any entry.

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.

Downdating a Cholesky factor against refactorising it, as the removed row's leverage approaches oneA rank-one term is removed from a 6×6 Gram matrix by hyperbolic rotations and, separately, by factorising the downdated matrix from scratch. The downdate's residual rises from 2.1·10⁻¹⁶ at h = 0.3 to 3.5·10⁻¹⁰ at h = 1 − 10⁻⁷, a slope of 1.04 against 1/(1 − h). The refactorisation is flat at 8.1·10⁻¹⁷, and the matrix both are producing has a condition number of 4.3 at every point on the axis — so the difficulty belongs to the route and not to the answer.110¹10²10³10⁴10⁵10⁶10⁷10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸1/(1 − h), the leverage of the removed row‖R̄ᵀR̄ − (G − aaᵀ)‖ / ‖G − aaᵀ‖downdatedrefactorisedat h = 1 − 10⁻⁷κ of the downdated matrix4.3κ of the matrix downdated9.3·10⁶rotation's amplification344downdate residual3.5·10⁻¹⁰a hyperbolic rotation is not orthogonaland that is exactly what it is for
Fig. 4 The amplification of a single downdate against the leverage of the row it removes, from the essay that measured it. Everything in this run sits comfortably on the left of that curve.
Leverage and the deleted residual for 40 observations of a 6-column fitThe upper panel is the diagonal of the hat matrix, one bar per observation, with the average p/m = 0.150 drawn through it; the leverages sum to 6.000000000, which is exactly the number of columns. The lower panel is the leave-one-out residual, computed in closed form as eᵢ/(1 − hᵢ) and, separately, by refitting the model 40 times without each observation; the two agree to 1.3·10⁻¹². The first observation carries a leverage of 0.5000 by construction, and 1 − h is the number a hyperbolic downdate takes the square root of.the diagonal of the hat matrix, hᵢ = aᵢᵀ(AᵀA)⁻¹aᵢaverage p/m = 0.15010the leave-one-out residual: eᵢ/(1 − hᵢ), and forty refitsbars: closed form · dots: refitted without that pointone number, two fieldsΣ hᵢ, exactly p6largest leverage0.5closed form against refits1.3·10⁻¹²1 − h of the first row0.5y appears in the residualand nowhere in the leverage
Fig. 5 And where the leverage comes from. A window of 24 rows on 6 columns has an average leverage of a quarter, which is why nothing here comes close to failing.

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.

How far a carried triangular factor drifts from the data, over 2976 steps of a sliding windowA window of 24 rows on 6 columns, moved one row at a time: each step folds a row in with Givens rotations and removes one with hyperbolic rotations, and the factor is never rebuilt from the rows except every 500 steps. The rows are integers times powers of two, so AᵀA is exact in a double and the vertical axis is a distance from the answer. No single step amplifies by more than 2.72, no downdate fails, and after 2976 steps the factor is 1.27·10⁻¹⁴ from the matrix it is supposed to factor, against 3.87·10⁻¹⁴ with no refresh at all.10²10³10⁴10⁻¹⁹10⁻¹⁸10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²steps taken‖RᵀR − AᵀA‖ ⁄ ‖AᵀA‖the bound, linear in the steps√k · uevery step safe, the chain notdrift after the run1.3·10⁻¹⁴the bound there3.3·10⁻¹³√k · u there6.1·10⁻¹⁵worst single amplification2.7worst leverage met0.69refactorisations6backward stable onceand three thousand times is a different claim
Fig. 6 The same run with the factor rebuilt from the rows every five hundred steps, which costs one per cent more work and halves the worst drift.

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.

Relative error of three summation algorithms in binary32A log–log plot of relative error against the number of terms for naive, pairwise and compensated summation, each measured against the exactly rounded sum.10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸
Fig. 7 The same shape in the place this collection first measured it: the error of a left-to-right sum against the number of terms, where the bound is n·u and the behaviour is √n·u.
Three accumulations, each against its own quantity, each divided by its own first pointA left-to-right sum against the number of terms, at a fitted slope of 0.486; a chain of rotations against the number of steps, at 0.554; and a conjugate gradient residual recurrence against the largest iterate, at 0.507. The three share no arithmetic and no vocabulary. Each has a standard bound that is linear in whatever it accumulates against, drawn here as the upper line, and each comes out at half of it. Nothing is rescaled except the division by each series' own first point, which is what makes three quantities of different sizes comparable in slope and in nothing else.110¹10²10³110¹10²10³the accumulating quantity, relative to its first valuethe error, relative to its first valuethe bounds: slope 1what all three do: slope ½three mechanisms, one exponenta left-to-right sum0.49a chain of rotations0.55a residual recurrence0.51every bound's slope1spread of the three0.067a bound is a sum of the roundingsand the roundings have signs
Fig. 8 And the three measurements side by side, which is what the essay on three walks is about. Nothing is rescaled except each series’ division by its own first point.

The repair, which is the field’s repair

Rebuild the factor from the rows every m steps.

refresh worst drift refactorisations extra work
never 6.6·10⁻¹⁴ 1 —
every 500 2.9·10⁻¹⁴ 6 +1%
every 100 5.7·10⁻¹⁵ 30 +6%
every 25 2.2·10⁻¹⁶ 120 +24%

At every twenty-five steps the drift is held at twice 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. 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.

How far a carried triangular factor drifts from the data, over 2976 steps of a sliding windowA window of 24 rows on 6 columns, moved one row at a time: each step folds a row in with Givens rotations and removes one with hyperbolic rotations, and the factor is never rebuilt from the rows except every 100 steps. The rows are integers times powers of two, so AᵀA is exact in a double and the vertical axis is a distance from the answer. No single step amplifies by more than 2.72, no downdate fails, and after 2976 steps the factor is 2.48·10⁻¹⁵ from the matrix it is supposed to factor, against 3.87·10⁻¹⁴ with no refresh at all.10²10³10⁴10⁻¹⁹10⁻¹⁸10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²steps taken‖RᵀR − AᵀA‖ ⁄ ‖AᵀA‖the bound, linear in the steps√k · uevery step safe, the chain notdrift after the run2.5·10⁻¹⁵the bound there3.3·10⁻¹³√k · u there6.1·10⁻¹⁵worst single amplification2.7worst leverage met0.69refactorisations30backward stable onceand three thousand times is a different claim
Fig. 9 The sawtooth at a hundred steps, where each tooth is cut off before it has climbed far.
How far a carried triangular factor drifts from the data, over 2976 steps of a sliding windowA window of 24 rows on 6 columns, moved one row at a time: each step folds a row in with Givens rotations and removes one with hyperbolic rotations, and the factor is never rebuilt from the rows except every 25 steps. The rows are integers times powers of two, so AᵀA is exact in a double and the vertical axis is a distance from the answer. No single step amplifies by more than 2.72, no downdate fails, and after 2976 steps the factor is 8.72·10⁻¹⁹ from the matrix it is supposed to factor, against 3.87·10⁻¹⁴ with no refresh at all.10²10³10⁴10⁻²²10⁻²⁰10⁻¹⁸10⁻¹⁶10⁻¹⁴10⁻¹²steps taken‖RᵀR − AᵀA‖ ⁄ ‖AᵀA‖the bound, linear in the steps√k · uevery step safe, the chain notdrift after the run8.7·10⁻¹⁹the bound there3.3·10⁻¹³√k · u there6.1·10⁻¹⁵worst single amplification2.7worst leverage met0.69refactorisations120backward stable onceand three thousand times is a different claim
Fig. 10 And at twenty-five, where the drift never leaves the unit roundoff at all.

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.

The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻¹⁴ and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 7.3·10¹⁰ apart, with the reported one at 6.92·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.01428425670849811210⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.1·10⁻¹⁰unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10¹³iterations drawn110the recurrence remembers every roundingand the stopping test is written in it
Fig. 11 The same distinction inside a single solve, from the essay on what a method reports: a residual carried by a recurrence, drifting from the residual of the answer it belongs to.
Iterative refinement of the inverse-and-multiply solve at κ = 10^8The backward error starts at 3.3·10⁻¹¹ and falls by about κu a step — 1.8·10⁻¹⁷, 2.7·10⁻¹⁷, 2.4·10⁻¹⁷, 2.1·10⁻¹⁷, 2.1·10⁻¹⁷ — reaching the LU route's 3.1·10⁻¹⁷ after 5 corrections, each costing 2n² flops against the 2n³ the inversion cost. The forward error, drawn above, does not improve: 1.6·10⁻⁸ to 1.2·10⁻⁹, against the LU route's 9.6·10⁻¹⁰. What refinement at working precision buys is stability, and the accuracy floor belongs to the problem.01234510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1correction steprelative errorLU route: η = 3.1·10⁻¹⁷LU route: forward 9.6·10⁻¹⁰forward errorbackward errorwhat a correction buysη before refinement3.3·10⁻¹¹η after four steps2.1·10⁻¹⁷forward, unchanged1.2·10⁻⁹cost of a step, flops2048the residual is repairableand the accuracy floor is the problem's
Fig. 12 And the opposite: an outer loop that asks the matrix again at every step, and forgives an inner object that is much worse than anything here.

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.

Downdating a Cholesky factor against refactorising it, as the removed row's leverage approaches oneA rank-one term is removed from a 8×8 Gram matrix by hyperbolic rotations and, separately, by factorising the downdated matrix from scratch. The downdate's residual rises from 1.8·10⁻¹⁶ at h = 0.3 to 4.1·10⁻¹⁰ at h = 1 − 10⁻⁷, a slope of 1.03 against 1/(1 − h). The refactorisation is flat at 1.3·10⁻¹⁶, and the matrix both are producing has a condition number of 10.3 at every point on the axis — so the difficulty belongs to the route and not to the answer.110¹10²10³10⁴10⁵10⁶10⁷10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸1/(1 − h), the leverage of the removed row‖R̄ᵀR̄ − (G − aaᵀ)‖ / ‖G − aaᵀ‖downdatedrefactorisedat h = 1 − 10⁻⁷κ of the downdated matrix10κ of the matrix downdated3.6·10⁷rotation's amplification421downdate residual4.1·10⁻¹⁰a hyperbolic rotation is not orthogonaland that is exactly what it is for
Fig. 13 A narrower window, where the leverage rises and the per-step danger with it. That is the other axis, and it is not the one this essay is on.
Sherman–Morrison against a direct solve, on 20×20 systems whose answer is the identityA is QΛQᵀ with Λ = (1, …, 1, ε) and the rank-one update takes ε back to 1, so A + uvᵀ is the identity and κ of the problem being solved is 1.000000 at every point on the axis. A direct solve returns 1.1·10⁻¹⁶. The update formula, which is exact algebra, returns 2.5·10⁻⁴ — a slope of 1.00 against κ(A), the matrix that was replaced. Its arithmetic is three solves against that matrix and there is nowhere else the error could have come from.10¹10³10⁵10⁷10⁹10¹¹10¹³10¹⁵10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A), the matrix that was updatedrelative forward errorSherman–Morrisondirect solve of A + uvᵀone answer, two routesκ of the answer's matrix1κ of the matrix replaced10·10¹³update formula's error2.5·10⁻⁴direct solve's error1.1·10⁻¹⁶the question's condition number is 1at every point on this axis
Fig. 14 And the standing warning that a cheap update to a hard problem inherits the hard problem — which is about conditioning, and is a third axis again.

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.

How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 15 The amplifier, applied to a drift rather than to a measurement error. It does not care where the perturbation came from.
Two perturbation bounds and the error that was measured, on a 8×8 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 4·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 16 And the sharper version of the same bound, from the field that measured it: what the amplification is entry by entry rather than in a norm.

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.

How far |r_nn| sits above σ_min on Kahan's matrix, against the size and the parameter3 curves of |r_nn| ÷ σ_min against n, one per Kahan parameter. Every curve rises without turning over, reaching 7·10⁴ at n = 30, c = 0.5. Column pivoting makes no interchange at any point on any of them, so the failure is not a poor choice — there is nothing to choose.813182328110¹10²10³10⁴10⁵10⁶size of the matrix|r_nn| ÷ σ_minthe two agreec = 0.2c = 0.35c = 0.5no ceilingratio at n = 1021ratio at n = 307·10⁴interchanges, anywhere0the greedy rule never had a choiceand the gap grows with every row
Fig. 17 The repair on the third of them, from the arithmetic field: compensated summation, which is the same trade — a fixed extra cost per step against an error that stops growing.
What residual replacement costs and what it buys, at three periodsEach row is the same conjugate gradient run with the recomputed residual assigned back into the recurrence every k steps. Without it, the reported residual reaches 6.92·10⁻²¹ and the answer's stalls at 5.07·10⁻¹⁰. Replacing every 5 steps costs 30 extra matrix–vector products on top of 150 — 20 per cent — and brings the answer's residual to 8.31·10⁻¹⁷, with the two residuals then agreeing to a factor of 2.8.never replace5.07·10⁻¹⁰replace every 251.5·10⁻¹⁶replace every 101.34·10⁻¹⁶replace every 58.31·10⁻¹⁷reported 6.92·10⁻²¹ · 0 extra productsreported 5.3·10⁻¹⁷ · 6 extra productsreported 1.08·10⁻¹⁶ · 15 extra productsreported 2.95·10⁻¹⁷ · 30 extra productsthe residual of the answer the run returnsone line, at three pricesnever: the answer's residual5.1·10⁻¹⁰never: what it reported6.9·10⁻²¹every 5: the answer's8.3·10⁻¹⁷every 5: what it reported3·10⁻¹⁷extra products for that30ask the matrix againand the recurrence forgets what it did
Fig. 18 And on the second, where the extra cost is a matrix–vector product every few steps.

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.

How far a carried triangular factor drifts from the data, over 2976 steps of a sliding windowA window of 24 rows on 6 columns, moved one row at a time: each step folds a row in with Givens rotations and removes one with hyperbolic rotations, and the factor is never rebuilt from the rows. The rows are integers times powers of two, so AᵀA is exact in a double and the vertical axis is a distance from the answer. No single step amplifies by more than 2.72, no downdate fails, and after 2976 steps the factor is 3.87·10⁻¹⁴ from the matrix it is supposed to factor — a fitted slope of 0.554 in the step count, against a bound whose slope is 1.10²10³10⁴10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²steps taken‖RᵀR − AᵀA‖ ⁄ ‖AᵀA‖the bound, linear in the steps√k · uevery step safe, the chain notdrift after the run3.9·10⁻¹⁴the bound there3.3·10⁻¹³√k · u there6.1·10⁻¹⁵worst single amplification2.7worst leverage met0.69refactorisations1backward stable onceand three thousand times is a different claim
Fig. 19 The run the claims are about, once more. Every one of the four numbers above is a measurement of this figure rather than a description of it.
Classical and modified Gram–Schmidt: the same subtraction, in a different orderTwo panels of pseudocode differing in one argument, with the resulting pairwise dot products of the computed columns listed beneath each.for each previous column i, subtract the projection of column j onto qᵢclassicalr[i][j] = qᵢ · a[j] ↑ the ORIGINAL columnv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 7 and 8: 1columns 6 and 8: 0.13columns 6 and 7: 0.13modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 8: 4.4·10⁻⁷columns 2 and 8: 2.7·10⁻⁷columns 3 and 8: 2.4·10⁻⁸The two R factors agree to 1.2·10⁻⁶ relative. The two Q factors do not.the 8×8 Hilbert matrixone word, eight orders
Fig. 20 Orthogonality lost by a chain of operations that are individually fine, from the essay that established the pattern for the collection.
‖QᵀQ − I‖ at every step of a Lanczos run, n = 40The departure from orthogonality of the Lanczos basis, against the step, on a logarithmic vertical axis. It sits at the level of rounding for 13 steps and then climbs by a factor of about seventeen a step for ten steps running, twelve orders of magnitude, before saturating. The dashed line marks the step at which the first Ritz value converged, which is step 11.051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹step‖QᵀQ − I‖first Ritz value converges√uno reorthogonalisationfullthe two events are one eventstep orthogonality crosses √u14step the first Ritz value converges11consecutive steps at ≥5× growth10a drift would grow like a square rootthis grows like a geometric series and then stops
Fig. 21 And the same collapse in a Krylov basis, where it is sudden rather than gradual — a third shape for a carried object to fail in.
What 24 members of a sequence cost, by what changes between themEvery bar is 24 solves of an 120×120 system, costed in multiplications with a dense Cholesky at n³/3 and a solve at 4n². When nothing changes the answer is reused and the cost is one factorisation. When only the right-hand side changes, one factorisation serves every member — 12.7 per cent of the independent cost, with every solve still at a relative residual of 1.54·10⁻¹³. When the matrix drifts, the factorisation is reusable for a while and the while has to be measured: 7 factorisations for 24 members. When everything changes, nothing carries.nothing4.2%the right-hand side12.7%the matrix, slowly83.2%everything100.0%what can be reused: the answerwhat can be reused: the factorisationwhat can be reused: the factorisation, for a whilewhat can be reused: nothingshare of the cost of a sequence that shares nothingwhat 24 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations7independent: factorisations24what changes between the membersdecides what may be carried
Fig. 22 Where this sits in the field: it is not a reuse decision at all, but the reason the decisions in this field have to be made against something that is recomputed.
The growth factor of a 12×12 elimination, against the condition number of the matrixThree curves against κ. Cholesky's growth factor is exactly 1 at every condition number drawn — the elimination never produces an entry larger than the matrix already had. Partial pivoting on the same matrices reaches the same growth, and reaches it by making up to 9 row interchanges, each of which destroys the symmetry that was the reason to use a symmetric factorisation. The bound the general theory allows, 2^11 = 2048, is drawn above them both.10¹10³10⁵10⁷10⁹10¹¹110¹10²10³10⁴condition number of the matrixgrowth factorbound 2^11partial pivotingCholeskyno pivot to gain fromCholesky growth, every κ1Cholesky interchanges0partial pivoting, worst9the bound, 2^112048both eliminations reach the same growthand only one of them had to swap to get there
Fig. 23 A Cholesky factorisation’s growth, from the essay that measured it. Nothing in this run is a growth problem, which is what makes the drift the only thing left to blame.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 24 And the two errors, once more: the drift lands in the first and the conditioning turns it into the second.

Where this goes next

10¹10³10⁵10⁷10⁹10¹¹10¹³10¹⁵10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A), the matrix that was updatedrelative forward errorSherman–Morrisondirect solve of A + uvᵀone answer, two routesκ of the answer's matrix1κ of the matrix replaced10·10¹³update formula's error2.5·10⁻⁴direct solve's error1.1·10⁻¹⁶the question's condition number is 1at every point on this axis A correction cheaper than the problem Sherman and Morrison's formula updates a solved system for a rank-one change to the matrix, at 4n² operations instead of (2/3)n³. It is exact algebra. On a problem whose updated matrix is the identity — condition number one, the easiest system there is — it returns a forward error of 2.5·10⁻⁴ where a direct solve returns 10⁻¹⁶. 110¹10²10³10⁴10⁵10⁶10⁷10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸1/(1 − h), the leverage of the removed row‖R̄ᵀR̄ − (G − aaᵀ)‖ / ‖G − aaᵀ‖downdatedrefactorisedat h = 1 − 10⁻⁷κ of the downdated matrix4.3κ of the matrix downdated9.3·10⁶rotation's amplification344downdate residual3.5·10⁻¹⁰a hyperbolic rotation is not orthogonaland that is exactly what it is for The observation that cannot be removed Removing a rank-one term from a Cholesky factor needs a rotation that is not orthogonal, and the number under its square root is 1 − h, where h is the leverage of the row being removed. The algorithm's breakdown condition and the statistician's warning are the same quantity, arrived at from opposite ends, and neither field states it in the other's language. 10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸ The order they are added in Addition is associative in the algebra and is not associative in the arithmetic. The same million numbers, added in a different order, give answers that differ in the third significant figure — and the fix is not a wider float, it is a different order.

Reads more easily once this is understood

Essays that name this one as worth reading first.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Backward errorCholesky factorisationExact ground truthGivens rotationHyperbolic rotationLeverageLow-rank updateRandom walkRecursive least squares