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

A factorisation kept past its date

One Cholesky factor can serve five members of a drifting sequence and save 44 per cent of the work. Kept for twenty it does not lose accuracy — it stops converging altogether. The optimum and the cliff are four members apart, both move with the drift, and a rule written in a ratio the iteration has already computed finds them without being told what the drift is.
29 min read 23 figures Where the algebra stops being the arithmeticMeasured, not assertedExact ground truthAssertions that reject

Worth reading first: The problem that arrives again · A factorisation with nothing to pivot for.

A continuation solves the same problem two hundred times with a parameter moved slightly between each one. A time-stepping code does the same thing with time as the parameter. In both, the matrix at member t is nearly the matrix at member t − 1, and the expensive object in the room is a factorisation of it.

Reusing that factorisation is the oldest trick in the subject. Newton’s method with a frozen Jacobian is the chord method, it is in every textbook, and what every textbook says about it is that it converges linearly rather than quadratically. That is true and it is not the interesting part. The interesting part is that a factorisation has a shelf life: a number of members it is worth serving, past which it should be thrown away, and the number is not a property of the algorithm.

The cost of a drifting sequence of nonlinear solves, against how often the Jacobian is refactorisedTwenty members, each warm-started from the last, each solved by a chord iteration on a factorisation that may be several members old. Refactorising at every member costs 16.42 MFlop; refactorising every 5 costs 9.27. Past a period of 20 the chord iteration stops converging altogether, which is drawn as an open circle on the ceiling rather than omitted. The filled square is the rule that refactorises when the observed contraction ratio exceeds 0.2: 4 factorisations, 9.22 MFlop, and it was never told the drift rate.024681012141618202210⁶10⁷10⁸members served by one factorisationmultiplications for the whole sequencedoes not convergethe contraction ruledrift 0.01 a memberevery member1.6·10⁷every 5 members9.3·10⁶contraction rule9.2·10⁶its factorisations4cliff at a period of20a factorisation has a shelf lifeand the cliff is past the optimum
Fig. 1 Twenty members of a drifting sequence, walked ten times with the factorisation rebuilt at ten different periods. The open circles on the ceiling are the periods at which the iteration stops converging.

What is measured

The family is a continuation: F_t(x) = Ax + λx³ − b_t, with the root of member t known exactly because b_t was built from it. Each member is warm-started from the previous member’s answer, which is what makes the whole thing worth doing, and solved by a chord iteration — the same triangular solves against the same stored factor, repeated until the residual is below 10⁻¹⁰.

The cost model is a count and not an estimate: a factorisation is n³/3 + n² multiplications and a chord step is 4n² — a residual evaluation and a pair of triangular solves. At n = 120 that makes a factorisation worth about 30 chord steps, which is the exchange rate the whole essay turns on.

Chord steps per member of a drifting sequence, with the factorisation rebuilt every 4Each point is one member of a continuation, warm-started from the previous member's answer and solved with a factorisation that is between zero and 3 members old. The lower series is the same sequence with a fresh factorisation at every member, which costs 6 to 5 steps throughout. The kept run costs 6 steps on a fresh factor and 12 at its worst, for 5 factorisations against 20.03691215182103691215member of the sequencechord steps to reach the tolerancea fresh factorisation each memberthe ringed points are the refactorisationsrebuilt every 4factorisations5chord steps160worst member's steps12multiplications1.2·10⁷with a fresh factor each time1.8·10⁷a factor is cheapest the member it was built forand dearest the member before it is replaced
Fig. 2 The same sequence seen per member: how many chord steps each one needs, with the ringed points marking the members where the factorisation was rebuilt.

The sawtooth

The per-member cost is the shape worth having, because it says what staleness costs in the currency the iteration is actually spending.

At a drift of 0.01 a member with the factor rebuilt at every member, every member takes 4 chord steps, twenty times over. Rebuilt every four, the members take

4 5 6 7 · 4 5 6 7 · 4 5 6 7 · 4 5 6 7 · 4 5 6 7

