The observation that cannot be removed

Adding a rank-one term to a Cholesky factor is easy and safe. A sequence of Givens rotations — orthogonal, each of them, and therefore each of them incapable of amplifying anything — takes R to the factor of RᵀR + xxᵀ in 2p² operations, and the result is as accurate as a fresh factorisation.

Removing one is a different problem, and the difference is not a matter of degree.

An orthogonal transformation preserves the Frobenius norm. Removing a term makes the matrix smaller. So no sequence of orthogonal transformations can perform a downdate, and the algorithm that does is built from rotations that are not orthogonal:

[ c  −s ]        with  c² − s² = 1   rather than   c² + s² = 1
[ −s  c ]

a hyperbolic rotation. It preserves a² − b² instead of a² + b², which is exactly what removing a term needs. And the property that made the orthogonal version safe — that ‖G‖ = 1, so nothing it touches can grow — is precisely the property it has given up. c and s are both unbounded.

The measurement #

A 40×6 design matrix, with the leverage of its first row set rather than measured — the row is scaled until h takes a chosen value, so the axis is an input. The Gram matrix G = AᵀA is factorised; then aaᵀ is removed from it two ways.

The fitted slope of the downdate against 1/(1 − h) is 1.04.

And the last column is the finding. The matrix being produced has a condition number of about four at every row of that table. It is not a hard matrix. A fresh factorisation of it is exact to rounding at every leverage. The difficulty belongs entirely to the route.

The number under the square root #

The mechanism is one line of algebra and it is the reason this essay exists.

Downdating G by aaᵀ begins by solving Rᵀp = a. The first hyperbolic rotation’s coefficients are 1/√(1 − pᵀp), and

pᵀp  =  ‖R⁻ᵀa‖²  =  aᵀ(RᵀR)⁻¹a  =  aᵀ(AᵀA)⁻¹a  =  h

which is the leverage of the row a — the diagonal of the hat matrix H = A(AᵀA)⁻¹Aᵀ, a quantity every regression package computes and every regression textbook plots.

So the algorithm’s amplification factor is 1/√(1 − h), its breakdown condition is h ≥ 1, and both are the statistician’s warning written in a different alphabet. A leverage near one means this observation is the only thing determining some direction of the fit; the downdate says the same thing by having no real square root.

Neither field states it in the other’s language. Numerical analysis texts derive the downdate and observe that it fails when the result is not positive definite. Statistics texts define leverage, prove it sums to p, and warn about points near one. They are the same number and the two derivations never meet.

The statistic that is the same downdate #

Leave-one-out cross-validation asks: refit without observation i, and how badly does the fit predict it? Doing that honestly is m refits. The closed form everybody uses instead is

e₍ᵢ₎  =  eᵢ / (1 − hᵢ)

— one division per observation, and it is a rank-one downdate written as arithmetic. Every package that reports a studentised or deleted residual is running this.

It agrees with refitting. Measured against forty actual refits on a forty-observation problem, the worst relative disagreement is 1.3·10⁻¹². That agreement is what licenses the shortcut, and it is worth having as a measurement rather than as a derivation, because the derivation is exact and the question is arithmetic.

And the divisor is 1 − h. At h = 0.999 the closed form multiplies the residual by a thousand, and the thousand is correct — that observation really does drag the fit onto itself, and removing it really does move the prediction that far. What is not correct at that point is any digit past the fourth, because the numerator was computed against a fit that observation dominated.

Why the closed form is everywhere and the downdate is not #

There is an asymmetry between the two halves of this essay worth accounting for. The leave-one-out formula is in every statistics package; the hyperbolic downdate is in LINPACK, in a handful of signal processing libraries, and almost nowhere else.

The reason is that the two are doing different amounts of work. eᵢ/(1 − hᵢ) produces one number — what the fit would have predicted at observation i without it — and needs the leverages, which cost one triangular solve each and are usually wanted anyway. The downdate produces a whole factor, which is what a caller needs when the reduced problem is going to be solved again: a sliding window, a recursive least squares with a forgetting factor, an active-set method releasing a constraint.

So the closed form is the cheap special case of the general operation, and it is the one that has a denominator a user can look at. The downdate has the same denominator buried inside a square root inside a loop, where nothing reports it.

That is worth naming as a general habit rather than a fact about these two routines. A shortcut whose dangerous quantity appears in the output is safer than one whose dangerous quantity appears in an intermediate, even when the arithmetic is identical — because a caller can see the first and cannot see the second. eᵢ/(1 − hᵢ) at h = 0.999 announces itself by returning a residual a thousand times the others. A downdate at h = 0.999 returns a factor that looks like every other factor.

The three quantities that agree, and one that does not #

Three routes to h, sharing no arithmetic:

The triangular solve. ‖R⁻ᵀa‖², from a forward substitution against the Cholesky factor. The inverse. aᵀG⁻¹a, forming G⁻¹ explicitly — one of the legitimate uses of an inverse, since the entries are the answer, and one where forming a Gram matrix has already been paid for. The hat matrix. The (i, i) entry of A(AᵀA)⁻¹Aᵀ, from the projection itself.

