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

An accuracy that is a backward error

Every backward error on this site is something an algorithm produced and somebody then measured. This one is a line in the program. Solving with a compressed matrix gives a residual that is the compression's own error, at a slope of 1.000 over ten decades, so the knob that sets the storage sets the backward error directly.
34 min read 25 figures Whose fault is itWhere the algebra stops being the arithmeticTwo routes to a number

Worth reading first: The exact answer to a nearby problem · The condition number is an amplifier · A block nobody can call sparse.

A backward error, everywhere else on this site, is a diagnosis. An algorithm ran, an answer came out, and the question was how far the problem would have to move for that answer to be exactly right. It is measured afterwards, it is a property of what the arithmetic did, and nobody chooses it.

This essay is about one that is chosen, in a line of code, before anything runs.

The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.994 and sit a constant 6.7× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸1.2·10⁻⁹backward error there1.3·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.5·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 1 The backward error of a solve against the error of the representation it used, at six accuracies over ten decades. The dashed line is where they would be equal.

The line of algebra

A hierarchical representation A_H is an approximation of A whose accuracy is set by a compression tolerance. Solve with it — exactly, by the recursion that never assembles anything — and the answer satisfies A_H x = b.

Ask what residual that answer has against the matrix that was actually wanted:

b − A x = A_H x − A x = (A_H − A) x .

There is nothing more to it. The residual is the representation’s error applied to the answer, so

‖b − Ax‖ ⁄ ‖A‖‖x‖ ≤ ‖A − A_H‖ ⁄ ‖A‖ ,

and the compression is the backward error. Not bounded by it, not related to it — it is the perturbation, written down in advance and applied on purpose.

The measurement

The algebra is one line and this collection’s habit is that a line of algebra is a hypothesis. Six accuracies, on a 256-square at κ = 24.39:

ε ‖A − A_H‖₂ ⁄ ‖A‖₂ ‖b − Ax‖ ⁄ ‖A‖‖x‖ ratio
10⁻² 1.94·10⁻³ 2.56·10⁻⁴ 0.13
10⁻⁴ 9.82·10⁻⁶ 1.54·10⁻⁶ 0.16
10⁻⁶ 1.28·10⁻⁷ 2.19·10⁻⁸ 0.17
10⁻⁸ 1.22·10⁻⁹ 1.27·10⁻¹⁰ 0.10
10⁻¹⁰ 1.21·10⁻¹¹ 1.43·10⁻¹² 0.12
10⁻¹² 1.20·10⁻¹³ 2.52·10⁻¹⁴ 0.21

The slope of one against the other, over ten decades, is 1.000. The ratio sits at about a seventh and does not drift, which is the difference between a norm of a matrix and a norm of that matrix applied to one particular vector — the inequality above is an inequality, and the residual is on the correct side of it at every point.

That is what makes the figure a check rather than an illustration. Two quantities that were both small would prove nothing; two quantities that move together decade for decade, on the correct side of a stated inequality, are the identity behaving.

The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.986 and sit a constant 5.7× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸5.7·10⁻¹⁰backward error there1.2·10⁻¹⁰representation at 10⁻¹²1.1·10⁻¹³backward error there2.3·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 2 The same identity at a quarter of the size, where both quantities are smaller and the line has not moved.

Both factors, known in advance

The identity this whole collection is organised around is

forward error ⪅ condition number × backward error ,

and it is normally used as an account of something that already happened. The condition number is a property of the problem, which can be estimated; the backward error is a property of what the algorithm did, which can be measured; and the product explains the answer that came out.

Here both factors are known before the solve starts. κ belongs to the problem. The backward error was typed. So the prediction is a prediction:

ε κ × backward measured forward error
10⁻² 4.7·10⁻² 2.55·10⁻³
10⁻⁴ 2.4·10⁻⁴ 1.86·10⁻⁵
10⁻⁶ 3.1·10⁻⁶ 2.71·10⁻⁷
10⁻⁸ 3.0·10⁻⁸ 1.59·10⁻⁹
10⁻¹⁰ 2.9·10⁻¹⁰ 1.73·10⁻¹¹
10⁻¹² 2.9·10⁻¹² 3.74·10⁻¹³