which is the sawtooth: the member the factor was built for takes four steps and the member three later takes seven. Rebuilt every eight, the teeth are longer and steeper —

4 5 6 7 8 9 11 12 · 4 5 6 7 8 9 10 11 · 4 5 6 7

— and the last member under one factor costs three times what the first did.

Nothing here is a loss of accuracy. Every one of those members reaches the same tolerance and returns an answer at the same distance from a root that is known exactly. What ageing costs is steps, and steps are cheap, which is exactly why the trade is worth making.

Chord steps per member of a drifting sequence, with the factorisation rebuilt every 1Each point is one member of a continuation, warm-started from the previous member's answer and solved with a factorisation that is between zero and 0 members old. The lower series is the same sequence with a fresh factorisation at every member, which costs 6 to 5 steps throughout. The kept run costs 6 steps on a fresh factor and 6 at its worst, for 20 factorisations against 20.0369121518210246810member of the sequencechord steps to reach the tolerancea fresh factorisation each memberthe ringed points are the refactorisationsrebuilt every 1factorisations20chord steps102worst member's steps6multiplications1.8·10⁷with a fresh factor each time1.8·10⁷a factor is cheapest the member it was built forand dearest the member before it is replaced
Fig. 3 The flat case, where the factor is rebuilt at every member and every member costs the same four steps. All of the cost is in the factorisations.
Chord steps per member of a drifting sequence, with the factorisation rebuilt every 8Each point is one member of a continuation, warm-started from the previous member's answer and solved with a factorisation that is between zero and 7 members old. The lower series is the same sequence with a fresh factorisation at every member, which costs 6 to 5 steps throughout. The kept run costs 6 steps on a fresh factor and 39 at its worst, for 3 factorisations against 20.036912151821071421283542member of the sequencechord steps to reach the tolerancea fresh factorisation each memberthe ringed points are the refactorisationsrebuilt every 8factorisations3chord steps272worst member's steps39multiplications1.7·10⁷with a fresh factor each time1.8·10⁷a factor is cheapest the member it was built forand dearest the member before it is replaced
Fig. 4 And the long teeth, where one factor serves eight members and the last of them costs three times the first.

Where the optimum is

Summing the sawtooth against the factorisations gives the crossover. At a drift of 0.01 a member, over twenty members:

rebuild every factorisations chord steps MFlop
1 20 80 16.42
2 10 116 11.09
3 7 138 9.84
4 5 110 9.29
5 4 129 9.27
6 4 155 9.62
8 3 144 10.07
10 2 168 11.26
15 2 219 14.72
20 1 — does not converge

The optimum is every five members at 9.27 MFlop against 16.42 for refactorising every time — 44 per cent of the work saved, on a decision that is one integer in a configuration file.

And then the cliff. At a period of fifteen the sequence still converges, at 59 per cent more work than the optimum, which is a bad choice and an available one. At twenty it does not converge at all.

The cost of a drifting sequence of nonlinear solves, against how often the Jacobian is refactorisedTwenty members, each warm-started from the last, each solved by a chord iteration on a factorisation that may be several members old. Refactorising at every member costs 14.11 MFlop; refactorising every 20 costs 4.85. Past a period of the range drawn the chord iteration stops converging altogether, which is drawn as an open circle on the ceiling rather than omitted. The filled square is the rule that refactorises when the observed contraction ratio exceeds 0.2: 1 factorisations, 4.85 MFlop, and it was never told the drift rate.024681012141618202210⁶10⁷10⁸members served by one factorisationmultiplications for the whole sequencethe contraction ruledrift 0.001 a memberevery member1.4·10⁷every 20 members4.9·10⁶contraction rule4.9·10⁶its factorisations1cliff at a period of∞a factorisation has a shelf lifeand the cliff is past the optimum
Fig. 5 At a tenth of the drift the curve flattens and the cliff moves out past the range drawn: with a sufficiently slow sequence, one factorisation serves all of it.

The cliff is not a gradual thing

The refusal this essay carries is aimed at the most natural wrong reading of the table above, which is that keeping a factorisation too long buys work with accuracy.