They agree to 10⁻⁸ at every leverage tested, which is the site’s two routes habit applied to a quantity that has three.

And one that does not agree, which is the sum rule: Σhᵢ = p, exactly. Not approximately — the trace of a projection onto a p-dimensional subspace is p, and the measurement returns 5.999999999999998 for a six-column fit. Adding a column adds exactly one unit of leverage to be shared out, whatever the data looks like, so the average leverage is p/m and depends on nothing.

That last fact is worth pausing on because it says what leverage is. The response y does not appear in it anywhere. Leverage is a property of the design — of where the observations were taken — and not of what was measured there.

When the downdate simply fails, which it does not #

The obvious expectation is that past a leverage the arithmetic can distinguish from one, a step’s discriminant a² − b² comes out negative and the algorithm stops — an honest outcome, better than a wrong answer, and the behaviour a caller can act on.

It does not happen. Swept from h = 0.9 down to 1 − h = 10⁻¹⁵, the discriminant stays positive at every step at every leverage, and the routine reports success every time:

The last row is a factor wrong in the second digit, returned without complaint, on a matrix a fresh Cholesky handles to 10⁻¹⁶. There is no refusal anywhere on the sweep, so there is no left-hand half of the refuse-or-lie picture to be in.

The mechanism is that the breakdown test is local and the amplification is global. The discriminant at step k compares two scalars, R[k][k] against x[k]. The quantity that decides the accuracy is 1 − ‖R⁻ᵀa‖², a sum over all p components. A vector can have norm within 10⁻¹⁵ of one with every individual component comfortably inside its own diagonal — so every step passes its own test, each contributes a modest amplification, and the product reaches 10⁷.

So the danger is not a band just before a refusal. It is twelve decades of leverage, from about 1 − h = 10⁻⁴ down, in which the routine answers, the answer degrades by one digit per decade, and nothing in the output distinguishes any of it from the h = 0.3 case. That is the essay’s own better sentence, four sections earlier, arriving as a measurement: a downdate at h = 0.999 returns a factor that looks like every other factor, and it turns out to be true all the way down.

There is one genuine boundary and it belongs to the problem. At 1 − h = 1.1·10⁻¹⁶ the refactorisation fails too, because G − aaᵀ is by then numerically indefinite and the removal really is not well-posed. Everything above it is a matrix with a condition number of four that one route returns to rounding and the other does not.

assertTheDowndateDoesNotRefuseItAnswersWrongly holds the whole sweep: that every leverage reports success, that the fresh factorisation is exact at every one of them, that the condition number stays a single digit throughout, that the worst residual is ten orders above the route it replaces, and that the growth is the 1/(1 − h) the first rotation predicts. The one thing it cannot hold is a refusal, because there is none to hold.

What follows for the guard #

The measurement changes what a caller has to do, and makes the advice at the end of this essay load-bearing rather than prudent.

If the routine refused near h = 1, a caller could treat the refusal as the guard: attempt the downdate, catch the failure, refactorise. That is the shape most numerical software is written in, and it is why the discriminant test exists at all. It does not work here, because the failure never arrives — so the leverage has to be computed before the downdate, not inferred from it, and the threshold has to be chosen rather than discovered.

Choosing it is straightforward once the table above exists. The downdate loses about one digit per decade of 1 − h, so a caller wanting k digits of the reduced factor needs 1 − h above 10^(k−16). Ten digits is 1 − h > 10⁻⁶, which on a forty-observation six-column fit is almost every row; the guard rejects the handful that are nearly determining the fit by themselves, and for those the honest answer is the refactorisation the essay recommends. One triangular solve buys the number that decides.

That is a better arrangement than a refusal would have been, and it is worth saying so rather than treating it as a repair. A refusal is a verdict computed by the routine on its own terms, with a threshold nobody chose; a leverage is a number the caller can look at and set a policy against, in the same units the statistician already uses. The defect is not that the guard is missing. It is that the routine looks as though it has one.

What to do instead, and what it costs #

The alternative is a refactorisation, and the cost is worth being exact about because it is smaller than it looks.

Downdating a p×p factor is 2p² operations. Refactorising G − aaᵀ from the Gram matrix is p³/3 — at p = 6 that is 72 against 72, which is to say nothing at all. The saving only becomes real when the Gram matrix itself has to be reassembled from a large A, at mp² operations — the same accounting a constraint imposed as a weight has to do when it changes the design, and that is a saving of a factor of m rather than of p.

So for a small model on a large dataset — the usual case in a sliding-window or leave-one-out setting — the downdate is worth having, and it is worth having with a guard: compute the leverage first. It costs one triangular solve, it is a number worth looking at for its own sake, and it decides whether the shortcut is safe before the shortcut is taken.

Where the leverage is high, the honest options are to refactorise, or to notice that the answer being asked for is a fit with one direction supported by a single observation, and that the numerical question is the smaller of the two problems.

