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

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.
30 min read 25 figures Where the algebra stops being the arithmeticMeasured, not assertedExact ground truthAssertions that reject

Worth reading first: The problem that arrives again · The rate the condition number predicts.

Newton’s method for F(x) = 0 does one thing per step: it replaces F by its tangent at the current point and solves the linear problem exactly. Everything anybody knows about the method — that it doubles the correct digits, that it needs a good starting point, that it converges quadratically near a simple root — is a statement about that replacement rather than about the linear solve.

At the sizes where any of this is expensive, the linear solve is iterative and it stops when a tolerance says so. So there is a number to choose, it is called the forcing term, and the natural instinct is that it should be small: the step is what the method is made of, and a step computed badly is a bad step.

The instinct is wrong, and it is wrong by a factor that can be measured rather than argued about.

A single Newton step solved to eleven inner tolerances, against how far the resulting point is from the rootThe iterate is 0.0372 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 1126 conjugate gradient iterations, and the resulting point is 0.005319 and 0.002497 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 — 0.001383 — however exactly it is computed.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
Fig. 1 One Newton step, solved to ten inner tolerances, against how far the resulting point is from a root that is known exactly. The number beside each point is the conjugate gradient iterations it cost.

The plateau

The measurement is direct and needs no theory to state. Take an iterate 3.719·10⁻² from the root of the model problem — a root that is known exactly, because the right-hand side was built from it. Solve the Newton system at that point to a relative residual of 10⁻³, and the point it steps to is 2.4966·10⁻³ from the root, at 354 conjugate gradient iterations. Solve the same system to 10⁻¹⁴, and the point it steps to is 2.4967·10⁻³ from the root, at 1,126.

Eleven decades of inner accuracy. Four parts in ten thousand of difference in the answer. A factor of 3.2 in work.

That is the whole essay and the rest of it is about why the floor is where it is, how far the effect extends, and what it costs to ignore.

A single Newton step solved to eleven inner tolerances, against how far the resulting point is from the rootThe iterate is 0.507 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 78 and 1005 conjugate gradient iterations, and the resulting point is 0.2069 and 0.1787 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 — 0.2572 — however exactly it is computed.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110⁻²10⁻¹110¹inner tolerance η, relative residual of the linear solvedistance from the root after the stepd² = 0.257distance before the step, 0.507782576211005the number by each point is the iterations it costeleven decades, one landing placedistance before the step0.51its square0.26where η = 10⁻³ lands0.18where η = 10⁻¹⁴ lands0.18iterations for the first257iterations for the second1005the accuracy that is thrown awaymeasured against a root that is known
Fig. 2 The same experiment three steps earlier, where the iterate is 5.071·10⁻¹ from the root. The plateau is at 1.7869·10⁻¹ and its left edge is in the same place.
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. 3 And one step later, where the iterate is 2.497·10⁻³ out. The plateau is at 1.1999·10⁻⁵ and the tolerance at which it begins has not moved.

Why there is a floor at all

A Newton step solves J(x)·s = −F(x), and s is the exact correction for the linearised problem. The real problem is not linear. Expanding F about x,

F(x + s) = F(x) + J(x)s + ½·F″(ξ)[s, s]

and the step is chosen to kill the first two terms exactly, so what is left is the second-order term — of size proportional to ‖s‖², which near the root is proportional to the distance from it, squared.

The step therefore cannot land closer than a constant times d², where d is where it started, however exactly the linear system is solved. Below that, the accuracy of the linear solve is being spent on a difference between two points that are both the wrong point by the same amount.

That is a derivation, and this collection’s habit is to measure the thing rather than take the derivation’s word for it. Five iterates, five plateaus:

distance in its square where the plateau sits
1.027 1.055 5.0714·10⁻¹
5.071·10⁻¹ 2.572·10⁻¹ 1.7869·10⁻¹
1.787·10⁻¹ 3.193·10⁻² 3.7187·10⁻²
3.719·10⁻² 1.383·10⁻³ 2.4967·10⁻³
2.497·10⁻³ 6.234·10⁻⁶ 1.1999·10⁻⁵

Fitted on log axes, the plateau’s height against the distance it started at has a slope of 1.82, against a predicted 2. Not 1, which is what a floor set by the linear solve’s own error would give, and not 0, which is what a floor set by the arithmetic would give. The constant is between 0.5 and 1.9 across the five, which is order one — so the floor is the linearisation and it is not an artefact of the conditioning, which runs to 10⁵ on this problem and does not appear.

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. 4 A Newton residual history from the essay on operators with no entries, where two very different derivatives produce the same sequence — the same phenomenon seen through the Jacobian rather than through the solve.

How much it costs to ignore