It does not. Every period from 1 to 15 returns the same answer — the same distance from the known root, to the tolerance the iteration was asked for — and differs only in what it cost. At 20 the chord iteration’s contraction factor passes 1 and the iteration diverges: the residual grows, the run is cut off, and there is no answer of any quality at all.

That is a different shape of failure from everything else in this field. A loose forcing term costs work. A stale preconditioner costs iterations. A stale factorisation costs steps until it costs everything, and there is no warning between the last good period and the first bad one that anything was about to change.

The cost of a drifting sequence of nonlinear solves, against how often the Jacobian is refactorisedTwenty members, each warm-started from the last, each solved by a chord iteration on a factorisation that may be several members old. Refactorising at every member costs 19.01 MFlop; refactorising every 3 costs 14.50. Past a period of 6 the chord iteration stops converging altogether, which is drawn as an open circle on the ceiling rather than omitted. The filled square is the rule that refactorises when the observed contraction ratio exceeds 0.2: 8 factorisations, 12.33 MFlop, and it was never told the drift rate.024681012141618202210⁶10⁷10⁸members served by one factorisationmultiplications for the whole sequencedoes not convergethe contraction ruledrift 0.032 a memberevery member1.9·10⁷every 3 members1.5·10⁷contraction rule1.2·10⁷its factorisations8cliff at a period of6a factorisation has a shelf lifeand the cliff is past the optimum
Fig. 6 At three times the drift the cliff has walked in to a period of four, and the optimum with it. The distance between them is a member or two.
Jacobi, Gauss–Seidel and SOR at ω = 1.000A semi-logarithmic plot of relative residual against iteration for three stationary methods, with dashed reference curves showing the rate each is predicted to contract at.02675348011068133510⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹iteration‖r‖ / ‖b‖JacobiGauss–SeidelSOR ω=1.00closed form vs measuredρ Jacobi, exact1ρ measured1ρ Gauss–Seidel, exact0.991D Laplacian, n = 32ω optimal = 1.826
Fig. 7 The mechanism, from the field that measured it: a stationary iteration converges when its contraction factor is below one and diverges when it is above, and there is nothing in between.

The shelf life shortens with the drift

The optimum is not a property of the method or of the matrix. It is a property of how fast the sequence is moving, and it moves the way it should:

drift a member best period its cost cliff at refactorise every member
0.01 5 9.27 MFlop 20 16.42
0.02 3 12.14 10 17.68
0.05 2 16.85 4 20.22
0.10 1 23.16 2 23.16

At the fastest drift the optimum and the cliff are adjacent: refactorising at every member is the only fixed policy that works, and it is also the best one. At the slowest, one factor serves five members and the failure is four periods further out.

The practical difficulty is now visible. The right period depends on a quantity the code does not know, and getting it wrong in one direction wastes 77 per cent of the work while getting it wrong in the other produces no answer at all.

The cost of a drifting sequence of nonlinear solves, against how often the Jacobian is refactorisedTwenty members, each warm-started from the last, each solved by a chord iteration on a factorisation that may be several members old. Refactorising at every member costs 15.26 MFlop; refactorising every 10 costs 6.71. Past a period of the range drawn the chord iteration stops converging altogether, which is drawn as an open circle on the ceiling rather than omitted. The filled square is the rule that refactorises when the observed contraction ratio exceeds 0.2: 2 factorisations, 7.46 MFlop, and it was never told the drift rate.024681012141618202210⁶10⁷10⁸members served by one factorisationmultiplications for the whole sequencethe contraction ruledrift 0.0032 a memberevery member1.5·10⁷every 10 members6.7·10⁶contraction rule7.5·10⁶its factorisations2cliff at a period of∞a factorisation has a shelf lifeand the cliff is past the optimum
Fig. 8 The middle of the range, where the optimum sits at four and the cliff at ten. Every drift rate has its own version of this figure and nothing in a single member says which one applies.

The rule that reads the ratio

The chord iteration already computes the quantity that decides this, and does not use it.