The bound holds at every point and over-predicts by a factor of about fifteen, which is the ordinary looseness of a worst-case statement applied to one right-hand side. The forward error tracks it decade for decade at a slope of 1.000.

And that licenses a sentence a code can act on:

Decide how many digits you want in the answer, divide by the condition number, and compress to that.

Nothing else on this site allows that sentence. Everywhere else the backward error is what it is — the algorithm is stable or it is not, and the analyst’s job is to find out which. Here it is a dial, and the dial is calibrated in digits of the final answer as soon as κ is known.

The forward error of a hierarchical solve, and κ times the backward error that was chosen for itForward error ⪅ condition number × backward error is the identity this collection is organised around, and it is normally used after the fact: the algorithm ran, somebody measured what it did, and the condition number explains the difference. Here both factors are known in advance — κ = 24.39 is a property of the problem, and the backward error was set at the top of the program — so the upper line is a prediction rather than an account. The measured error tracks it at a slope of 0.981 and sits 14× below it throughout, which is the usual looseness of a worst-case bound applied to one right-hand side and not a failure of the prediction. What the figure licenses is the sentence a code needs: decide how many digits the answer requires, divide by κ, and compress to that.10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹‖A − A_H‖ ⁄ ‖A‖, chosen‖x − x*‖ ⁄ ‖x*‖κ × the chosen backward errorthe error in the answerboth factors known firstκ24chosen at 10⁻⁸1.2·10⁻⁹predicted forward3·10⁻⁸measured forward1.6·10⁻⁹bound ⁄ measured19the amplifier, used forwardsfor once
Fig. 3 The identity used forwards for once: the amplifier’s prediction and the measured error, on one figure, with κ known before either was computed.
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. 4 The identity itself, from the field it organises: the two factors, and the fact that only their product is what anyone sees.

What the perturbation is, which the previous essays could not answer

the-exact-answer-to-a-nearby-problem establishes what a backward error means: the answer that came out is the exact answer to a problem some distance from the one posed. a-nearby-problem-of-the-wrong-kind is the awkward follow-up — the nearby problem a backward error promises may not be a problem of the kind that was posed at all. A symmetric matrix’s backward error may be attained only by an unsymmetric perturbation, and then the promise is about a problem nobody would have asked.

This case answers that question cleanly and in the good direction. The perturbation is A_H − A, and A_H is available: it is a hierarchical matrix, built from the same points, with the same partition and the same block structure. So the nearby problem is a kernel matrix of the same kind as the one posed, with slightly different entries.

More than that, it is a perturbation whose structure is known blockwise. Every dense block of the partition is exact and contributes nothing; the whole perturbation lives in the compressed blocks, and in each of those it is the discarded tail of a spectrum. If the original matrix is symmetric the perturbation is not quite symmetric — the (1, 2) and (2, 1) blocks are truncated independently — and that is a real and small defect worth naming rather than hiding, and one a symmetric implementation removes by compressing one triangle and copying.

So this is a structured backward error in the sense the structure field means, with the structure mostly preserved, and the mostly has a name.

Backward error with and without the structure, for Levinson and for elimination, n = 8Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 4.7·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 9.478·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 perturbation1.9·10⁻¹⁷Levinson, Toeplitz only1.8·10⁻¹²the ratio between them9.5·10⁴diagonal defect of the first0.89the number that is reportedand the number that is asked about
Fig. 5 The question this answers, from the structure field’s essay on it: whether the nearby problem a backward error promises is a problem of the kind that was posed.

Where the digits actually go

It is worth walking one case end to end, because the chain has three links and each of them loses something.

