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

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.
31 min read 25 figures Where the algebra stops being the arithmeticMeasured, not assertedExact ground truthTwo routes to a number

Worth reading first: The exact answer to a nearby problem · The condition number is an amplifier · The accuracy that is thrown away.

Every essay in this collection so far has solved one system. A matrix arrives, an answer is wanted, the answer comes back wrong by some amount, and the whole of the work is deciding whether the fault is the algorithm’s or the problem’s. That framing is what the collection is for and it has held up through a hundred and thirty essays.

It also describes almost nothing that anybody runs.

A solve is nearly always one step of an outer loop. A Newton iteration produces a linear system per step and each one is built out of the answer to the last. A time-stepping code produces one per step for a hundred thousand steps. A continuation walks a parameter and solves at each stop; an optimisation solves a subproblem per iteration; a design sweep solves the same structure under two hundred load cases. In every one of those the vector the solver returns is not the answer to anything. It is an input, and the next member of the sequence will be assembled from it and handed back.

What 12 members of a sequence cost, by what changes between themEvery bar is 12 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 — 16.5 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: 4 factorisations for 12 members. When everything changes, nothing carries.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
Fig. 1 Twelve members of a sequence, costed four ways. The bars differ only in what changes between one member and the next, and every one of them solves twelve systems of the same size.

That changes what the accuracy of a solve is for, and the change is larger than it sounds. This field is about the change, and this essay is about the four situations it comes in.

The classification, which is not a taxonomy for its own sake

Between one member of a sequence and the next, one of four things has happened.

Nothing changed. The same matrix, the same right-hand side. This is rarer than it sounds and it is not a joke: it is what a badly written outer loop does, and the reason to name it is that the cost of the whole sequence is then the cost of one member and the correct implementation is a cache.

The right-hand side changed. One matrix, many loads. One factorisation serves the entire sequence, and every further member costs a pair of triangular solves. On a 120×120 problem the twelve members together cost 16.5 per cent of what twelve independent problems cost, and the reused factorisation is still solving every one of them to a relative residual of 1.5·10⁻¹³ — reuse here costs nothing at all, and the only thing to get wrong is failing to do it.

The matrix drifted. This is the interesting case and it is most of them: a Jacobian at a new iterate, an operator at a new time step, a stiffness matrix at a slightly deformed configuration. The factorisation from an earlier member is not a factorisation of this one, and it is not useless either. It is an approximation, and everything in this field follows from asking how good an approximation and for how long.

Everything changed. Independent problems in a loop. Nothing carries and the sequence is a sequence only in the scheduler.

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. 2 The same four kinds at twice the length. The independent bar is the one that grows in proportion to the members; the others grow in proportion to how often something has to be rebuilt.

Read as a ladder, each level licenses exactly one more kind of reuse than the level below it. Read as a decision, each level has a different way of being got wrong, and the expensive mistakes are not symmetric: reusing where reuse was not licensed costs correctness, and declining to reuse where it was costs a constant factor that nothing will ever report.

The problem this field measures on

Every claim in this field is made against an answer that is known rather than approximated, which is this collection’s standing habit and is unusually cheap to arrange here.

The family is F(x) = Ax + λx³ − b, with A symmetric positive definite at a prescribed condition number and b constructed from a chosen x*. So the root is not estimated, it is an input. The Jacobian is A + 3λ·diag(x²), which is a formula rather than a difference quotient — the essay on operators with no entries measures what happens when it is not — and it is symmetric positive definite, so the inner solve can be conjugate gradients and the inner tolerance is a knob with a known iteration count behind it.

Error of a difference quotient for J(x)v against ε, on the Bratu problem at n = 64The Jacobian of this problem is a formula, so the error of each quotient is measured against a derivative that is exact rather than against a better quotient. The forward difference falls with a slope of 1.01 — first order — reaches 1.28·10⁻¹⁰ at ε = 10^-6, and rises again with a slope of -1.00 as cancellation takes over. The central difference falls with a slope of 2.00 and bottoms at 1.11·10⁻¹² for two residual evaluations instead of one. Neither gets near the unit roundoff at any ε.10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice
Fig. 3 The other way to get a Jacobian, priced elsewhere in the collection: a difference quotient, with truncation on one side and cancellation on the other and a floor in between.

