A correction cheaper than the problem

A matrix changes a little. One row of a design matrix arrives; one entry of a stiffness matrix is edited; one constraint is switched on; one observation is dropped. The factorisation already in hand cost (2/3)n³ operations and the change is rank one, so the incentive not to do it again is enormous — and there is a formula for exactly that.

(A + uvᵀ)⁻¹  =  A⁻¹ − (A⁻¹u)(vᵀA⁻¹) / (1 + vᵀA⁻¹u)

Sherman and Morrison’s identity — an equation whose unknown is a matrix has the same shape one dimension up. In the form anybody should use — never forming an inverse, per the previous essay — it is two triangular solve pairs against the factorisation already in hand, one inner product and one axpy: 4n² against (2/3)n³, a factor of n/6, which at n = 1,000 is 167.

It is an identity. It is exact. And it charges for the shortcut in a currency nobody looks at.

The construction, which removes every excuse #

The temptation with a result like this is to blame the test problem, so the test problem is chosen to leave nowhere to stand.

A is QΛQᵀ with Λ = (1, 1, …, 1, ε): an orthogonal matrix, a spectrum of ones, and one eigenvalue at ε. Its condition number is 1/ε. The rank-one update is exactly the correction that takes that last eigenvalue back to one — u = (1−ε)q and v = q, with q the offending eigenvector — so

A + uvᵀ  =  QQᵀ  =  I

The updated matrix is the identity. Its condition number is 1.0000000000000029. There is no easier linear system in existence; the solution is the right-hand side.

A direct solve of it returns a forward error of 1.1·10⁻¹⁶. Sherman–Morrison, at ε = 10⁻¹⁴, returns 2.5·10⁻⁴.

The slope names the culprit #

The forward error of the update formula, fitted against κ(A) across four decades, has a slope of 1.00. Not against κ of the answer, which is 1 throughout; against the condition number of the matrix the formula routed its arithmetic through.

That is the whole explanation and it is one sentence long. Every quantity in the identity — A⁻¹b, A⁻¹u, vᵀA⁻¹b, vᵀA⁻¹u — is a solve against A. The arithmetic never touches A + uvᵀ. So the errors it makes are A’s errors, and A is the matrix nobody is asking about.

The identity is not wrong. Nothing about it is even slightly approximate. What it is not is a computation of the object it names, and the distinction between those two things is the difference between algebra and numerical analysis.

The slope is the claim, so it is the thing to check at more than one size.

At n = 6, 12, 20, 30 and 40 the fitted slope reads 1.005, 0.986, 0.995, 1.004 and 1.035. That is the exponent of a mechanism rather than a coincidence of one size: the error is κ(A)·u, to the first digit, at every size the generator will draw, on a family whose answer is the identity throughout.

The control, without which this is a trick #

A construction chosen to break a formula will break it, and the finding is worth nothing without the other half.

Take an ordinary A with κ = 10 and an ordinary rank-one update, and Sherman–Morrison returns 6.7·10⁻¹⁴ against a direct solve’s 4.9·10⁻¹⁴ — within a factor of 1.4, on a problem whose own condition number is 2,885. The formula is fine. It is worth a factor of n/6, it is exact algebra, and on a well-conditioned base it costs nothing at all.

So the rule is not do not use Sherman–Morrison. The rule is that its accuracy is a statement about the base matrix and not about the answer, and a reader who checks κ of the thing they wanted has checked the wrong number.

That is a sharper and more useful conclusion than a blanket prohibition, and it comes with a decision procedure: where a direct solve with A would have been acceptable, the update is safe. Where A was being updated because it had gone bad, it is not.

And the repair is one more application of the same identity #

The decision procedure above is sound and it gives up more than it needs to, because the defect it is avoiding is repairable by the cheapest means this collection has.

Refine with the updated matrix. Compute r = b − (A + uvᵀ)x̂, which is one matrix–vector product, and solve for the correction with the same Sherman–Morrison identity. On the construction above:

at ε = 10⁻⁶ the error goes 3.2·10⁻¹² → 0, exactly, in one step. At 10⁻⁸, 3.4·10⁻¹⁰ → 0. At 10⁻¹⁰, 2.9·10⁻⁸ → 6.5·10⁻¹⁴ → 0. At 10⁻¹², 4.7·10⁻⁶ → 1.3·10⁻⁹ → 2.9·10⁻¹³ → 0. At 10⁻¹⁴, 2.5·10⁻⁴ → 3.1·10⁻⁶ → 1.1·10⁻⁸ → 1.2·10⁻¹⁰, still falling.

The contraction per step is κ(A)·u — about a hundredfold at ε = 10⁻¹⁴, where κ(A)u is 10⁻² — so the iteration converges whenever κ(A) < 1/u, which is every base matrix that has a usable factorisation at all. The number of steps needed is one per decade of κ(A) beyond the accuracy wanted, and that is a quantity a code can estimate.