Ask for 10⁻⁶. The blocks are truncated at 10⁻⁶ of their own leading singular value; the assembled matrix comes out at 1.28·10⁻⁷ of ‖A‖ — a factor of eight better than asked, because the off-diagonal blocks carry less of the matrix’s norm than the diagonal ones do, which is the previous field’s measurement. The solve then returns an answer whose residual is 2.19·10⁻⁸, another factor of six better, because the perturbation applied to one vector is smaller than its norm. And the forward error is 2.71·10⁻⁷, which is the residual multiplied by something like the condition number.

Net: a tolerance of 10⁻⁶ produced an answer good to about 3·10⁻⁷, on this problem, and the coincidence that those two numbers are close is a coincidence. Two of the three factors happen to cancel the third at κ = 24. At κ = 10⁴ they do not, and the essay’s refusal is that case.

The general shape is worth extracting because it is the whole reason this essay is in the error field rather than the hierarchy one. There are three separate ratios between the number typed and the number obtained, they are computed from three unrelated things — where the matrix keeps its norm, what the perturbation does to one vector, and how sensitive the problem is — and a code that assumes any of them is one is a code that has guessed.

The accuracy asked for, against the accuracy obtained, for two kernels on one partitionε is applied to each block against that block's own largest singular value; the error is then reported against the whole matrix's norm. The two are not the same number and the dashed diagonal is where they would be. For 1/r the obtained error is 29 times smaller than the tolerance asked for; for log r on the identical partition it is 13 times smaller. Both curves are straight and parallel to the diagonal, so the knob does what a knob should — a decade in buys a decade out — and neither of them sits on it. The factor is how much of the matrix's mass lives off the diagonal, which is a property of the kernel; it is not the number of blocks, which is 112 and would put the curves on the other side of the diagonal.10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³ε asked for, per block‖A − A_H‖ ⁄ ‖A‖ obtainedasked = obtainedlog r1 ⁄ rtwo numbers, not one1/r, obtained at 10⁻⁸3.4·10⁻¹⁰log r, obtained at 10⁻⁸7.8·10⁻¹⁰1/r, obtained ⁄ asked0.034log r, obtained ⁄ asked0.078blocks in the partition112the tolerance is per blockand the error is per matrix
Fig. 6 The first of the three ratios, from the hierarchy field: the accuracy asked of each block against the accuracy the assembled matrix has, which is thirty times better and belongs to the kernel.
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. 7 And the third, from this field’s oldest essay: a residual and an error are two numbers and the ratio between them is the condition number.

What it means for the rest of the collection

A backward error that is chosen rather than discovered changes the shape of several arguments this site has already made, and it is worth saying which.

It makes the compression a member of a family. This collection already has three approximations that turn out to be backward errors rather than errors: iterative refinement’s residual, computed at a lower precision, perturbs the problem rather than the answer; a matrix-free Jacobian, formed by difference quotients, is the exact Jacobian of a nearby function; and a static pivot order, kept from an earlier member of a sequence, changes the matrix rather than the arithmetic. Each of those was found by analysis after the fact. This one is the first that was designed that way.

It changes what “approximate” means for a solver. A hierarchical solve is usually described as an approximate direct method, which sounds like a hedge between two categories. It is not: it is an exact direct method for a nearby problem, and the distance to that problem is a parameter. That is a completely standard object — every backward-stable algorithm on this site is one — and the only unusual thing is that the distance is large and deliberate rather than the size of the unit roundoff.

And it explains why the outer loop repairs it. The sequence field’s central finding is that an outer iteration recomputes its residual from the matrix at every step, so whatever the inner solve got wrong is measured again rather than carried. A perturbation to the matrix is exactly the kind of wrong that gets measured again — a preconditioned iteration on the true A, preconditioned by a solve with A_H, converges to the true answer however bad A_H is. That is the subject of the next essay in the hierarchy field, and this identity is the reason it works.

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. 8 The first member of the family, from the arithmetic field: a residual computed in the wrong precision, which perturbs the problem rather than the answer.
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. 9 And the sequence field’s finding that makes the repair work: eleven decades of inner accuracy landing an outer step at the same four digits, because the loop measures again rather than carrying.

The one thing this does not buy

A backward error known in advance is not a backward error that can be made small for free, and the temptation is to read it as one.

