Numerical Linear Algebra
where the algebra stops being the arithmetic
Two errors, and whose fault they are

The zero you are allowed to write

A deflation criterion sets a subdiagonal entry to zero because it is small. A drop tolerance discards an entry of a factor because it is small. A truncation discards a singular value because it is small. Three fields, three vocabularies, no shared arithmetic — and plotted as work saved against error accepted, one curve.
28 min read 26 figures Whose fault is itMeasured, not assertedIdentical algebra, different arithmeticTwo routes to a number

Worth reading first: The exact answer to a nearby problem · The algorithm the libraries actually run · Changing the condition number on purpose.

Every essay in this collection’s current run has been about a quantity going to zero on its own: a divisor vanishing, a curvature changing sign, an eigenvalue leaving for infinity, a determinant that was identically zero before anybody looked.

This one is about the other half, and it is the half a program actually contains. A zero is never met in floating point. It is declared. Every one of those events reaches the code as a comparison against a constant, and somebody chose the constant.

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. 1 Three tolerances, from three fields, on one pair of axes. Nothing has been rescaled.

Three of them, in three vocabularies

The deflation criterion. The QR algorithm splits its problem when a subdiagonal entry is small compared with its neighbours, and continues on the smaller block. Deflation is what makes the whole algorithm affordable — without it every iteration would cost the full n² however much had already converged — and it works by writing a zero into a matrix entry that is not zero.

Driving the subdiagonal to zero, with λ₄/λ₃ = 0.60A semi-logarithmic plot of the magnitude of the subdiagonal entry against iteration count for three shift strategies. The unshifted curve is a straight line; the two shifted curves plunge to the bottom of the plot within a few steps.0918273645546310⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration|subdiagonal entry|no shiftRayleighWilkinsontwo routes to one raterate, from the spectrum0.6rate, measured0.6iterations, none / Wilkinson8.6symmetric 4×4, spectrum 8, 4, 2, 1.2the dashed line is the prediction
Fig. 2 The iteration whose cost that decision governs, from the field’s own essay.
After reflector 3 of 5: the bulge in column 2A square matrix drawn as a grid with one nonzero subdiagonal. A small triangle of entries below that subdiagonal is highlighted, and moves down the matrix from panel to panel.0.46-1.1-1.60.19-0.047-0.181.90.69-0.440.026-0.041-0.01601.6-0.120.078-0.049-0.11000.192.60.920.04800-0.0610.912.11.500-0.230.0491.52.3after reflector 3the highlighted entries are the bulge — the only thing that is not Hessenbergentries below the subdiagonal3reflectors used3reflectors in a whole step5the shifts are never formed — only their sum and productand both of those are real
Fig. 3 And the step it performs between deflations, which is where the work is.

The drop tolerance. An incomplete factorisation keeps a factor sparse by discarding entries that are small. The exact factor of a sparse matrix is not sparse — this collection has a whole field about that — so an incomplete one keeps a chosen subset, and the entries not kept are set to zero.

Fill growth under natural: the factor rises as n^1.49A log-log plot of nonzero count against matrix dimension. The matrix's own count is a straight line of slope one; the factor's is steeper; a dense factor is steeper still.10²10².³10²10³10⁴n (dimension)nonzerosdense factorfactor, n^1.49matrix, n^1.05grid Laplacians from 5×5 to 12×12fitted, not quoted
Fig. 4 Why an incomplete factorisation exists: the fill an exact one produces, from the sparsity field.
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 12, 12, 12, 12, 12 — 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 κ, 6.5·10⁻⁶ at the worst. The SAME sketch used as an answer instead of as a preconditioner is 10 away.10²10⁴10⁶10⁸10¹⁰110¹10²10³condition number of the matrixLSQR steps to 10⁻¹²no preconditionersketchedthe cost, not the answersteps, preconditioned12steps, without150answer gap, preconditioned6.5·10⁻⁶κ(A)·u at the worst κ1.1·10⁻⁶answer gap, sketch-and-solve10one sketch, two things to do with itand only one of them moves the answer
Fig. 5 And what the resulting preconditioner is worth, in the iterations it saves.