One parameter moves the whole family, so a sequence is a continuation: member t has root t·x*, and how fast the members drift apart is a number the experiment sets rather than inherits.

What an outer loop forgives

Here is the sentence the field turns on, and it is worth stating before any of the measurements that support it.

An outer iteration recomputes its residual from the matrix at every step. So whatever the inner solve got wrong is not carried forward — it is measured again, from scratch, against the problem that is actually being solved.

That is why a Newton step can be solved to two digits and still produce the same next iterate as one solved to sixteen. It is why a factorisation several members old can drive an iteration to full accuracy. It is why a pivot order chosen for a different matrix, with pivots deliberately replaced by made-up numbers, can be repaired to the working precision in two steps.

And it is exactly why the things that are not recomputed are dangerous, which is the other half of this field and is measured in the essays on what a carried factor accumulates and on what a recurrence remembers.

Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 4 The distinction this collection is built on, in the figure that opened it: the backward error is what the algorithm did and the forward error is what the problem did with it. In a sequence, only the first of the two is recomputed at every step.

Four assets, and each of them ages

Once a computation is a sequence, four objects that had looked like properties of a solve turn into objects with a lifetime.

The inner tolerance. In a single solve, a tighter tolerance is never wrong; it is only wasteful. In a sequence it is not even wasteful in an interesting way — it is wasteful in a measurable way, because the next outer step throws away most of what it bought. Solving every step of one Newton run to 10⁻¹⁴ costs 9,358 conjugate gradient iterations. A rule that reads the last two residuals and asks for what it can use costs 1,009, for the same final residual.

The factorisation. Reusing one Cholesky factor across a drifting sequence saves a factor of 1.9 at its best fixed period and stops converging one period later. The optimum is a period, the failure past it is a cliff rather than a slope, and a rule written in the observed contraction ratio finds both without being told the drift rate.

The preconditioner. The same object one level less exact, and it fails differently: a stale preconditioner never stops converging, it just costs more iterations, so nothing announces the problem. What decides how fast it ages is not how much the matrix changed but where the change landed in the spectrum, and the two differ by the condition number.

The pivot order. The most valuable of the four, because what it saves is a symbolic phase rather than arithmetic, and the least forgiving: a reused pivot order is safe or useless depending on the units the rows are written in, and there iterative refinement — the repair every other essay here leans on — does not work.

Iterative refinement of the inverse-and-multiply solve at κ = 10^8The backward error starts at 3.3·10⁻¹¹ and falls by about κu a step — 1.8·10⁻¹⁷, 2.7·10⁻¹⁷, 2.4·10⁻¹⁷, 2.1·10⁻¹⁷, 2.1·10⁻¹⁷ — reaching the LU route's 3.1·10⁻¹⁷ after 5 corrections, each costing 2n² flops against the 2n³ the inversion cost. The forward error, drawn above, does not improve: 1.6·10⁻⁸ to 1.2·10⁻⁹, against the LU route's 9.6·10⁻¹⁰. What refinement at working precision buys is stability, and the accuracy floor belongs to the problem.01234510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1correction steprelative errorLU route: η = 3.1·10⁻¹⁷LU route: forward 9.6·10⁻¹⁰forward errorbackward errorwhat a correction buysη before refinement3.3·10⁻¹¹η after four steps2.1·10⁻¹⁷forward, unchanged1.2·10⁻⁹cost of a step, flops2048the residual is repairableand the accuracy floor is the problem's
Fig. 5 The repair the whole field runs on, in the essay that introduced it: a factorisation that is not good enough, an outer loop that recomputes the residual, and an answer at the working precision.
LSQR steps to a normal residual of 10⁻¹², with and without a sketched preconditionerTwo curves against κ(A). The unpreconditioned solve takes 15, 22, 32, 59, 150 steps as the conditioning worsens, and at the worst one it does not reach the tolerance inside 150. The preconditioned solve takes 10, 10, 10, 10, 10 — flat across eight decades. Both return the same vector, to the accuracy the problem allows: the preconditioned answer and the direct one differ by about κ(A)·u at every κ, 5.7·10⁻⁶ at the worst. The SAME sketch used as an answer instead of as a preconditioner is 4.5 away.10²10⁴10⁶10⁸10¹⁰110¹10²10³condition number of the matrixLSQR steps to 10⁻¹²no preconditionersketchedthe cost, not the answersteps, preconditioned10steps, without150answer gap, preconditioned5.7·10⁻⁶κ(A)·u at the worst κ1.1·10⁻⁶answer gap, sketch-and-solve4.5one sketch, two things to do with itand only one of them moves the answer
Fig. 6 And the object whose ageing is hardest to see, in the essay that measured what it is worth when it is fresh.