At each step it evaluates the residual, so it has ‖F_k‖ and ‖F_{k−1}‖. Their ratio is the observed contraction factor of the iteration as it is actually running — not a bound, not an estimate from the matrix, the thing itself. A fresh factorisation gives a ratio around 0.05 on this problem; a factor several members old gives 0.3; a factor about to fail gives one close to 1.

So: refactorise the moment the observed ratio exceeds ρ*, and use nothing else.

At ρ* = 0.2 the rule costs 9.22 MFlop at a drift of 0.01 with four factorisations, against the best fixed period’s 9.27 with four. At 0.02 it costs 10.80 with six, against the best fixed period’s 12.14. At 0.05 it costs 13.26 with nine, against 16.85. At 0.10 — where every fixed period but one has fallen off the cliff — it costs 15.51 with thirteen factorisations, against 23.16 for the only fixed policy that works, a third less.

It beats the best fixed period at every drift rate measured, by more as the drift gets faster, and it never fails.

Total inner iterations for a whole inexact Newton solve, against the constant forcing termEvery point is a complete solve of the same problem to the same outer tolerance of 10^-10. At η = 10⁻¹⁴ it takes 9358 conjugate gradient iterations across 9 Newton steps; at η = 0.1 it takes 980 across 10. The left arm is oversolving and the right arm is too many outer steps. The open circle is the adaptive rule, which is not a constant: it costs 1009 iterations and reaches 3.37·10⁻¹¹.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110²10³10⁴10⁵constant forcing term ηtotal inner iterationsthe adaptive ruleoversolvingstarving the outer loopone problem, twelve budgetsη = 10⁻¹⁴, inner iterations9358cheapest constant0.1its inner iterations980adaptive rule, iterations1009adaptive final residual3.4·10⁻¹¹outer tolerance asked for10⁻¹⁰a tolerance is a cost decisionand its optimum is not machine precision
Fig. 9 The same shape one object earlier: a rule that reads the residual beating every constant, on the tolerance rather than on the factorisation.
What each Newton step asked of its linear solve, and what the solve cost, under the adaptive policyThe lower series is the tolerance handed to the inner solve at each step and the upper one is the number of conjugate gradient iterations it took. Under the adaptive rule the first step asks for 0.9 and costs 1 iteration, and the last asks for 0.0042 and costs 271. The whole solve costs 1009 inner iterations across 10 Newton steps and ends at a relative residual of 3.37·10⁻¹¹.0123456789101110⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110³Newton steptolerance asked for, and iterations paiditerations paidtolerance asked forouter residualthe adaptive policy, step by stepNewton steps10inner iterations, total1009first step's cost1last step's cost271final outer residual3.4·10⁻¹¹the rule reads the last two residualsand asks for nothing it cannot use
Fig. 10 And that rule’s own trail, for comparison. Both rules read a ratio of consecutive residuals; only the quantity they set differs.

Why it wins where it wins

The reason the rule does better as the drift gets faster is worth being explicit about, because it is not that it is cleverer there.

A fixed period is a claim that every factorisation is worth the same number of members. It is not: the drift is not uniform along the sequence, the nonlinearity means the Jacobian changes faster at some members than at others, and a factorisation built at an easy member genuinely lasts longer than one built at a hard member. A fixed period has to be short enough for the worst factorisation in the run, so it wastes work on all the others.

The rule spends each factorisation until it is used up. At the fastest drift, where the variation between members is largest, that is worth a third of the run; at the slowest, where every factorisation is worth about five members, it is worth almost nothing, and the two policies land within half a per cent of each other.

A rule is worth what the variation is worth, which is a sentence about when to bother rather than about what is optimal.

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. 11 A cheap correction to a hard problem, from the essay that priced it. A stale factorisation is the same bargain run backwards: the correction is free and the problem is what it always was.
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⁵00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10000worst found 76206×6, 200 directionsmedian reaches 0.29 of κ
Fig. 12 And the amplifier that does not move. A factorisation four members old solves the same problem as a fresh one; what changes is how many steps it takes to get there.