The truncation. A low-rank approximation keeps the k largest singular values and discards the rest, which is setting n − k singular values to zero.

Error of the best rank-k approximation to a 10×10 matrixApproximation error against k on a logarithmic axis, with the measured error and the next singular value drawn as separate curves lying exactly on top of one another.12345678910⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 4.3·10⁻⁹worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁4.3·10⁻⁹worst Frobenius discrepancy4.3·10⁻⁹κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁
Fig. 6 And the theorem that prices it exactly: the error of the best rank-k approximation is the next singular value.

What they have in common

Each of them is a perturbation of the problem, chosen deliberately.

The deflated matrix is not the matrix that was handed over; it differs from it by whatever was written into that subdiagonal entry. The incomplete factor is the exact factor of nothing; ‖A − LLᵀ‖ is what was discarded and it is not rounding. The truncated matrix is not the matrix; it differs by σ_{k+1}.

So all three are backward errors that somebody chose, in a subject where every other backward error is something the arithmetic did. And that reframing is what makes them comparable, because a backward error is a backward error whatever produced it.

The other axis is what the choice bought: iterations not taken, entries not stored, storage not used.

Which of them is the odd one out

Before the comparison, one distinction that the comparison would otherwise flatten.

Two of the three tolerances are chosen by a library author and one is chosen by the caller. The deflation criterion is buried in an eigensolver — a caller can rarely reach it and would not know what to set it to. The drop tolerance is an argument to every incomplete factorisation routine there is, with a default the caller usually keeps. The truncation is entirely the caller’s: choosing a rank is what a low-rank approximation is.

That is a distinction about visibility rather than about mathematics, and the mathematics is what this essay compares. It matters afterwards, though, because it says where the finding is actionable. A caller who understands that a rank is a backward error can price it; a caller who understands that a deflation criterion is one has learned something about a number they cannot change.

What they can do with it is read a result differently. An eigenvalue returned by a library is the exact eigenvalue of a matrix perturbed by the deflation criterion as well as by the rounding, and the first of those is usually the larger.

Householder reduction to Hessenberg form, on a symmetric 7×7Two matrices side by side. The first is full; the second has zeros everywhere below the first subdiagonal, and being symmetric is tridiagonal.2.7-1.10.0332.4-4.21.12.6-1.15.21.70.27-0.66-0.36-2.50.0331.7-0.16-1.9-1.4-1.1-0.682.40.27-1.93.31-1.2-2.6-4.2-0.66-1.41-4.1-1.8-0.0911.1-0.36-1.1-1.2-1.81.90.652.6-2.5-0.68-2.6-0.0910.653.6A, symmetric→2.75.7000005.7-1.63.9000003.93.34.6000004.60.72-2.500000-2.51.61.9000001.93.71.5000001.52H = QᵀAQ, tridiagonal‖A − QHQᵀ‖/‖A‖10⁻¹⁵below the subdiagonal0worst eigenvalue movement1.4·10⁻¹⁴a similarity, so the spectrum is untouched — and every later step is O(n²) rather than O(n³)one reduction, then every iteration is cheapthe eigenvalues did not move
Fig. 7 The reduction that precedes it, which changes no eigenvalue at all — the contrast that makes the deflation’s perturbation visible as a choice.
How far a perturbation of size ε moves an eigenvalue, 10×10A log–log plot of eigenvalue movement against perturbation size. The symmetric case lies on a line of slope one; the non-symmetric case lies on a line of slope one eighth, and at a perturbation of ten to the minus sixteen it has already moved by a hundredth.10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/10)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision
Fig. 8 And what a perturbation of a matrix does to its eigenvalues, which is the amplification the criterion’s error is subject to.

One curve

Sweeping each tolerance over its useful range and plotting the fraction of the work not done against the relative backward error accepted:

  • deflation criterion: from 6.7·10⁻¹³ of error and no saving, to 6.6·10⁻² of error and 55 per cent of the iterations not taken;
  • drop tolerance: from the unit roundoff and no saving, to 0.43 of error and 90 per cent of the factor’s entries not stored;
  • truncation: from 6.2·10⁻⁹ and a negative saving, to 7.8·10⁻² and 75 per cent of the storage not used.

