Numerical Linear Algebra
where the algebra stops being the arithmetic
Least squares, and the road not to take

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⁻¹⁶.
22 min read 16 figures Whose fault is itIdentical algebra, different arithmeticMeasured, not assertedAssertions that reject

Worth reading first: The inverse that is never formed · The condition number is an amplifier · The exact answer to a nearby 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. 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.

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. 1 Sherman–Morrison against a direct solve, on systems whose updated matrix is the identity. κ of the problem being solved is 1 at every point on the axis; the horizontal axis is the condition number of the matrix that was replaced. Drag the size.

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

Sherman–Morrison against a direct solve, on 40×40 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.9·10⁻¹⁶. The update formula, which is exact algebra, returns 0.0031 — a slope of 1.04 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 error0.0031direct solve's error1.9·10⁻¹⁶the question's condition number is 1at every point on this axis
Fig. 2 The same construction at n = 40, where the shortcut is worth a factor of 6.7 in arithmetic and loses the same digits. The size changes what the formula saves and not what it costs.

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.

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. 3 What a condition number amplifies, and what the slope of one above is measuring: the update’s answer is a function of quantities computed against a matrix whose amplification factor is 1/ε.
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: κ · u30×30, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 4 The two errors again. Here the backward error is with respect to a system the arithmetic never saw, which is why the honest measurement is the forward error against a solution that is an input.

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.

The condition number a sketched preconditioner leaves, against the condition number it was givenTwo curves against κ(A), both axes logarithmic. The matrix's own condition number climbs the diagonal from 100 to 10¹⁰; κ(AR⁻¹), where R comes from a QR of a 4n-row sketch, is 2.2284 at every one of them — the same number to ten digits, not a similar one. The reason is two lines of algebra: with G = SU the preconditioned singular values are those of (GᵀG)⁻¹, which has no spectrum of A in it at all.10²10⁴10⁶10⁸10¹⁰110²10⁴10⁶10⁸10¹⁰condition number of the matrixcondition number seen by the iterationκ(A), unchangedκ(AR⁻¹)a bound with no κ(A) in itκ(AR⁻¹), every κ(A)2.2κ(SU), the other route2.2κ(A), across the sweep10⁸κ(AR⁻¹), across the sweep1the sketch never sees the spectrumand the spectrum cancels out of the answer
Fig. 5 The general shape of the question — which matrix’s condition number governs an answer. Here the arithmetic goes through one matrix and the question is about another, and the first wins.

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.

Two penalties on the same problem, against the offset in the signalBest relative error for each penalty at four offsets. With no offset the two are within 7% of each other. At an offset of 10 the derivative penalty is 1.85 times better, because a constant lies in its null space and costs it nothing, while the norm penalty pays for the whole offset at every λ.best relative error at each offset‖x‖, offset 00.2035‖L₁x‖, offset 00.1901‖x‖, offset 20.0572‖L₁x‖, offset 20.0524‖x‖, offset 50.0333‖L₁x‖, offset 50.0245‖x‖, offset 100.0240‖L₁x‖, offset 100.0129what the null space buysadvantage at offset 01.1advantage at offset 21.1advantage at offset 51.4advantage at offset 101.9the norm penalty pays for a constantthe derivative penalty does not
Fig. 6 The constrained case, where the base matrix is badly conditioned by construction: a penalty method multiplies one block by a large number on purpose, and a rank-one edit to that system routes its arithmetic through the large number.
The same first pivot decision, with the corner entry at 10⁻¹²Three sparsity patterns side by side: the matrix, the factor obtained by eliminating the corner first, and the factor obtained by pivoting on the largest entry instead.the matrix105 entriescorner first — sparsest227 entries, growth 1.9·10¹¹largest first — safe242 entries, growth 1.19the middle factor is the smaller one, and its answer has no correct digitsboth factorisations reproduce the matrix‖PA − LU‖/‖A‖, sparsest3.8·10⁻¹⁷‖PA − LU‖/‖A‖, pivoted5.4·10⁻¹⁷forward error, sparsest3·10⁻⁵forward error, pivoted4.8·10⁻¹⁶red marks are entries elimination createdthe fill argument and the stability argument disagree
Fig. 7 And the assembled physical model, whose conditioning comes from the modelling rather than from the edit. An update to it is a solve against it.

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.