Two ways to be wrong, and they are not the same size

The asymmetry in this decision is worth stating as a number, because it is what makes the rule worth having rather than merely elegant.

At a drift of 0.01 a member, being too eager — refactorising every member when every five would do — costs 16.42 MFlop against 9.27, which is 77 per cent more work for the same answer. Annoying, invisible, and survivable for years.

Being too patient — keeping a factor for twenty members when five was right — costs the run. Not an answer with fewer digits: no answer. The residual grows, the iteration is cut off, and whatever the outer code does with a member that failed to converge is now the behaviour of the program.

Those two are three periods apart at one drift rate and adjacent at another. A configuration file holding a number tuned on last month’s problem is therefore sitting somewhere on a line whose safe end and cliff have both moved since it was written, and nothing about a successful run reports how close to the edge it was.

The contraction rule removes the asymmetry rather than balancing it. It cannot be too patient, because the quantity it watches is exactly the one that says the iteration is in trouble, and it notices within a single step. It can be slightly too eager — ρ* = 0.2 is a choice — and being too eager is the cheap direction.

The cost of a drifting sequence of nonlinear solves, against how often the Jacobian is refactorisedTwenty members, each warm-started from the last, each solved by a chord iteration on a factorisation that may be several members old. Refactorising at every member costs 14.11 MFlop; refactorising every 20 costs 4.85. Past a period of the range drawn the chord iteration stops converging altogether, which is drawn as an open circle on the ceiling rather than omitted. The filled square is the rule that refactorises when the observed contraction ratio exceeds 0.2: 1 factorisations, 4.85 MFlop, and it was never told the drift rate.024681012141618202210⁶10⁷10⁸members served by one factorisationmultiplications for the whole sequencethe contraction ruledrift 0.001 a memberevery member1.4·10⁷every 20 members4.9·10⁶contraction rule4.9·10⁶its factorisations1cliff at a period of∞a factorisation has a shelf lifeand the cliff is past the optimum
Fig. 13 The slow case again, where the two ways of being wrong are fifteen members apart and a fixed period is a perfectly reasonable thing to write down.
Chord steps per member of a drifting sequence, with the factorisation rebuilt every 6Each point is one member of a continuation, warm-started from the previous member's answer and solved with a factorisation that is between zero and 5 members old. The lower series is the same sequence with a fresh factorisation at every member, which costs 6 to 5 steps throughout. The kept run costs 6 steps on a fresh factor and 20 at its worst, for 4 factorisations against 20.03691215182104812162024member of the sequencechord steps to reach the tolerancea fresh factorisation each memberthe ringed points are the refactorisationsrebuilt every 6factorisations4chord steps202worst member's steps20multiplications1.4·10⁷with a fresh factor each time1.8·10⁷a factor is cheapest the member it was built forand dearest the member before it is replaced
Fig. 14 And a period one past the optimum, which is not a failure and is not free either: the last member under each factor costs twice the first.
Jacobi sweeps on the 6×6 Hilbert matrixA plot of the off-diagonal norm against sweep number, falling from about one to ten to the minus seventeen in five sweeps, with the matrix shown at three stages beneath it.01234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweep‖off-diagonal‖eigenvalues, largest first1.6190.24240.016326.157·10⁻⁴1.257·10⁻⁵1.083·10⁻⁷κ = 1.5·10⁷‖VᵀV − I‖ = 1.9·10⁻¹⁵1.000.500.330.250.200.170.500.330.250.200.170.140.330.250.200.170.140.130.250.200.170.140.130.110.200.170.140.130.110.100.170.140.130.110.100.09as given1.620.050.00-0.040.00-0.000.050.240.000.010.00-0.000.000.000.000.00-0.00-0.00-0.040.010.000.020.000.000.000.00-0.000.000.00·-0.00-0.00-0.000.00·0.00after one sweep1.62······0.24······0.00······0.02······0.00······0.00after 8 sweepseach sweep is n(n−1)/2 rotationsthe trace is conserved
Fig. 15 A contraction factor measured rather than bounded, from the field that established the practice. The rule in this essay is that measurement used as a decision instead of as a report.