One step is not a solve. The forcing term is chosen once and applies to every step, so the interesting quantity is the total inner work for a complete run to a stated outer tolerance.

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. 5 Total conjugate gradient iterations for a whole Newton solve to a relative residual of 10⁻¹⁰, against the constant forcing term chosen for it. Both arms are real and the bottom is nowhere near machine precision.

The curve is a U and both of its arms are worth naming.

The left arm is oversolving. At η = 10⁻¹⁴ the run takes 9 Newton steps and 9,358 inner iterations. At η = 10⁻⁶ it takes the same 9 Newton steps and 4,382. The outer iteration is doing exactly the same thing — the same number of steps, landing in the same places to within the plateau — and one of the two runs is paying twice as much for it.

The right arm is starving the outer loop. At η = 0.5 each step is nearly free, but the steps stop being Newton steps: the quadratic convergence goes, the run takes 30 outer steps instead of 9, and the total climbs back to 1,178.

The bottom is at η ≈ 0.1, at 980 iterations, which is a factor of 9.5 below the fully solved run. Nine and a half times the work, for a final answer whose residual is below the same tolerance.

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^-12. At η = 10⁻¹⁴ it takes 9358 conjugate gradient iterations across 9 Newton steps; at η = 0.1 it takes 1256 across 12. 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 1206 iterations and reaches 4.8·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 iterations1256adaptive rule, iterations1206adaptive final residual4.8·10⁻¹³outer tolerance asked for10⁻¹²a tolerance is a cost decisionand its optimum is not machine precision
Fig. 6 The same sweep at a tighter outer tolerance. The whole curve lifts and the minimum moves, which is the argument against a constant tuned once.
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^-6. At η = 10⁻¹⁴ it takes 8230 conjugate gradient iterations across 8 Newton steps; at η = 0.1 it takes 459 across 6. 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 593 iterations and reaches 4.47·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 iterations8230cheapest constant0.1its inner iterations459adaptive rule, iterations593adaptive final residual4.5·10⁻⁷outer tolerance asked for10⁻⁶a tolerance is a cost decisionand its optimum is not machine precision
Fig. 7 And at a looser one, where the left arm is even more wasteful because there are fewer steps to spread the waste over.

The one thing the tight run does buy

An honest accounting has to say what the extra work is not wasted on, and there is one thing.

The run at η = 10⁻¹⁴ ends with a forward error of 2.79·10⁻¹⁰. The run at η = 0.1 ends with 2.43·10⁻⁸ — eighty-five times worse. Both satisfy the stopping test, which asks for a relative residual below 10⁻¹⁰, and neither is more correct than the other by that test.

The reason for the difference is not that the loose run was less careful. It is that a quadratically converging iteration overshoots its stopping test: the tight run’s last step lands at a residual of about 10⁻¹⁶, five orders of magnitude below what was asked for, because the step before it was at 10⁻⁸ and the next one squares that. All of that accuracy is free in the sense that it was not asked for, and it is expensive in the sense that the last inner solve is the most expensive one in the run.

So the tight constant is not buying nothing. It is buying digits nobody requested, at the highest price the run has to offer, on every step rather than on the last one.

What each Newton step asked of its linear solve, and what the solve cost, under the tight 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 10⁻¹⁴ and costs 1775 iterations, and the last asks for 10⁻¹⁴ and costs 1128. The whole solve costs 9358 inner iterations across 9 Newton steps and ends at a relative residual of 1.43·10⁻¹³.01234567891010⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110³Newton steptolerance asked for, and iterations paiditerations paidtolerance asked forouter residualthe tight policy, step by stepNewton steps9inner iterations, total9358first step's cost1775last step's cost1128final outer residual1.4·10⁻¹³the rule reads the last two residualsand asks for nothing it cannot use
Fig. 8 The tight policy step by step: the same demand every time, and a cost that climbs as the Jacobian gets harder rather than as the answer gets better.
What each Newton step asked of its linear solve, and what the solve cost, under the loose 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.1 and costs 4 iterations, and the last asks for 0.1 and costs 133. The whole solve costs 980 inner iterations across 10 Newton steps and ends at a relative residual of 4.09·10⁻¹¹.0123456789101110⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110³Newton steptolerance asked for, and iterations paiditerations paidtolerance asked forouter residualthe loose policy, step by stepNewton steps10inner iterations, total980first step's cost4last step's cost133final outer residual4.1·10⁻¹¹the rule reads the last two residualsand asks for nothing it cannot use
Fig. 9 And the loose one, which is cheap per step and pays for it in the number of steps — the right arm of the U seen from inside.

That observation is what the essay on a rule that reads its own residual is about, and the shape of the answer is already visible here: what is wanted is a loose tolerance early and a tight one at the end, which no constant can be.