The cost of a smaller compression tolerance is more columns per block, and the second essay in the hierarchy field measures the exchange rate: about half a column per decade, per block, at a fixed geometry. So ten decades of backward error costs roughly five columns everywhere, which is why the storage figures across this whole phase move by a factor of two or three across ten decades of tolerance rather than by a factor of ten.

That is a very favourable exchange rate and it is the reason the whole format works. It is not a free one, and there are two places where it stops being favourable: a kernel with a length scale of its own, where the columns per decade is a function of the geometry’s size, and a tolerance below the working precision, where the count stops selecting columns and starts counting rounding error.

How many columns a decade of accuracy costs, measured and predicted, at q = 0.500The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 1. The measured curve is a straight line at 0.55 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 3.32 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 5.3× about the constant, which is the safe direction for a quantity you have to allocate storage from.0369121501122334455digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.55bound, a decade3.3rank at 10⁻⁸5bound at 10⁻⁸28q0.5the shape is rightand the constant is not
Fig. 10 The exchange rate, from the hierarchy field: columns against digits, straight, at about half a column a decade.
Two kernels on identical geometry, as the picture is scaled up at a fixed wavelengthTwo parallel segments of length L a distance 2L apart, so q = 0.50 at every point on this sweep and nothing about the arrangement changes. 1/r needs 6, 6, 6, 6, 6 columns — one number. cos(κr)/r at κ = 40 needs 12, 16, 22, 33, 53, growing like a power of the size with no size at which it stops, because what decides it is the number of wavelengths across the pair — 3.2, 6.4, 12.7, 25.5, 50.9 — and no ratio of lengths can see that. The segments are parallel rather than collinear on purpose: in one dimension cos(κ|x − y|) obeys an addition formula and is exactly rank two, which makes the collinear version of this look like a confirmation of the smooth case and is an artefact of the arrangement.-1012301020304050log₂ of the segment lengthcolumns above 10⁻⁸cos(κr) ⁄ r at κ = 401 ⁄ r, same pointsthe case with no answerq, held fixed0.51/r at L = 0.561/r at L = 86cos(κr)/r at L = 0.512cos(κr)/r at L = 853the geometry did not moveand the rank did
Fig. 11 And where it stops being favourable, from the same field: a kernel whose columns-per-decade is a function of how large the picture is.

Two routes to the perturbation

The identity above is checked one way in the table and there is a second way, which shares no arithmetic with the first and is worth taking because a single route has been wrong every time this collection has published one.

Route one is the residual. Solve, form b − Ax with the true dense matrix, and divide by ‖A‖‖x‖. That touches the answer and the original matrix and never looks at A_H.

Route two is the representation. Assemble A_H — which no real code does, and which is affordable at these sizes — subtract it from A, and take the norm. That touches neither the answer nor the right-hand side.

The two are the columns of the first table and their ratio is a seventh, held to within a factor of two across ten decades. They cannot agree exactly: one is a matrix norm and the other is that matrix applied to one vector, and the second is smaller for the same reason a matrix norm is a maximum.

There is a third route hiding in the second and it is worth naming because it is the only one a real code can afford. ‖A − A_H‖₂ is estimated here by a power iteration on the difference, which is forty products with a matrix that is never formed — and the check that it is allowed to be used at all is a comparison against the decomposition at a size where the decomposition is affordable. Two of the three numbers on that check agree to eight digits; the third, the norm of the difference, agrees to four, because the difference is the discarded tail of a truncation and its leading singular values are nearly equal, so a power iteration on it converges at the ratio of two numbers that are almost the same. That is a real limitation of the cheap route, it is written into the assertion’s tolerance rather than into a comment, and it is four digits below the smallest claim anything here rests on.