What the warm start is doing underneath

There is a quantity that does not appear anywhere above and without which none of it works, which is worth making explicit because it is the reason a factorisation can be stale at all.

Each member is started from the previous member’s answer. The correction it needs is therefore of the size of the drift rather than of the size of the solution, and a chord iteration converges linearly with a factor governed by how far the frozen Jacobian is from the current one. Those two facts multiply: a small correction contracted by a factor near 0.05 reaches the tolerance in four steps, and the same factor applied to a cold start would need dozens.

So the shelf life measured here is not a property of the factorisation alone. It is a property of the factorisation and the warm start together, and a code that solves each member from scratch would find a much shorter one — the frozen Jacobian would have to carry the iteration all the way from nothing rather than across a gap the continuation already closed.

That is the sense in which this field’s four decisions are not independent of each other. Warm starting makes reuse worth having; reuse makes the per-member cost small enough that the drift rate becomes the binding constraint; and the drift rate is what decides the period. Change any one of the three and the arithmetic above has to be redone, which is the argument for measuring rather than inheriting a number from somebody else’s problem.

What is not measured here

Two limits, both real.

The cost model is dense. A factorisation at n³/3 and a solve at 4n² is right for a dense Cholesky and wrong for everything a large code actually holds. The essay on what a rebuild is worth takes the exchange rate out of the model and makes it the axis, because it is the only thing that moves the answer.

The problem is symmetric positive definite, so there is no pivoting to reuse and no question of whether the factorisation is still valid — only of whether it is still good. The essay on a pivot order kept is the other case, and there the reuse can fail in a way that no number of extra steps repairs.