Warm starting, which is the reuse nobody argues about

There is a fifth thing carried between members and it is so obviously right that it is worth saying why it belongs on the list: the previous answer, used as the starting point for the next member.

It is free, it is never wrong, and it is the reason the other four decisions exist at all. A cold start on member two hundred would cost as much as member one; a warm start turns each member into a correction, and a correction is small. Every measurement in this field warm-starts, which is why the chord iteration in the next essay converges at all from a factorisation four members old — the point it starts from is already close.

What warm starting does not do is make the sequence a smaller problem. The conditioning of member two hundred is the conditioning of member two hundred whatever it is started from, and an ill-conditioned member met at the end of a long sequence is exactly as ill-conditioned as it would have been met alone. The essay on a condition number scaling cannot move is the standing statement of that, and it is worth carrying into this field intact: a warm start changes the distance to the answer and not the sensitivity of the answer to the data.

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⁶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. 7 The amplifier the whole collection is organised around. A warm start moves where an iteration begins and moves nothing on this figure.
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. 8 And the same point made about a rank-one correction: a cheap update to a hard problem inherits the hard problem’s conditioning, however cheap the update was.

And one asset that must not be reused at all

The counterweight, and the field needs it early because the four above make reuse sound like a uniformly good idea with a schedule attached.

A randomised method draws a random matrix and applies it. The guarantee attached to it is a probability statement over that draw — over a random object chosen independently of what it is applied to. Use it once and deflate what it found, and the next input is a function of the draw. The independence is gone, and it does not degrade gracefully: the second sketch of the deflated matrix is the zero matrix, exactly, and the measurement is 9.0·10⁻¹⁵ where the first round saw 2.6.

So the field’s rule is not keep things. It is that every object carried between members of a sequence has a lifetime, that the lifetime is a measurement rather than a policy, and that for one of them the lifetime is a single use.

The randomised SVD against the optimum it cannot beatA semi-logarithmic plot of approximation error against target rank. A shaded band shows the spread across seeds, a solid line the optimal error from the exact singular values, and a dashed line the published probabilistic bound well above both.04812162010⁻¹10⁻⁰.⁵1target rank k‖A − A_k‖₂published boundrandomisedσ_{k+1}, optimalhow far apart the three areworst seed spread1.6bound / median at k = 125.9median / optimum at k = 121.960×60, 6 seeds, oversampling p = 5band is best to worst
Fig. 9 The randomised bound in the essay that priced it, where the guarantee is a band rather than a line because the answer depends on a draw.

What a sequence’s backward error even is

There is a question underneath all of this that is worth asking out loud, because the collection’s own vocabulary does not answer it.

A backward-stable solve returns the exact answer to a nearby problem. What does a backward-stable sequence return? Not the exact answer to a nearby sequence — the members are not independent, and a perturbation to member three changes what member four is. The honest statement is narrower and it is the one this field uses: each member is solved backward-stably for the problem it was actually handed, and the problem it was handed was built from an answer that was itself only approximate.