Fitted slopes, in fraction of the work saved per decade of error accepted: 0.039, 0.051 and 0.225. Three fields, no shared arithmetic, no shared vocabulary, and the widest pair of those is a factor of 5.7 apart.

Nothing is rescaled to make that happen. The two quantities are defined the same way for all three because they are the same two quantities, and the agreement is what comes out.

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.188, 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.19slope, rank truncation0.23widest apart, as a ratio5.7a tolerance is an offerand the three offers are one offer
Fig. 9 On a smaller grid the drop curve shifts and its slope does not: the relationship survives a change of problem.
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.057, 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.057slope, rank truncation0.23widest apart, as a ratio5.7a tolerance is an offerand the three offers are one offer
Fig. 10 And on a larger one, where there is more fill to discard and the same tolerance saves a larger share of it.

The negative end nobody plots

The low-rank curve goes below zero at the tight end, and that is worth a sentence because it is the part of the trade that gets left off.

A symmetric n×n matrix is n(n+1)/2 numbers. A symmetric rank-k factored form is k(n+1). So the factored form costs more than the matrix once k passes n/2, and at a tolerance tight enough to keep most of the spectrum the “low-rank approximation” is a more expensive object than the thing it approximates.

Every plot of approximation error against rank stops before that point. It is not wrong to stop there — nobody would use a rank-37 approximation of a 40×40 matrix — but the omission is what makes the trade look like a free lunch at the tight end rather than a bad deal.

What a rank-10 approximation can achieve, by spectrumA semi-logarithmic plot of singular value against index for three spectra — geometric decay, algebraic decay, and flat — with the rank-ten approximation error marked on each.1112131415110⁻³10⁻²10⁻¹1index jσⱼσ11geometric, 0.85ʲalgebraic, j^−0.90flatmeasured, and equal to σ₁₁rank-10 error, geometric0.2rank-10 error, algebraic0.12rank-10 error, flat160×60, spectrum chosen rather than the entriesthe matrix decides, not the method
Fig. 11 The spectra where truncation is worth anything, from the randomised field: the decay decides how far left the useful part of the curve extends.
The randomised SVD against the optimum it cannot beat, with 2 power iterationsA 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.1bound / median at k = 1211median / optimum at k = 12160×60, 6 seeds, oversampling p = 5band is best to worst
Fig. 12 And what a randomised method’s error looks like on the same axis, which is a band rather than a line.

Why the slopes are not identical, and why that is the right result

A factor of six between the slopes is not agreement to three digits, and reading it as one would be overclaiming.

What the three share is the sign and the order of magnitude: in each field, one decade of accepted error buys somewhere between four and twenty-two per cent of the work. What differs is how much structure each problem has for the tolerance to exploit.

The truncation slope is the steepest because a decaying spectrum is the most compressible thing on the list — one decade of tolerance discards many singular values when the decay is geometric. The deflation slope is the shallowest because the subdiagonal entries the criterion is testing are falling quadratically once convergence starts, so a decade of tolerance moves the deflation point by only a step or two.

Both of those are statements about the problem. Which is the useful form of the finding: the trade exists in every field, its slope is set by how much structure the particular problem has, and neither of those facts is visible from inside any one of the three vocabularies.

Growth factor under partial pivoting to n = 24: 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.0510152025110¹10²10³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 = 24: bound 8.4·10⁶, worst 3.6
Fig. 13 A quantity that is bounded by theory and never attained in practice, which is the shape a tolerance’s worst case usually has.

A fourth, which does not fit

It is worth naming one tolerance that belongs to the same family and is not on the curve, because the reason it is not is instructive.

The stopping tolerance of an iterative solver looks like the others: a constant, compared against a residual, deciding when to stop. It is not a backward error somebody chose, though. Stopping early does not perturb the problem — it returns an approximate answer to the exact problem, which is a forward error, and the whole point of the site’s error identity is that those are different positions.