The cheapest rebuild period for one drifting sequence, against what a rebuild costsOne sequence, one drift rate, one preconditioner — and six answers, because the answer is not a property of the sequence. Costed in iteration-equivalents, the optimal period runs from every 1 member at a setup worth 5 iterations to every 12 at a setup worth 200. A dense Cholesky at this size is worth 6.7 iterations, which is at the left of the axis, so for it the answer is always to rebuild. The rule that rebuilds when a solve takes 1.5 times what the last fresh one did beats the best fixed period at 1 of the 6 ratios and loses at the rest, because it cannot see the setup cost at all.10¹10²10³03691215what one rebuild is worth, in preconditioned iterationscheapest number of members between rebuildsa dense Cholesky here is worth 6.7 iterationsringed: where the growth rule beats every fixed periodcheapest period, by setup costsetup worth 51setup worth 102setup worth 203setup worth 506setup worth 1008setup worth 20012how long to keep itis a question about what it cost to build
Fig. 16 The exchange rate as an axis, in the essay that measures it. Everything on this page is one column of that figure.
The backward error of a solve along a sequence, with the pivot order chosen fresh, kept, and kept after equilibrationA 10×10 conflict grid, 100 unknowns, with three of its rows scaled down through 9 decades as the sequence runs — a row that was ordinary becoming a row that is small in the matrix's own units. The sparsity pattern is identical at every member, so a symbolic phase computed once stays valid throughout. Choosing a fresh order each time holds the backward error at the working precision. Keeping the first member's order costs 2 replaced pivots and a backward error of 4.82·10⁻⁹. Keeping the same order after dividing each row by its largest entry costs nothing at all: 6.55·10⁻¹⁷, with no pivot replaced anywhere in the run. The ringed points on the upper curve are the members at which another pivot fell below the floor and was replaced.02468101210⁻¹⁸10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸member of the sequencebackward error of the solvefresh order · equilibrated kept orderringed: another pivot replaced hereone order, two sets of unitsfresh order, last member9.5·10⁻¹⁷kept order, last member4.8·10⁻⁹equilibrated, last member6.5·10⁻¹⁷pivots replaced, kept2pivots replaced, equilibrated0the perturbation floor √u1.1·10⁻⁸a reused pivot order is safe or notdepending on what the rows are measured in
Fig. 17 And the reuse that can be invalid rather than merely stale.
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. 18 Why there is no pivoting question here: a symmetric positive definite matrix has nothing to pivot for, from the essay that established it.
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. 19 And where this sits in the field’s ledger: the drifting bar, which is the one whose height is a decision rather than an arithmetic fact.
An incomplete Cholesky kept while the operator turns anisotropic, against one rebuilt at every memberThe five-point operator's two directions stop being equally weighted, ε running from 1 down to 0.0331, with the sparsity pattern unchanged throughout — so the factorisation from the first member stays applicable for the whole run, which is the situation in which it gets kept. The kept one goes from 18 iterations to 52. The rebuilt one goes from 18 to 13, because an anisotropic operator is an easier problem for a factorisation that knows about the anisotropy. The two start at the same point by construction and never meet again.0246810121416182001020304050member of the sequencepreconditioned conjugate gradient iterationskept from the first memberrebuilt every memberone operator, two policieskept, first member18kept, last member52rebuilt, first member18rebuilt, last member13ε at the last member0.033the problem got easierand the kept preconditioner got worse at it
Fig. 20 The next object along, whose ageing is gentler and much harder to see.
Preconditioned iterations against the size of the drift, for a drift on each half of the spectrumThe preconditioner is the exact Cholesky factorisation of the first matrix, so it solves that matrix in one step and every iteration counted here is bought by the drift. Both curves are drifts of the same relative Frobenius norm; the upper one sits on the half of the spectrum with the small eigenvalues and the lower on the half with the large. At a relative drift of 10⁻² they cost 38 and 9 iterations, because the drift is divided by the eigenvalue it lands on: the preconditioned matrix has eigenvalues 1 + w/λ, and λ runs over 10⁴.10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10112233445566relative drift ‖E‖ / ‖A₀‖preconditioned conjugate gradient iterationson the small eigenvalueson the large onessame drift, two placessmall end at 10⁻²38large end at 10⁻²9κ(M⁻¹A), small end36κ(M⁻¹A), large end1.3‖E‖/‖A₀‖, both0.01how much the matrix changedis not what the preconditioner cares about
Fig. 21 And the quantity that decides how fast it ages, which is not the one a rule of thumb would use.
Growth factor under partial pivoting to n = 20: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.048121620110¹10²10³10⁴10⁵10⁶10⁷10⁸matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 20: bound 5.2·10⁵, worst 2.8
Fig. 22 The growth factor of an elimination, from the field that measured it. Nothing in this essay’s sequence moves it, because a symmetric positive definite matrix has nothing to pivot for.
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. 23 And the habit the whole measurement rests on: an answer that is known rather than approximated, so “the same answer” is a statement rather than a comparison of two computations.

Where this goes next

nothing8.3%the right-hand side16.5%the matrix, slowly85.9%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 12 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations4independent: factorisations12what changes between the membersdecides what may be carried The problem that arrives again A hundred and thirty essays have solved a system once and measured how wrong the answer was. Almost no computation is shaped like that. A solve is one step of an outer loop, its answer is an input rather than a deliverable, and four quantities this collection has treated as accuracy requirements turn out to be assets with a shelf life. 10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110⁻⁴10⁻³10⁻²10⁻¹1inner tolerance η, relative residual of the linear solvedistance from the root after the stepd² = 0.00138distance before the step, 0.03721093547231126the number by each point is the iterations it costeleven decades, one landing placedistance before the step0.037its square0.0014where η = 10⁻³ lands0.0025where η = 10⁻¹⁴ lands0.0025iterations for the first354iterations for the second1126the accuracy that is thrown awaymeasured against a root that is known The accuracy that is thrown away A Newton step is the exact answer to a linearised problem, and the linearisation is wrong at second order. So there is a floor under how close the step can land, the floor is the square of where it started, and eleven decades of inner tolerance below it buy the same four digits at four times the price. 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⁻¹⁶.

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.

Cholesky factorisationChord methodContraction factorExact ground truthFlop countLinear convergenceNewton iterationTriangular solveWarm start