That distinction is what makes the essay on a chain of rotations necessary. There, one object is carried across three thousand steps without ever being recomputed from the data, and per-step backward stability says nothing at all about where it ends up. Here, the residual is recomputed from the data at every outer step, and per-step backward stability is exactly what is wanted.

Two situations that look the same from a single step and are not the same situation.

What the solver reported, and what it achievedFinal relative residual and final relative error against the working precision, on a logarithmic vertical axis. The residual sits at the stopping tolerance at every precision — the method converged, by every test available to it. The error runs from 5.26·10⁻¹⁴ at 53 bits to 0.00804 at 8, tracking the unit roundoff, which is drawn beside it.614223038465410⁻¹³10⁻¹⁰10⁻⁷10⁻⁴working significand bitsrelative sizeerrorresidualuthe gap the solver cannot seeerror ÷ residual at 24 bits4.1·10⁵error ÷ residual at 16 bits3.6·10⁷error ÷ residual at 8 bits1.9·10¹⁰every convergence test passesand the answer is wrong in proportion to u
Fig. 10 A small residual is not a small error, in the essay that established it. Within a sequence the residual is recomputed at every step, which is what makes the inner solve’s error survivable and does nothing whatever about the problem’s conditioning.
How far a carried triangular factor drifts from the data, over 2976 steps of a sliding windowA window of 24 rows on 6 columns, moved one row at a time: each step folds a row in with Givens rotations and removes one with hyperbolic rotations, and the factor is never rebuilt from the rows. The rows are integers times powers of two, so AᵀA is exact in a double and the vertical axis is a distance from the answer. No single step amplifies by more than 2.72, no downdate fails, and after 2976 steps the factor is 3.87·10⁻¹⁴ from the matrix it is supposed to factor — a fitted slope of 0.554 in the step count, against a bound whose slope is 1.10²10³10⁴10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²steps taken‖RᵀR − AᵀA‖ ⁄ ‖AᵀA‖the bound, linear in the steps√k · uevery step safe, the chain notdrift after the run3.9·10⁻¹⁴the bound there3.3·10⁻¹³√k · u there6.1·10⁻¹⁵worst single amplification2.7worst leverage met0.69refactorisations1backward stable onceand three thousand times is a different claim
Fig. 11 And the other situation, from the essay on a carried factor: an object that is never recomputed from the data, drifting away from it over three thousand steps.

What the four bars actually cost

The hero figure is the classification with numbers on it. Twelve members of a 120×120 problem, costed in multiplications with a dense Cholesky at n³/3 and a solve at 4n²:

  • nothing changes — one factorisation, one solve, 8.3 per cent of the independent cost, and eleven of the twelve members are free;
  • the right-hand side changes — one factorisation, twelve triangular solves, 16.5 per cent, with the worst relative residual across all twelve at 1.5·10⁻¹³;
  • the matrix drifts — four factorisations and seventy-five chord steps under a rule that decides when to rebuild, 85.9 per cent;
  • everything changes — twelve factorisations, twelve solves, and the definition of 100 per cent.

The drifting bar is the one worth staring at. At twelve members it costs 86 per cent of the independent sequence, which is a poor advertisement; at six members it costs more than the independent sequence, because the rule has not had time to keep anything and it pays for the measurement anyway. Reuse is a claim about a run and not about a step, and a field that only ever draws one step cannot see that.

What 6 members of a sequence cost, by what changes between themEvery bar is 6 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 — 24.1 per cent of the independent cost, with every solve still at a relative residual of 1.41·10⁻¹³. When the matrix drifts, the factorisation is reusable for a while and the while has to be measured: 3 factorisations for 6 members. When everything changes, nothing carries.nothing16.7%the right-hand side24.1%the matrix, slowly100.4%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 6 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations3independent: factorisations6what changes between the membersdecides what may be carried
Fig. 12 Six members, where reuse has not yet paid for itself. A policy that is right at twenty is wrong here, and nothing about a single solve says which regime it is in.

Where the accuracy goes