The distinction has a consequence that is easy to state and easy to get wrong. Loosening a drop tolerance changes the preconditioner and not the answer; loosening a stopping tolerance changes the answer and not the problem; loosening a deflation criterion changes the problem. Three constants that a caller meets in the same routine, in the same units, doing three different things to the error identity.

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: κ · u40×40, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 14 The gap the stopping tolerance actually controls, and the reason it is not the one it appears to control.
Conjugate gradients on an ill-posed problem at 1.0% noiseTwo curves against the step count on a logarithmic vertical axis. The relative residual falls at every one of the 120 steps without exception. The error against the true signal falls to 0.1426 at step 20 and then climbs, reaching 6.02 by the end — 42.2 times its best value.015304560759010512010⁻²10⁻¹110¹steprelative sizeleast error: 20discrepancy stop: 7errorresidualthe knob is an integerleast error, at step20error there0.14error at step 1206the residual falls at every stepthe error turns and keeps rising
Fig. 15 And the case where stopping early is not an approximation but the answer: an iteration whose error turns round while its residual keeps falling.

The knee, which is the only thing anybody actually wants

A trade curve is used by looking for its knee — the place where more error stops buying much work, or more work stops buying much error — and the three curves here have one in the same place for the same reason.

Below about 10⁻¹⁰ of accepted error, all three are flat: tightening further saves nothing because there is nothing left to keep. Above about 10⁻³ all three are steep and the answer is being visibly damaged. The useful band is the four decades in between, and the default in every library this collection has looked at sits inside it.

That is a mildly reassuring result and it is worth stating as such, because the essays around it are not reassuring. Defaults chosen decades ago by three different communities, for three different reasons, land in the same part of a curve none of them plotted.

Fill against growth as the pivot threshold moves, on the 8×8 gridTwo curves against the pivot threshold on a logarithmic horizontal axis. One falls steeply from left to right; the other rises gently. A vertical line marks the value libraries default to.10⁻³10⁻²10⁻¹1110¹10²10³10⁴pivot threshold τgrowth factor · entries in L+U, ÷ entries in Agrowthfillthe library default‖PA − LU‖/‖A‖ at τ = 0.12.1·10⁻¹⁶growth at τ = 0.138entries in L+U at τ = 0.1875one knob, two measurements, opposite directionsand the default is most of both
Fig. 16 The same knee in the sparse-pivoting field, where the tolerance trades fill against growth.
Convergence factor against the strength threshold, ε = 0.01A curve of convergence factor against the strength threshold on a logarithmic axis. It is flat at a high value to the left of a marked position and flat at a low value to the right of it.10⁻³10⁻²10⁻¹00.20.40.60.8strength threshold θresidual reduction per cycleθ = εsemi-coarseningkept whole rowsfull coarseningone parameter, two methodsbest factor above ε0.047best factor below ε0.46the ratio across the switch9.731×31 anisotropic operatora switch, not a dial
Fig. 17 And in the algebraic multigrid essays, where a threshold decides which of two behaviours a method has.

Where the defaults came from

Since the three defaults land in the same band, it is worth asking what put them there, and the answers are unrelated.

The deflation criterion in every implementation descended from the original QR algorithm is a small multiple of the unit roundoff times the sum of two neighbouring diagonal entries. That form has a reason: it makes the perturbation relatively small compared with the local scale, so a deflation does not damage a small eigenvalue in a matrix with a large one. It was chosen to be as tight as it could be while still firing, which puts it at the flat left end of its curve.

The drop tolerance’s usual default of 10⁻³ or 10⁻⁴ came from a different argument entirely: it is what was found, empirically, to produce a factor that fits in memory while still preconditioning usefully. That is a statement about the other axis — the work — and it was tuned by watching iteration counts rather than accuracy.

The truncation rank is usually chosen by a caller from a plot of singular values, by eye, at a visible gap or a visible knee. That is neither axis; it is a judgement about the data.

Three constants arrived at by three arguments — tightness, memory, and a picture — landing within four decades of each other on a curve none of the three communities drew. That is the mildly cheerful result at the centre of this essay, and the fact that it is a coincidence is the reason it is worth recording rather than assuming.