How close Hager's estimate is to the true κ₁, over 200 seeded 8×8 matricesFive bars. The estimate is exactly the true condition number on 81% of the sample and inside ten per cent on 87%; the worst underestimate in the whole sample returns 38% of the truth. The last bar is the matrix built to defeat it, at 7.7%, well below anything the random sample reached.each bar is a percentage — of the sample, or of the true condition numberexactly right80.5%inside 10%87.0%inside a factor of 287.0%worst in the sample, ×10037.7%the constructed matrix, ×1007.7%usually exactexact share0.81worst of the sample0.38the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 12 Estimating a norm without forming what it is a norm of, from the field that measured how well it can be done. The route this essay’s second column uses is of exactly that kind.
The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻¹⁴ and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 7.3·10¹⁰ apart, with the reported one at 6.92·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.01428425670849811210⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.1·10⁻¹⁰unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10¹³iterations drawn110the recurrence remembers every roundingand the stopping test is written in it
Fig. 13 And the collection’s standing warning about two numbers that look like the same number, from the iterative field.

The refusal

The claim under test is the one anybody would make after reading a compression tolerance in a manual: that compressing a matrix to ε gives an answer accurate to ε.

Its verdict is wrong blame rather than false, because at κ = 24 it is very nearly true, for reasons that have nothing to do with it being right. The refusal is therefore fed a case where the reason is visible: the same experiment with the matrix’s shift enlarged so that κ = 1.07·10⁴, at a tolerance of 10⁻⁶.

The representation’s error is 1.7·10⁻⁷. The residual is 2.9·10⁻⁸. And the answer’s forward error is 8.5·10⁻⁵ — a hundred times the tolerance that was asked for. The assertion that the answer is accurate to 10⁻⁶ is fed that number, and it fails.

Nothing went wrong. The compression did exactly what it promised, the solve did exactly what it promised, and the answer is exactly as accurate as the identity says it should be. What was wrong was the reading, and the reading is the one every user of every tolerance in numerical software makes at least once.

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 The same misreading in the arithmetic field, where the parameter is a mantissa rather than a tolerance: how many digits a computation keeps is not how many digits it was given.
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 1.4·10⁻¹³ at 53 bits to 0.00597 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 bits1.5·10⁶error ÷ residual at 16 bits2.3·10⁷error ÷ residual at 8 bits2.4·10¹⁰every convergence test passesand the answer is wrong in proportion to u
Fig. 15 And the sharpest version of it on the site: a residual at the unit roundoff on an answer with no correct digits at all.
The accuracy worth paying for in a preconditioner, against the condition number of the problemFour complete cost sweeps, each reduced to the ε that minimised the total work, against the κ of the problem it was run on. The shift that moves κ over 2.7 decades leaves every distance between every pair of clusters exactly where it was, so the ranks the partition asks for stay at 10, 9, 7, 5 throughout and nothing geometric is moving. The optimum walks 0.5 → 0.5 → 10⁻⁶ → 10⁻⁸ — six decades of accuracy, bought because the problem got harder and for no other reason. On the easiest problem the best hierarchical preconditioner in the comparison has rank one; on the hardest it is the tightest one on the sweep. How accurate an approximate inverse should be is a question with an answer, and the answer is not in the matrix's structure.10¹10²10³10⁴10⁵-9-7-5-3-11condition number of the problemlog₁₀ ε at the cheapest totalrank 1rank 1rank 6rank 5the knob answers to κκ = 210.5κ = 1.1·10⁴10⁻⁸rank at the easy end1rank at the hard end5ranks at 10⁻⁸, spread5the geometry did not moveand the right accuracy did
Fig. 16 What follows from taking the identity seriously, from the hierarchy field: the accuracy worth paying for, against the condition number, and nothing else.

The identity, at other settings