And it stays cheap, which is the whole reason the shortcut exists. A refinement step is a matrix–vector product with the updated matrix plus another 4n² identity solve. At n = 1,000, four refinements cost about 2·10⁷ multiplications against the direct route’s 6.7·10⁸ — thirty times cheaper, against the unrefined shortcut’s one hundred and sixty-seven.

Which changes the rule, and says why the other repair failed #

One caution on the stopping rule, because the residual it reads is computed through the updated matrix and the correction is computed through the base. The residual is trustworthy — it is a matrix–vector product with A + uvᵀ and a subtraction, and nothing in it goes near A⁻¹. The correction is not, and inherits κ(A) exactly as the first solve did. That is why the iteration converges at a rate rather than in one step: each pass reduces the error by 1/(κ(A)u) and no further, because the corrector is as inaccurate as the original solve was.

The consequence for a code is that the residual is the right thing to watch and the correction’s size is not. A correction that has stopped shrinking means the iteration has reached its floor; a correction that is small on the first step means nothing, because a small correction computed through an ill-conditioned route is a small wrong number.

So the rule is narrower than the section above states. It is not use the update where a direct solve with A would have been acceptable. It is:

Use the update, refine it, and stop when the correction stops shrinking. The residual against the updated matrix is computable, it costs one product, and it is the number that says whether the answer is any good.

That restores the thing the identity had given away — a quantity a caller can check — and it does so without giving back the factor of n/6 that made the shortcut worth having.

It is worth putting beside the explicit inverse, where refinement is measured on the same kind of defect and does not fully repair it: there four steps take the backward error from 4.5·10⁻⁵ to 3.1·10⁻¹⁷ and leave the forward error at 2.5·10⁻³, wandering around κu and staying there.

Same repair, two outcomes, and what decides between them is whose conditioning caused the trouble. For the explicit inverse the answer is ill-conditioned — κ = 10¹⁴ is a property of the problem — and no refinement at working precision sees past it. Here the answer is the identity’s, its condition number is 1.0000000000000029, and the only ill-conditioned object in the room is the route the algorithm chose. Refinement repairs a bad route and cannot repair a bad question, and the distinction is exactly the one the whole collection is organised around.

Where the shortcut is actually invoked #

It is worth being concrete about who reaches for this formula, because the diagnosis above is only useful if a reader can tell which case they are in.

Adding an observation to a fit. A new row arrives and the Gram matrix AᵀA gains aaᵀ. Here the base matrix is the Gram matrix of everything seen so far, it is generally the better conditioned of the two, and the update is safe. This is the benign case and it is the most common one.

Switching a constraint on. An active-set method adds a constraint and the KKT matrix gains a rank-one term. Here the base is a matrix that was already at the boundary of what the method could handle — that is why the constraint was activated — and the update is exactly the dangerous case.

A leave-one-out or a sliding window. A row is removed, which the formula handles by taking u to be its negative. The base is the full Gram matrix, which is fine; the answer is the reduced one, which may not be. The next essay is entirely about this case, and finds that the number governing it is one a statistician already plots.

A rank-one edit to a physical model. A spring stiffness changes, a boundary condition switches. The base is the assembled stiffness matrix, which is usually badly conditioned for reasons that have nothing to do with the edit — a graded mesh, a penalty term, mixed units — and the update inherits all of it.

Three of those four are cases where somebody would happily solve with A. One is not, and it is the one where the incentive to avoid a refactorisation is strongest.

Woodbury, which is the same sentence at rank k #

The rank-one identity generalises:

(A + UCVᵀ)⁻¹  =  A⁻¹ − A⁻¹U (C⁻¹ + VᵀA⁻¹U)⁻¹ VᵀA⁻¹

and the generalisation is where it earns its keep, because the inverse in the middle is k×k. A Kalman filter’s measurement update, a Gaussian process with k inducing points, a domain decomposition with k interface unknowns — all of them are this formula, and all of them are worth the factor of (n/k)³ it buys.

And all of them inherit the same defect, from the same place: every term is a solve against A. The extra term makes it worse rather than better, because C⁻¹ + VᵀA⁻¹U is a k×k matrix assembled from quantities that already carry A’s error, and its own conditioning multiplies on top.

The Kalman filter’s practitioners found this the hard way and their answer is instructive: the square-root filters propagate a Cholesky factor rather than the covariance, so the object that could lose definiteness is never formed. It is the same move again — do not build the object — reached independently in a field that had no reason to phrase it that way.

What the denominator is #

The formula has one obviously dangerous place: 1 + vᵀA⁻¹u, in a denominator. It is zero exactly when A + uvᵀ is singular, which is honest — a formula for an inverse ought to break where the inverse does not exist — and near-zero when the update nearly destroys invertibility.

That case is well known and it is not the case this essay is about. It is worth separating them carefully, because conflating them lets the formula off.

When the denominator is small, the answer is genuinely sensitive: A + uvᵀ is nearly singular, κ of the question is large, and a direct solve would be in trouble too. That is conditioning and nobody is to blame.

When the denominator is perfectly healthy — as it is throughout the construction above, where the updated matrix is the identity — and the answer is still wrong by ten orders of magnitude, that is the algorithm. Checking the denominator is a real safeguard against a real failure and it says nothing at all about this one.