Nonzeros in the Cholesky factor of the 10×10 grid Laplacian, by orderingA horizontal bar chart comparing the number of nonzero entries in the Cholesky factor under four elimination orderings, with reference marks for the matrix itself and for a dense factor.natural1009reverse Cuthill–McKee805minimum degree656nested dissection867matrix: 280 entries · dense factor: 5050bandwidth 10 · 3.60× the matrixbandwidth 10 · 2.88× the matrixbandwidth 82 · 2.34× the matrixbandwidth 75 · 3.10× the matrixn = 100, five-point stencilevery ordering fills in; none avoids it
Fig. 18 The memory argument in its own field, where an ordering rather than a tolerance is what decides whether a factor fits.
Singular values of a rank-4 matrix with noise of relative size 10⁻⁴Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 820decisivetrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 19 And the picture the third of them is chosen from, which this collection has already argued is a judgement rather than a measurement.

What this is not

Two things it would be easy to take from the above and should not be.

It is not that the three tolerances are interchangeable. They perturb different objects and the perturbations propagate differently. A deflation criterion’s error enters as a perturbation of the matrix whose eigenvalues are wanted, so it is amplified by the eigenvalue conditioning. A drop tolerance’s error enters as a perturbation of a preconditioner, which affects the rate and not the answer — the iteration converges to the solution of the original system whatever the preconditioner is. A truncation’s error enters as a perturbation of the answer itself, with no amplification at all.

Three positions in the error identity, and only the first of them is amplified by a condition number.

Backward error, forward error and the condition number, with measured valuesTwo boxes at the top — the problem posed and the nearby problem the algorithm answered exactly — and two answers below them, with the distances between all four labelled by numbers from a Hilbert solve.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
Fig. 20 The identity those three positions sit in, from the error field.
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 2n-row sketch, is 6.1286 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)6.1κ(SU), the other route6.1κ(A), across the sweep10⁸κ(AR⁻¹), across the sweep1the sketch never sees the spectrumand the spectrum cancels out of the answer
Fig. 21 And the reason the middle one is the mild case: a preconditioner changes the rate, not the fixed point.

And it is not an argument for tuning them. The flat left end of every curve says that tightening a tolerance below the knee costs work and buys nothing measurable, and the steep right end says loosening it past the knee buys work and costs the answer. The defaults are in the band. The value of the picture is understanding what kind of quantity is being chosen, not choosing it differently.

What the rest of this site does with the same idea

There is a fourth instance on this site that is not a tolerance and is the same trade, and it is worth finishing on because it makes the shape general rather than particular to a comparison.

Every mixed-precision result in this collection is a point on this curve. Computing a preconditioner at three significand bits rather than fifty-three is accepting a backward error in the preconditioner — a large one — in exchange for work not done. The essays that measured it found the preconditioner roundable to three bits and the residual not roundable at all, which is the same statement as the distinction above: the preconditioner’s error is in the mild position of the error identity and the residual’s is not.

So the general form of what this essay has been about:

Any decision to compute something less exactly than possible is a backward error accepted in exchange for work not done, and it can be plotted. Precision, sparsity, rank, deflation, and the whole of iterative refinement are the same trade in different units, and a code contains dozens of them that were each chosen by somebody looking at one of them alone.

Preconditioned conjugate gradients with one part of it roundedIteration count against significand bits, on the 100-unknown model problem. Rounding the preconditioner's output takes the count from 18 to 35 and leaves the answer correct to 8.8·10⁻¹³ throughout. Rounding the working arithmetic instead leaves the count at 18–300 and takes the answer to 0.159. The horizontal line is the 39 steps the unpreconditioned method takes.21018263442500102030405060708090100110120130140150160170180190200210220230240250260270280290300310320330340significand bitsiterationspreconditioner roundedarithmetic roundedno preconditionersame bits, different casualtyerror, 3-bit preconditioner8.8·10⁻¹³error, 3-bit arithmetic0.16‖A − LLᵀ‖/‖A‖ of the factor0.083a direction may be roundeda measurement may not
Fig. 22 The measurement that made that concrete on this site: which parts of a solver may be rounded and which may not, and how far apart the two answers are.
Iterative refinement from a 24-bit factorisation, κ = 10⁴A semi-logarithmic plot of forward error against refinement step. One curve falls steeply to the level of a double-precision solve; the other is nearly flat.012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in24-bitresidual indoubleone argument apartκ·u of the factorisation6·10⁻⁴double residual, final3.2·10⁻¹³same-precision, final1.3·10⁻⁴30×30, κ = 10⁴, same factors in both runsidentical cost
Fig. 23 And the machinery that buys the accuracy back afterwards, which is the same curve travelled in the other direction.