The last thing to establish before the field’s measurements begin is what the outer loop actually does with the inner solve’s accuracy, and it is worth doing directly rather than by argument.

Take one iterate of the Newton sequence, 2.50·10⁻³ from a root that is known exactly. Solve its linear system to a relative residual of 10⁻⁴: the resulting point is 1.200·10⁻⁵ from the root, and it cost 426 conjugate gradient iterations. Solve the same system to 10⁻¹⁴: the resulting point is 1.200·10⁻⁵ from the root, and it cost 1,129.

Eleven decades of inner accuracy, the same four digits, and a factor of 2.6 in work. The floor is not the arithmetic and it is not the conditioning: it is the square of the distance the step started at, 6.23·10⁻⁶, because a Newton step is the exact solution of a linearised problem and the linearisation is wrong at second order.

A single Newton step solved to eleven inner tolerances, against how far the resulting point is from the rootThe iterate is 0.0025 from a root that is known exactly by construction. The same linear system is solved to relative residuals from 0.3 down to 10⁻¹⁴, at 109 and 1129 conjugate gradient iterations, and the resulting point is 2.837·10⁻⁴ and 1.2·10⁻⁵ from the root. The curve is flat below about 10⁻³: the linearisation is wrong at second order, so the step cannot land closer than the square of the distance it started at — 6.234·10⁻⁶ — however exactly it is computed.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²inner tolerance η, relative residual of the linear solvedistance from the root after the stepd² = 6.23·10⁻⁶distance before the step, 0.00251093397131129the number by each point is the iterations it costeleven decades, one landing placedistance before the step0.0025its square6.2·10⁻⁶where η = 10⁻³ lands1.2·10⁻⁵where η = 10⁻¹⁴ lands1.2·10⁻⁵iterations for the first339iterations for the second1129the accuracy that is thrown awaymeasured against a root that is known
Fig. 13 The measurement that the rest of this field is built on, at the iterate closest to the root. Eleven decades along the bottom, one landing place.

That is the whole of the field’s first claim, and everything else in it is the same claim about a different object: the tolerance, the factorisation, the preconditioner and the pattern are all approximations to something, the outer loop measures the difference from the matrix at every step, and what may be approximated is decided by what is measured again rather than by what looks important.

Backward error with and without the structure, for Levinson and for elimination, n = 12Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 3.1·10⁻¹⁷ at every ρ, and the two curves are indistinguishable. Insisting the perturbation be a symmetric Toeplitz matrix — the kind of object the problem was posed with — lifts both by orders, and at ρ = 0.999 Levinson's rises by a factor of 2.48·10⁶. Neither number is alarming in absolute terms; what is worth carrying is that the reassuring one is the one that is reported, and the two are not measuring the same thing.10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation3.1·10⁻¹⁷Levinson, Toeplitz only7.6·10⁻¹¹the ratio between them2.5·10⁶diagonal defect of the first0.97the number that is reportedand the number that is asked about
Fig. 14 The measurement the whole collection is named for, from the essay that established it. In a sequence it is computed once per outer step, which is what makes the inner object’s error forgivable.
Newton on the Bratu problem: analytic Jacobian, differenced Jacobian at ε = 10^-6, and matrix-freeThree residual sequences, all starting from zero. The analytic Jacobian gives 8, 0.052, 3·10⁻⁶, 5.2·10⁻¹³; a Jacobian differenced at ε = 10^-6, whose entries are correct to about 1.3·10⁻¹⁰, gives 8, 0.052, 3·10⁻⁶, 4.6·10⁻¹³. They stop at the same residual. The matrix-free run, where GMRES sees only a closure and no entry exists anywhere, lands 8.1·10⁻¹⁶ from the analytic answer in 438 products.0123456710⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1Newton step‖F(x)‖analytic Jdifferenced Jmatrix-freeone fixed point, three derivativesanalytic floor4.2·10⁻¹³differenced floor3.9·10⁻¹³matrix-free floor3.8·10⁻¹³products used, matrix-free438the derivative chooses the stepand the residual decides the answer
Fig. 15 A Newton residual history from the matrix-free essay, where the same sequence is run with two different derivatives and stops at the same place.
Conjugate gradients at κ = 105, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound168 steps40×40, spectrum spread evenly in logbound permits 4479
Fig. 16 And the inner solver, whose iteration count is the currency every reuse decision in this field is priced in.