The arrowhead matrix, eliminated from each endThree sparsity plots. The first shows an arrowhead matrix with a dense first row and column. The second shows its Cholesky factor, completely dense. The third shows the factor obtained after moving the dense row to the end, which has no fill at all.the matrix43 entriestip eliminated first253 entriestip eliminated last43 entries‖A − LLᵀ‖/‖A‖, tip first1.4·10⁻¹⁶‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 253 · sparse factor is 2n − 1 = 43one row swapped to the endnothing numerical chose between them
Fig. 8 The pattern a Woodbury update usually has: a large easy block and a small dense correction. Solving it as a whole and solving it as a correction are two routes to one answer, and they are not equally accurate.
Two-level convergence factor for three interpolations, at ε = 0.001Three bars per operator: the classical formula, the unconstrained energy minimiser, and the minimiser constrained to reproduce a constant. On the isotropic and aligned operators the constraint costs a factor of 4.8; on the rotated one it buys a factor of 1.33.lower is better; each bar is the factor the error falls by per two-level cycleisotropic, classical0.0609isotropic, minimiser0.0609isotropic, constrained0.0848aligned, classical0.0612aligned, minimiser0.0615aligned, constrained0.2933rotated 45°, classical0.3464rotated 45°, minimiser0.2843rotated 45°, constrained0.2131what the constraint is worthconstraint at isotropic1.4constraint at aligned4.8constraint at rotated 45°0.75the constraint costs where the method worksand buys where it does not
Fig. 9 And the same shape as a constrained problem, where the correction is the constraint. The choice between assembling the whole system and correcting a solved one is the same choice.

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.

The widest box a uniqueness proof is available forVerified half-width against the distance between two roots, both axes logarithmic. The measured widths are 0.05, 0.005, 5·10⁻⁴, 5·10⁻⁵ at separations of 0.1, 0.01, 0.001, 10⁻⁴ — a ratio of 0.5000 at every one of them. The operator verifies exactly up to the point where the box would reach the second root, and a box containing both is refused.10⁻⁴10⁻³10⁻²10⁻¹10⁻⁴10⁻³10⁻²10⁻¹distance between the two rootswidest verified half-widththe separationverified widthuniqueness is a claim about a regionwidth ÷ separation at δ = 0.10.5width ÷ separation at δ = 0.010.5width ÷ separation at δ = 0.0010.5verified up to the second rootand refused past it
Fig. 10 The genuinely ill-posed case, where a range of answers fits equally well and no arithmetic can choose. That is the failure a small denominator warns about, and it is a different failure from the one measured here.

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.

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.

The same arithmetic, twice, with M = 144 words of fast memoryThree counters against the matrix size on logarithmic axes. The operation count is one curve, because both orderings perform exactly the same operations — 338,120 of them at n = 80. The words moved between fast and slow memory are two curves a factor of 3.24 apart. The factorisations they produce are identical to the last bit.10²10³10⁴10⁵matrix size ncountoperations, bothwords, unblockedwords, blocked (b = 6)the answer does not move‖PA − LU‖/‖A‖, unblocked2.8·10⁻¹⁶‖PA − LU‖/‖A‖, blocked2.8·10⁻¹⁶difference between them0the dashed curve is both orderings' operation countthe solid pair is what they cost
Fig. 11 Counting what an algorithm costs, which is the discipline this trade needs on both sides. An operation count is a fact; the accuracy it buys or spends is a measurement.

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.