What it is worth, counted #

The saving is real and it should be stated as precisely as the cost, because the conclusion is a trade rather than a prohibition.

At n = 20 the update is 1,600 operations against 5,333 — a factor of 3.3. At n = 100 it is 40,000 against 667,000, a factor of 17. At n = 1,000 it is a factor of 167, and at that size a refactorisation is 667 megaflops that somebody is going to notice.

The far end of the sweep is where the trade is most favourable and the accuracy is worst, and the near end is where there is no trade at all.

So the trade has a size below which there is nothing on one side of it. The saving is n/6, so it is 1.0 at n = 6, 2.0 at 12, 3.3 at 20, 5.0 at 30 and 6.7 at 40 — while the update’s error over the same five sizes runs 3·10⁻⁴, 3.1·10⁻⁴, 2.5·10⁻⁴, 8.1·10⁻⁴ and 3.1·10⁻³ against a direct solve that stays between 3.4·10⁻¹⁷ and 1.9·10⁻¹⁶ throughout. Both columns move the same way, which is the uncomfortable shape: the formula is worth more exactly where it costs more, and worth nothing at the size where it costs least.

Against that: a forward error of κ(A)·u instead of κ(A + uvᵀ)·u. If the two condition numbers are comparable — which is the usual case, since most updates are small perturbations — the trade is free. If they are not, the trade is a factor of n/6 in time against a factor of κ(A)/κ(A + uvᵀ) in accuracy, and that second factor is unbounded.

There is a middle option that is usually available and rarely taken: update, then refine. One step of iterative refinement against the updated matrix costs 2n² — the same order as the update itself — and it is a residual computed against A + uvᵀ, which is the matrix nobody’s arithmetic had touched. The previous essay measured what refinement does to a route whose backward error is κu: it takes it down by about κu a step. Here it does the same thing, and the cost is a constant factor on a shortcut that was already a factor of n/6 ahead.

That the option exists and is not standard is worth noticing. The reason is that the update formula is usually reached for in a context — a filter, an inner loop, an active-set iteration — where the result is immediately consumed rather than inspected, so nothing ever computes the residual that would show the problem.

The shape, again #

This is the fourth essay in this phase to reach the same sentence from a different algorithm, and the second to reach the second sentence.

The first: the object a question names — det A, A⁻¹, (A + uvᵀ)⁻¹ — is usually not the object worth computing.

The second, which this essay states most sharply: a shortcut’s error is set by the object it routed through and not by the problem it answers. The next essay finds the identical statement in an algorithm that shares no arithmetic with this one, where the object routed through is a Cholesky factor and the number that decides everything is a statistic every regression package already prints.

The filter that had to stop forming its own covariance #

The clearest case of an update formula being reached for under pressure, and of the field working out what to do about it, is the Kalman filter — and the answer it reached is this phase’s sentence arrived at independently.

The filter carries a covariance P and updates it at every measurement. The textbook form is P ← (I − KH)P with K the gain, which is a rank-k correction and is exactly a Woodbury identity written out. It is also numerically catastrophic in a specific way: P has to stay symmetric positive semidefinite, the update is a difference of two positive quantities, and in floating point the difference can come out with a negative eigenvalue. A covariance with a negative eigenvalue is not a covariance, and the filter that produced it will happily keep running and produce a gain built from an imaginary standard deviation.

The repairs are a catalogue of increasing seriousness. Joseph’s form rewrites the same update as (I − KH)P(I − KH)ᵀ + KRKᵀ, which is a sum of two positive semidefinite terms and cannot go negative — more arithmetic, and structurally safe. Symmetrising each step, P ← (P + Pᵀ)/2, removes the drift that the asymmetry of the arithmetic introduces. And the square-root filters — Potter’s, Bierman and Thornton’s U-D factorisation, and the modern QR-based ones — never form P at all: they propagate a factor S with SSᵀ = P, so a matrix that could lose definiteness does not exist anywhere in the computation.

That last is the phase’s sentence, reached in the 1960s by people who had no interest in stating it generally: the object the equations name is a covariance, and the object worth computing is a factor of it. The condition number of the factor is the square root of the covariance’s, which is the same saving the least-squares field gets from QR over the normal equations, and it arrives here for the same reason — a factor is one square root closer to the data than the matrix it multiplies out to.

What is worth carrying #

An update formula is exact algebra and its arithmetic goes through the old matrix. Sherman–Morrison computes four quantities and all four are solves against A, so its error is κ(A)u whatever the updated matrix looks like — measured here at a slope of 1.00 against κ(A) on a problem whose own condition number is one.

Check the base, not the answer. Where a direct solve with A would have been acceptable, the update is safe. Where A is the matrix being updated because it had gone bad, refactorise.

And the denominator is a different warning. 1 + vᵀA⁻¹u near zero means the problem is nearly singular; it is a genuine safeguard against a genuine failure and it does not catch this one.

The next essay removes a row instead of adding one, and finds the same mechanism in an algorithm that has no formula in it at all: the observation that cannot be removed.