Where the plateau’s left edge is

The tolerance at which the inner solve stops mattering is not machine precision and it is not a fixed number either. It sits at about the square root of the plateau — which is to say at about the distance to the root — because that is where the linear solve’s own error stops being smaller than the linearisation’s.

That has a consequence worth stating plainly: the correct inner tolerance is proportional to how far from the root the iteration currently is, and the iteration knows that quantity to within a constant, because it is roughly the ratio of consecutive residuals.

At the first step of the run drawn here, the iterate is 2.95 away from the root and the plateau is at 1.80. Any tolerance below about 0.1 is spent on nothing. At the last step the iterate is 2.5·10⁻³ away and the plateau is at 1.2·10⁻⁵; there a tolerance of 10⁻² lands at 1.4669·10⁻⁵ against the plateau’s 1.1999·10⁻⁵, a genuine 22 per cent shortfall, and one more decade closes it entirely.

The demand a step can use rises by two decades for every decade the outer iteration gains, and a constant is by construction wrong at one end of that or the other.

A single Newton step solved to eleven inner tolerances, against how far the resulting point is from the rootThe iterate is 1.8 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 33 and 700 conjugate gradient iterations, and the resulting point is 1.138 and 1.027 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 — 3.238 — however exactly it is computed.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110⁻¹110¹inner tolerance η, relative residual of the linear solvedistance from the root after the stepd² = 3.24distance before the step, 1.833167416700the number by each point is the iterations it costeleven decades, one landing placedistance before the step1.8its square3.2where η = 10⁻³ lands1where η = 10⁻¹⁴ lands1iterations for the first167iterations for the second700the accuracy that is thrown awaymeasured against a root that is known
Fig. 10 The far-field case, where the whole picture is nearly flat. A Newton step taken this far from the root is barely worth computing carefully at all.
A single Newton step solved to eleven inner tolerances, against how far the resulting point is from the rootThe iterate is 0.179 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 91 and 1099 conjugate gradient iterations, and the resulting point is 0.04787 and 0.03719 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 — 0.03193 — however exactly it is computed.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110⁻²10⁻¹1inner tolerance η, relative residual of the linear solvedistance from the root after the stepd² = 0.0319distance before the step, 0.179913086891099the number by each point is the iterations it costeleven decades, one landing placedistance before the step0.18its square0.032where η = 10⁻³ lands0.037where η = 10⁻¹⁴ lands0.037iterations for the first308iterations for the second1099the accuracy that is thrown awaymeasured against a root that is known
Fig. 11 And the middle of the run, where the plateau has dropped by four decades and its left edge has moved by two.

What the sweep costs to make

It is worth saying what the U above is, because a curve of twelve points where each point is a whole solve is a different kind of object from a curve where each point is an evaluation.

Every point on it is a complete inexact Newton run on the same 200×200 problem to the same outer tolerance: nine to thirty outer steps, each with its own conjugate gradient solve, and the vertical coordinate is the sum of every inner iteration in the run. Nothing is estimated and no step count is extrapolated. The adaptive point drawn beside the curve is a thirteenth complete run under a different rule, and it is a point rather than a curve because it is not a constant and there is no axis it belongs on.

The reason to say so is that the two arms of a U are usually a fit through a model. These are measurements of runs, and the asymmetry between the arms is a measurement too: the left arm rises smoothly, because oversolving costs a predictable number of extra iterations per decade, and the right arm rises in jumps, because a Newton step either survives being solved that loosely or the run needs another whole step.

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. 12 The thirteenth run, opened up: what each step asked for and what it paid. The staircase is the shape no constant has.
Adding 1 to 256, a thousand times, at 24 significand bitsThree lines against the number of additions: the exact total rising steadily, a staircase following it, and a flat line at the starting value.02505007501000250500750100012501500additionsrunning totalexactstochasticnearesta thousand additionshalf an ulp at 2561.5·10⁻⁵moves, round to nearest1000moves, stochastic1000relative error, nearest0relative error, stochastic024 significand bits, unbounded exponenta flat line is not a small error
Fig. 13 An iteration that stops improving, from the arithmetic field. The right arm of the U is the same failure arriving from the opposite direction — not too little precision, but too little demand.
Forward error of a 6×6 Hilbert solve at eight precisionsA bar for each significand width from 12 to 53 bits showing the relative error in the computed solution, with the condition number times the unit roundoff marked as a prediction.κ = 1.5·10⁷ · the exact answer is (1, 2, …, 6)12 bits3.316 bits0.8620 bits4.324 bits0.09330 bits7.6·10⁻⁴36 bits6.2·10⁻⁶43 bits2·10⁻⁷53 bits4.5·10⁻¹¹dashed: κ · unit roundoffone matrix, eight arithmeticsmeasured against a known answer
Fig. 14 And the knob this collection is named for. The forcing term is a second knob on the same computation, and unlike the mantissa it can be moved between one step and the next.