Iterative refinement of the inverse-and-multiply solve at κ = 10^10The backward error starts at 6.4·10⁻⁹ and falls by about κu a step — 1.5·10⁻¹⁶, 2.2·10⁻¹⁷, 2.4·10⁻¹⁷, 2.4·10⁻¹⁷, 2·10⁻¹⁷ — reaching the LU route's 2.9·10⁻¹⁷ after 5 corrections, each costing 2n² flops against the 2n³ the inversion cost. The forward error, drawn above, does not improve: 1.1·10⁻⁶ to 4.9·10⁻⁸, against the LU route's 2·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: η = 2.9·10⁻¹⁷LU route: forward 2·10⁻⁸forward errorbackward errorwhat a correction buysη before refinement6.4·10⁻⁹η after four steps2·10⁻¹⁷forward, unchanged4.9·10⁻⁸cost of a step, flops1800the residual is repairableand the accuracy floor is the problem's
Fig. 12 What a correction step against the right matrix buys, from the previous essay. The same arithmetic applies to an updated system, and for the same reason.
Products with A to return all 2 copies, against the block widthOne bar per block width. Widths below 2 return fewer copies than the eigenvalue has, at every step count up to a basis of half the problem's dimension, and are marked rather than drawn. The cheapest width that works is 2, at 32 products with A; the widest drawn costs 84.an eigenvalue repeated 2 times, in a 200-dimensional problemblock of 1never returns them allblock of 232 products, 16 stepsblock of 351 products, 17 stepsblock of 460 products, 15 stepsblock of 570 products, 14 stepsblock of 684 products, 14 stepsthe narrowest that workscheapest width2its products with A32the widest drawn84narrower than the multiplicity is not slowit is absent
Fig. 13 And the general form of a saving that is bought with a structural choice. The question is always what the structure charges as well as what it saves.

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.

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. 14 The other instance. A hyperbolic downdate’s residual against the leverage of the row it removes, with a refactorisation flat beneath it and the answer’s condition number constant at four.
Backward and forward error of two routes to x, on 30×30 systems across eight decades of κBoth routes start from the same LU factorisation. The LU solve's backward error is 2.2·10⁻¹⁷ at every conditioning — a flat line at the unit roundoff. Multiplying by the explicitly formed inverse gives 4.5·10⁻⁵ at κ = 10¹⁴, a slope of 0.94 against κ. The two forward errors, drawn above them, are 2.8·10⁻⁴ and 0.015 — within a factor of 54, which is why the difference between the two methods is invisible to anyone measuring how wrong the answer is.10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the backward error
Fig. 15 And the instance from the previous essay, where the object routed through was an explicitly formed inverse.

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.

The arrowhead matrix, eliminated from each endThree sparsity plots. The first shows an arrowhead matrix with a dense first row and column. The second shows its Cholesky factor, completely dense. The third shows the factor obtained after moving the dense row to the end, which has no fill at all.the matrix15 entriestip eliminated first36 entriestip eliminated last15 entries‖A − LLᵀ‖/‖A‖, tip first3.1·10⁻¹⁷‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 36 · sparse factor is 2n − 1 = 15one row swapped to the endnothing numerical chose between them
Fig. 16 The two ways of solving a system with a small dense correction. The choice between assembling the whole thing and correcting a solved one is the same choice a filter makes at every measurement.

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.

Where this goes next

10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the backward error The inverse that is never formed x = A⁻¹b is how the solution of a linear system is written and it is not how it is computed. The usual reason given is cost — three times the arithmetic. The real reason is that one of the two routes is backward stable and the other is not, and at κ = 10¹⁴ they differ by twelve orders of magnitude in the number that says whose fault a wrong answer is. 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. 0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it The units the matrix is measured in One linear system, written twice. The rows of the second are the rows of the first in different units, the solution is identical to the last bit, and the condition number has moved by eight orders of magnitude. One of those two numbers is a fact about the problem and the other is a fact about the notation.

What links here

Computed from the collection, not written here: the essays that point at this one.

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.

Condition numberFlop countForward errorLeverageLow-rank updateLU factorisationMatrix inverse