The refusal

The assertion is fed a drop tolerance of one.

At that value every off-diagonal entry of the incomplete factor is discarded, so LLᵀ is diagonal and the “incomplete factorisation” is a scaling. On the model problem the diagonal is constant, so the preconditioner is a multiple of the identity and buys exactly nothing: the preconditioned iteration takes the same number of steps as the unpreconditioned one.

A trade curve that reported it as an ordinary point would be reporting a preconditioner that does not precondition, at a plausible-looking place on the axis, with a backward error of 0.43 that looks like the honest right-hand end of a sweep. The refusal is what keeps the last point on a curve from being a different kind of object from the rest of it.

CG iterations on the model problem, with the two preconditioners a matrix-free code can and cannot haveAt 256 unknowns, conjugate gradients takes 55 iterations unpreconditioned. Diagonal preconditioning — n probes, or the analytic diagonal, and therefore available without any entries — takes 55, which is the same number, because this operator's diagonal is constant and scaling by it is scaling by a scalar. An incomplete Cholesky takes 21, a factor of 2.62, and it is defined by the sparsity pattern: there is no way to ask a subroutine for it. That factor is the price of an operator with no entries.36861361862360102030405060unknownsCG iterations to 10⁻¹⁰no preconditionerdiagonal (available)incomplete Choleskyat 256 unknownsunpreconditioned55diagonal55incomplete Cholesky21what the pattern buys2.6a Krylov method needs only productsand a preconditioner needs the entries
Fig. 24 The measurement that says why the diagonal buys nothing here, from the matrix-free essay: this operator’s diagonal is constant, so scaling by it is scaling by a scalar.
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.040, 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.04slope, rank truncation0.23widest apart, as a ratio5.7a tolerance is an offerand the three offers are one offer
Fig. 25 At eight points across the grid, between the two extremes above.
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.054, 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.054slope, rank truncation0.23widest apart, as a ratio5.7a tolerance is an offerand the three offers are one offer
Fig. 26 And at twelve, where the drop curve has shifted again and kept its slope.

What is next

Every tolerance in this essay is an offer: accept this much error, save this much work, and the trade is monotone in both directions so there is no way to be wrong about it except by choosing a point nobody wanted.

The last essay is about the other kind, where both directions are failures — where the constant does not buy anything, it decides which of several true answers the computation returns, and at a low enough precision it produces a proof of something false from a comparison that was performed correctly.

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. -0.16-0.63-0.493.43-0.025-0.632.3-0.69-0.95-1.70.076-0.49-0.694.3-1.6-1.41.93.4-0.95-1.65.40.0141.63-1.7-1.40.0144.60.055-0.0250.0761.91.60.0553.1A, symmetric→-0.164.600004.65.92.500002.51-1.20000-1.24.90.3800000.385.80.2400000.242H = QᵀAQ, tridiagonal‖A − QHQᵀ‖/‖A‖1.1·10⁻¹⁵below the subdiagonal0worst eigenvalue movement7.1·10⁻¹⁵a similarity, so the spectrum is untouched — and every later step is O(n²) rather than O(n³)one reduction, then every iteration is cheapthe eigenvalues did not move The form that makes it affordable One Householder reduction, done once, turns every subsequent iteration of the eigenvalue algorithm from cubic to quadratic cost. It changes no answer at all, which is why it is easy to describe as an optimisation and wrong to. 12345678910⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 4.3·10⁻⁹worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁4.3·10⁻⁹worst Frobenius discrepancy4.3·10⁻⁹κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁ The best approximation there is The error of the best rank-k approximation is not bounded by the next singular value. It is equal to it. That is an unusually sharp theorem, and it makes the theorem itself usable as an independent check on the computation.

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 errorDeflationDrop toleranceEckart–YoungIncomplete factorisationLow-rank approximationPreconditioningThe QR algorithmTolerance