What this is not

Three readings to rule out, because each of them is available from the figure and each is wrong.

Not that the inner solve does not matter. It matters enormously; it is where all the work is. What does not matter is the part of its accuracy that lies below the linearisation error, and that part is most of the accuracy a default tolerance asks for.

Not that this is about conjugate gradients. The plateau is a property of the outer method. Any inner solver put in its place meets the same floor; what changes is how much a decade of tolerance costs, and that only moves the vertical distance between the points on the plateau, not the plateau.

Not that the effect is small on well-conditioned problems. It is invisible on well-conditioned problems, which is different: there the inner solve costs a handful of iterations at any tolerance, so the waste is real and cheap. The effect scales with how expensive the inner solve is, which is to say it is largest exactly where it is worth having.

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. 15 Why the inner solve is expensive here: the conjugate gradient rate the condition number predicts, at the conditioning this problem is built to.
The residual basis of conjugate gradients, at κ = 105A semi-logarithmic plot against iteration count showing the loss of orthogonality among the residual vectors rising while the relative residual falls.05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here
Fig. 16 And why it is more expensive than that: the loss of orthogonality that stops conjugate gradients from being a direct method, which is the reason 1,126 iterations are possible on a 200×200 system.

Two numbers a solver could report and does not

The measurement suggests two things an inexact Newton implementation could print, neither of which is expensive and neither of which is usual.

The linearisation error of the step just taken, which is ‖F(x + s)‖ compared against what the linear model predicted. It costs nothing extra — the residual at the new point is evaluated anyway on the next step — and it is exactly the floor this essay is about. A solver that printed it would be telling its user how much inner accuracy the next step can use.

And the share of the inner solve that was spent below that floor, which is the difference between the tolerance asked for and the tolerance the floor justified. On the runs here that share is most of the work, and it is invisible in any log a solver currently writes.

Neither is a new algorithm and neither changes an answer. They are both the collection’s standing habit applied to somebody else’s code: a quantity that decides something should be printed beside the thing it decided, so that a reader can see whether the decision was a good one.

The refusal

The claim that a Newton step landing short of the root was solved too loosely is a claim about which of two quantities is the binding one, and it is fed here the case where the answer is the other one.

At the iterate 3.719·10⁻² from the root, the seven tolerances from 10⁻³ down to 10⁻¹⁴ land between 2.4966·10⁻³ and 2.4967·10⁻³. The assertion that the shortfall belongs to the inner solve is fed those numbers and required to fail, which it does. Then the same assertion is fed the first Newton step from a cold start — where the tolerance genuinely does move the answer, because the linearisation error there is the size of the iterate — and it passes, which is what makes the first result a measurement rather than a tautology.

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. 17 The distinction underneath all of this, from the essay that established it: the two errors, and which of them a stopping test can see.
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. 18 And the standing warning that a small residual is not a small error — here in its sequence form, where the residual being small is what licenses the next step rather than what ends the computation.

What it means for the rest of the field

The plateau is the first of the four assets this field is about, and it is the one with no downside at all. A tolerance is not an object that has to be built, so choosing it badly costs nothing but work; there is no risk of a wrong answer, no cliff, and no case where being too loose destroys anything that cannot be recovered on the next step.

That is exactly why it is first. The other three — the factorisation, the preconditioner, the pivot order — are objects with a construction cost, and each of them adds a failure mode that a tolerance does not have.

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. 19 The next object along: a factorisation, which is cheap to keep and has a member past which keeping it stops working.
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. 20 And the shape that decision takes, where the failure is a cliff rather than a slope.
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 The preconditioner, whose ageing depends on where the drift landed rather than on how large it was.
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. 22 And the pivot order, which is the one this field’s repair does not reach.
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. 23 All four in one picture, as the cost of a sequence by what changes between its members.
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. 24 And the collection’s earlier finding about tolerances, which this one sits beside: three thresholds from three fields tracing a single curve of work saved against error accepted.
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. 25 Four knobs with one floor under them, from the field about methods designed apart. The forcing term is a fifth knob and its floor is not the problem’s — it is the linearisation’s.

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⁻⁹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. 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 A small residual is not a small error Substituting the answer back and finding that it fits is the most natural check there is, and it verifies the wrong thing. A residual of 10⁻¹⁷ is entirely compatible with an answer whose second digit is wrong.

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.

Conjugate gradientsExact ground truthFlop countForcing termInexact newtonKrylov subspaceLinearisationNewton iterationQuadratic convergence