Refactorisation returning exactly zero at three columns is the cleanest possible statement of the comparison, and it is not a fluke of that width — every other width puts it at the rounding level:

The slope is the claim rather than the size of the error, so it is worth reading down the whole slider rather than off two frames. Two more widths, and the middle of the range:

Six is where the essay’s headline numbers come from, so the two widths either side of it are the ones that say whether the headline is about the method or about the width it was measured at:

The slope is one, and now it is one five times. 0.958, 1.125, 1.038, 1.022 and 0.963 — a mean of 1.02 with the extremes 12% either side, across a factor of three in the number of columns. The error’s growth is exactly first order in 1/(1 − h), which is the statement the caption above makes and which two frames could only suggest.

And the answer is well conditioned throughout, which is what makes the result a defect rather than a difficulty. κ of the answer reads 2.3, 3.5, 4.3, 4.9 and 6.7 — single figures at every width. The downdate loses six orders of magnitude on a problem whose conditioning would forgive it nothing, so the loss belongs to the algorithm and not to the question. The refactorisation on the same problems sits at 0 to 1.5·10⁻¹⁶.

The trade the essay is about turns out to be a fixed one. The downdate costs 18, 32, 72, 128 and 200 flops against 369, 661, 1,512, 2,731 and 4,333 — so it is 20.5, 20.7, 21.0, 21.3 and 21.7 times cheaper, rising very slowly. About twenty-one times cheaper at every width, and about six orders of magnitude worse near leverage one at every width. Neither half of that bargain moves with the model’s size, which is the useful form of the recommendation: a code choosing between them is choosing between two fixed numbers and can decide once.

The mechanism, twice #

This is the second algorithm in two essays whose error is governed by the object it routed through rather than by the problem it answers, and the two share no arithmetic.

Sherman–Morrison’s four quantities are solves against A, so its error is κ(A)u while the answer’s condition number is one. A hyperbolic downdate’s rotations are built from R, so its error is proportional to 1/(1 − h) — a property of the original Gram matrix and the row being removed — while the matrix it produces has a condition number of four.

One formula, one sweep of rotations; one governed by a condition number, one by a leverage. The general statement is the phase’s second finding: where a method reaches its answer through an intermediate, the intermediate’s conditioning is what it charges, and the answer’s conditioning is what a reader checks.

The sliding window, which is a downdate at every step #

The setting where downdating is not an occasional operation but the whole algorithm is recursive least squares, and it is worth describing because it turns the leverage question from a diagnostic into a design parameter.

A recursive least-squares filter maintains a fit over the most recent w observations. At each step one row arrives and one leaves: an update and a downdate, both rank one, both at 2p² operations against a refactorisation’s mp². For a filter running at kilohertz over a model of six parameters, that is the difference between running and not.

The exponentially weighted variant is the same thing done continuously. Instead of dropping a row, every existing row is multiplied by √λ with λ slightly below one — a forgetting factor — so the influence of old data decays geometrically. Written as a factorisation update it is R ← √λ R followed by an ordinary Givens update for the new row, and there is no downdate in it at all. The hyperbolic rotations are gone, the algorithm is orthogonal throughout, and the whole of this essay’s difficulty disappears.

What it costs is that the window is soft. The effective number of observations is 1/(1 − λ), which is a statement about a geometric weighting rather than a count, and a filter that must forget an exact number of samples cannot use it. So the choice between a sliding window and a forgetting factor is a choice between a hard boundary computed unstably and a soft one computed stably — which is a modelling decision being made by a numerical constraint, and worth recognising as one.

And the leverage is the design parameter. In a sliding window the row about to leave has a leverage determined by how much of the fit it was holding up, and a window short enough that individual rows carry high leverage is a window in which the downdate is unstable and the fit is dominated by individual observations. Those are the same statement. A filter whose downdate is in trouble is a filter whose window is too short for its model, and the numerical diagnostic and the statistical one agree because they are the same number.

What is worth carrying #

An orthogonal sweep cannot make a matrix smaller, and the sweep that can has no bound on its coefficients. That is not an accident of one algorithm; it is what removing a term costs, and every downdating method pays it in some form.

The breakdown condition is the leverage, exactly, which means the diagnostic already exists and is already plotted in a different field. Compute h before downdating; it is one triangular solve and it is the number that decides.

A dangerous quantity in the output beats the same quantity in an intermediate. The leave-one-out formula divides by 1 − h where a reader can see it; the downdate takes its square root inside a loop where nobody does. Same number, same danger, and only one of them announces itself.

And a leverage near one is not ill-conditioning of the answer. The downdated matrix has a condition number of four at every leverage measured here. What a high leverage says is that one observation was holding a direction of the fit up, which is a statement about the design and a much more interesting one than a numerical warning.

The next essay leaves updates behind and asks a question the phase has been circling: what is a function of a matrix, and is the definition a method? A function of a matrix is not a function of its entries.