What the field will not claim

Two things are worth ruling out at the start, because both are easy readings of the essays that follow and both are wrong.

Not that a rule beats a constant. A rule that reads the residual beats every constant forcing term, and a rule that reads the contraction ratio beats every fixed refactorisation period. A rule that reads the iteration count of a preconditioned solve does not beat every fixed rebuild period: it wins in the middle of the range and loses at both ends, because what it cannot see is what a rebuild cost, and at the ends that is the whole decision.

Not that an outer loop repairs everything. It repairs a loose tolerance, a stale factorisation and a stale preconditioner. It does not repair a pivot replaced by a number nine orders of magnitude larger than the row it came from, and it does not repair a random projection that the input has already been made orthogonal to. Both of those are changes to the problem rather than errors in solving it, and a residual computed from the problem cannot see them.

The rest of this field is those measurements, in that order.

Three tolerances from three fields: work not done, against backward error acceptedA deflation criterion writes a zero into a subdiagonal entry of a Hessenberg matrix; a drop tolerance discards an entry of an incomplete Cholesky factor; a truncation discards a singular value. The three are described in different vocabularies, live in different fields and share no arithmetic. Plotted as the fraction of the work not done against the relative backward error the choice introduces, they are three samples of one curve: fitted slopes of 0.039, 0.051, 0.225 of the work a decade of error. Nothing is rescaled to make that happen. The low-rank curve goes below zero at the tight end, which is the part nobody plots: a rank-k factored form of an n×n symmetric matrix costs more to store than the matrix once k passes n/2.10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²-1-0.500.51relative backward error acceptedfraction of the work not doneQR deflation criterionincomplete Cholesky droplow-rank truncationthree fields, no shared arithmeticslope, deflation criterion0.039slope, drop tolerance0.051slope, rank truncation0.23widest apart, as a ratio5.7a tolerance is an offerand the three offers are one offer
Fig. 17 The nearest thing this collection already had to a field about choices: three tolerances from three fields on one pair of axes. Every one of them is a single-solve decision, and every one of them changes when the solve is inside a loop.
Four knobs on one problem at 1.0% noiseRelative error against the fraction of each method's own range, on a logarithmic vertical axis. A truncation, a Tikhonov parameter, a conjugate gradient step count and a randomised rank each have an interior minimum, and the four minima are 0.1445, 0.1406, 0.1426, 0.1449 — a spread of 3%. The horizontal line is the best of them.00.250.50.75110⁻¹110¹10²10³fraction of the method's own rangerelative errorfloor 0.141truncation KTikhonov λCGLS steprandomised rankfour methods, one floortruncation K0.14Tikhonov λ0.14CGLS step0.14randomised rank0.14four knobs from four fieldsand one obstruction underneath them
Fig. 18 And the essay that found four knobs with one floor under them. This field is what happens to those knobs when the problem they are set for arrives two hundred times.
Where the Arnoldi recurrence stops, against the tolerance it calls a zero, n = 27The spectrum has 3 coarse clusters of width 10^-4 and 9 fine ones inside them, so this one matrix has three near-invariant subspaces and all three are real. The tolerance decides which the recurrence stops at: it stops at 3 at the loose end, at 9 in the middle, and at 27 at the tight end. The residual it leaves runs from 1.51·10⁻⁵ to 5.48·10⁻¹⁷. Neither end is a mistake and neither is a trade: stopping too early throws away an answer that was available, and stopping too late means dividing by rounding error.10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²0612182430tolerance called a zerostep the recurrence stops atsubspace of 3subspace of 9subspace of 27too late: dividing by roundingtoo early: an answer thrown awayone matrix, three answerssteps at the tightest tolerance27at the loosest3residual left, tightest5.5·10⁻¹⁷residual left, loosest1.5·10⁻⁵both directions are failuresand the constant decides which
Fig. 19 A detection threshold, which is the one kind of constant this field does not reach: both of its directions are failures, so there is nothing for an outer loop to forgive.
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. 20 And the shape the essay on a kept factorisation opens with: what a factor costs the members it is kept for, with the ringed points marking where it was rebuilt.