The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.983 and sit a constant 7.4× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.98representation at 10⁻⁸1.4·10⁻⁹backward error there1.4·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.3·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 17 The same identity at half the size, where both quantities are smaller and the slope is unchanged.
The forward error of a hierarchical solve, and κ times the backward error that was chosen for itForward error ⪅ condition number × backward error is the identity this collection is organised around, and it is normally used after the fact: the algorithm ran, somebody measured what it did, and the condition number explains the difference. Here both factors are known in advance — κ = 20.88 is a property of the problem, and the backward error was set at the top of the program — so the upper line is a prediction rather than an account. The measured error tracks it at a slope of 0.971 and sits 16× below it throughout, which is the usual looseness of a worst-case bound applied to one right-hand side and not a failure of the prediction. What the figure licenses is the sentence a code needs: decide how many digits the answer requires, divide by κ, and compress to that.10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹‖A − A_H‖ ⁄ ‖A‖, chosen‖x − x*‖ ⁄ ‖x*‖κ × the chosen backward errorthe error in the answerboth factors known firstκ21chosen at 10⁻⁸1.4·10⁻⁹predicted forward2.9·10⁻⁸measured forward1.6·10⁻⁹bound ⁄ measured26the amplifier, used forwardsfor once
Fig. 18 And the amplifier’s prediction at the same size, where κ has barely moved and the gap to the measurement has not.
The accuracy asked for, against the accuracy obtained, for two kernels on one partitionε is applied to each block against that block's own largest singular value; the error is then reported against the whole matrix's norm. The two are not the same number and the dashed diagonal is where they would be. For 1/r the obtained error is 29 times smaller than the tolerance asked for; for log r on the identical partition it is 13 times smaller. Both curves are straight and parallel to the diagonal, so the knob does what a knob should — a decade in buys a decade out — and neither of them sits on it. The factor is how much of the matrix's mass lives off the diagonal, which is a property of the kernel; it is not the number of blocks, which is 112 and would put the curves on the other side of the diagonal.10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³ε asked for, per block‖A − A_H‖ ⁄ ‖A‖ obtainedasked = obtainedlog r1 ⁄ rtwo numbers, not one1/r, obtained at 10⁻⁸3.4·10⁻¹⁰log r, obtained at 10⁻⁸7.8·10⁻¹⁰1/r, obtained ⁄ asked0.034log r, obtained ⁄ asked0.078blocks in the partition112the tolerance is per blockand the error is per matrix
Fig. 19 The first of the three ratios between the number typed and the number obtained, which is a property of where the matrix keeps its norm.
How many columns a decade of accuracy costs, measured and predicted, at q = 0.333The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 2. The measured curve is a straight line at 0.50 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 2.11 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 4.3× about the constant, which is the safe direction for a quantity you have to allocate storage from.0369121507142128digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.5bound, a decade2.1rank at 10⁻⁸4bound at 10⁻⁸18q0.33the shape is rightand the constant is not
Fig. 20 What the chosen backward error costs in storage, from the hierarchy field: about half a column a decade, per block.
The residual of a formatted Cholesky, against the accuracy its blocks were compressed atA leaf of 16 means 34 truncations inside the factorisation, and the ranks of its blocks run to 4, 7, 9, 12, 14 across the sweep. The residual tracks the tolerance at a slope of 1.027 over ten decades: the knob does what a knob should, and nothing in the approximate arithmetic bends it. The ratio to the representation's own error stays at 0.93, 1.83, 0.88, 0.91, 0.91 — one excursion above one, at 10⁻⁴, where the tolerance sits between two singular values of one block and the rank it takes is a rounding of a decision rather than a measurement.10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²ε the blocks were compressed at‖A − LLᵀ‖ ⁄ ‖A‖the factorisationthe knob, on a factorisationtruncations34residual at 10⁻⁴2.2·10⁻⁵residual at 10⁻¹⁰9.9·10⁻¹²slope1worst ratio to the representation1.8ten decadesand a slope of one
Fig. 21 The same knob applied to a factorisation rather than to a solve, tracking at a slope of one over ten decades.
‖A − LLᵀ‖ ⁄ ‖A‖ for a Cholesky computed in the format, against how many truncations it performedThe recursion is the textbook one and every step stays in the format: the off-diagonal factor is exact, because a triangular solve against one factor of a rank-k block leaves a rank-k block, and the Schur complement is a rank-k addition to a hierarchical matrix, which is the only approximate step there is. So the number of truncations is a count — one per admissible block in the subtree being updated, at every level — and it runs from 0 at a leaf of 128 to 98 at a leaf of 8. The residual goes from 3.856·10⁻¹² to 1.035·10⁻¹¹, and all of that movement is the representation, whose own error is the upper curve and moves by the same factor: the ratio between them is 1.00, 0.93, 0.93, 0.91, 0.86, never above one. A hundred approximate operations contributed nothing measurable to the answer.02040608010010⁻¹²10⁻¹¹10⁻¹⁰truncations performed by the factorisationrelative errorthe representation's own error‖A − LLᵀ‖ ⁄ ‖A‖no decomposition without its residualresidual, 0 truncations3.9·10⁻¹²residual, 98 truncations10⁻¹¹representation, deepest1.2·10⁻¹¹residual ⁄ representation0.86levels, deepest5a hundred approximate stepsand an exact-looking factorisation
Fig. 22 And the measurement that says the approximate arithmetic inside that factorisation does not add a second term to the budget.
What a hierarchical preconditioner costs, against how accurate it is, at κ = 1.53e+3The iteration is conjugate gradients on the true matrix and the preconditioner is a solve with the representation, so what is measured is how much of the compression the outer loop repairs and what the repair costs. The count falls monotonically — 17, 17, 8, 4, 2, 2 steps against 105 with no preconditioner — and it keeps falling all the way to the tightest accuracy, which is the number a solver prints and the one everybody watches. The total work, one dense product with A plus one solve with A_H per step, is the other curve, and its minimum is at ε = 10⁻⁶ where the preconditioner has rank 6. Past that point every extra digit in the preconditioner is bought and not needed: at the tightest setting the run costs 1.34 times the cheapest, having taken 0 fewer steps to get there.-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 105 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ1527steps, no preconditioner105cheapest ε10⁻⁶its rank6tightest ⁄ cheapest1.3the count is what is printedand the cost is what is spent
Fig. 23 What the identity licenses a code to do with the knob, from the hierarchy field, once κ is known.
Multiplications in a hierarchical solve and in a dense factorisation, and where they crossBoth counted rather than estimated: the recursion carries a counter and reports what it actually did. The dense factorisation is n³/3 and rises at exactly three per doubling. The hierarchical solve rises at 2.25, 1.95, 1.78 — falling, because the cost is n log²n and its exponent is on its way to one. The two cross between n = 128 and n = 256: below it the format is the more expensive way to solve the system, at 2.36 times the dense count, and at n = 512 it is 3.4 times cheaper. Every point returns an answer at a backward error of about 2.4·10⁻¹⁴, so the comparison is between two ways of getting the same thing.567891010⁴10⁵10⁶10⁷10⁸10⁹log₂ nmultiplicationsdense factorisation, n³⁄3the recursion, countedwhere the format starts payingratio at n = 642.4ratio at n = 5120.29exponent, first doubling2.2exponent, last doubling1.8backward error2.4·10⁻¹⁴cheaper is a sizenot a property
Fig. 24 And what a tight backward error costs in arithmetic, which is where the crossover goes when more digits are asked for than the problem needs.
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. 25 The question of whether a nearby problem is a problem of the right kind, from the structure field, which this case answers in the good direction.

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. 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. the smallest perturbation of any kind — 1.89·10⁻¹⁷the smallest Toeplitz one — 2.65·10⁻¹³both exact for the same x̂smallest of any kind1.9·10⁻¹⁷smallest Toeplitz one2.7·10⁻¹³the price of the constraint1.4·10⁴diagonal defect, unconstrained0.97an exact answer to a nearby problemof a kind nobody posed A nearby problem of the wrong kind A good algorithm returns the exact answer to a nearby problem. A hundred and eighteen essays have measured the distance and not one has asked what the nearby problem looks like. On a Toeplitz system it is a rank-one matrix that is constant along none of its diagonals — and the smallest one that is Toeplitz is two and a half million times larger.

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.

AdmissibilityBackward errorBackward stabilityCondition numberForward errorHierarchical matrixResidualStructured perturbationToleranceWoodbury identity