Where this field touches the others

It is worth naming the essays this one is built on, because the field is not a new subject — it is every existing field’s object, met a second time.

‖QᵀQ − I‖ at every step of a Lanczos run, n = 40The departure from orthogonality of the Lanczos basis, against the step, on a logarithmic vertical axis. It sits at the level of rounding for 13 steps and then climbs by a factor of about seventeen a step for ten steps running, twelve orders of magnitude, before saturating. The dashed line marks the step at which the first Ritz value converged, which is step 11.051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹step‖QᵀQ − I‖first Ritz value converges√uno reorthogonalisationfullthe two events are one eventstep orthogonality crosses √u14step the first Ritz value converges11consecutive steps at ≥5× growth10a drift would grow like a square rootthis grows like a geometric series and then stops
Fig. 21 Orthogonality measured rather than assumed, from the field that made it a number. A basis carried across restarts is the same kind of object as a factorisation carried across members.
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. 22 What a preconditioner is worth when it is built for the matrix it is used on — the baseline every measurement of a stale one is against.
The 12×12 grid Laplacian and its Cholesky factor, ordered by reverse Cuthill–McKeeTwo square sparsity plots side by side. The left shows the nonzeros of the matrix; the right shows the nonzeros of its Cholesky factor, with the entries created by elimination marked in a second colour.the matrix, lower triangle408 entriesits Cholesky factor1354 entries · 946 created→‖A − LLᵀ‖/‖A‖1.1·10⁻¹⁶fill, symbolic946fill, numeric946n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding
Fig. 23 A sparsity pattern and the order that produced it. In a sequence the pattern is the same at every member and the order that was right for the first one may not be.
Sketch distortion against sketch width, for 200 vectorsA log-log plot of the worst relative change in vector length against the number of rows in the sketch, for vectors of two dimensions a factor of four apart. The two curves lie almost on top of one another and both fall steadily.10²10².³10².⁵⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁶10⁻¹10⁻⁰.⁵1rows in the sketchworst relative distortiondimension 64dimension 2565 seeds per point, band is best to worstthe dimension does not appear
Fig. 24 And the randomised object whose whole guarantee is about a draw the input has not seen.
Separating two eigenvalues, against how close they areIterations against the gap between the two largest eigenvalues, on a logarithmic gap axis. The single-vector method needs 17, 18, 20, 22 steps as the gap closes through four decades. The block of two needs 11 at every gap — and spends 22 products with A doing it, which is no less arithmetic. What it saves is synchronisations.10⁻⁴10⁻³10⁻²10⁻¹0510152025gap between the two eigenvaluesiterationsone vectorblock iterationsblock productsiterations, not arithmeticsingle-vector steps at 0.117single-vector steps at 0.000122block iterations, every gap11the same products with Aand half the synchronisations
Fig. 25 The price of a restart, measured where it was chosen rather than forced — the closest thing the collection already had to a lifetime for a carried object.

Where this goes next

the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ The exact answer to a nearby problem A good algorithm does not give an approximate answer to your problem. It gives the exact answer to a problem very close to yours — and once that is the definition, a wrong result has two possible authors and they can be measured apart. 10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice An operator with no entries At the sizes where linear algebra is expensive the matrix does not exist. What exists is a subroutine that returns Av. Every Krylov method survives that unchanged; every algorithm that reads an entry disappears. And the derivative such a code computes is accurate to ten digits instead of sixteen, which turns out to cost nothing at all. 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.

Backward errorCholesky factorisationCondition numberExact ground truthFlop countNewton iterationPreconditioningTriangular solveWarm start