An accuracy that is a backward error

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 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:

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.

And the size is a third axis, which decides whether the ratio is a constant of the identity or of this matrix.

The slope is one at every size and the ratio is a seventh at every size. At n = 64, 128 and 256 the fitted slope is 0.986, 0.983 and 0.994, and the eighteen ratios across the three sizes and six tolerances span 0.09 to 0.21 with no trend in either variable — not falling as the matrix grows, not drifting as the tolerance tightens, just scattered inside a factor of two about an eighth.

So the ratio is a property of the identity rather than of the 256-square it was first measured on: whatever ‖A − A_H‖ is, the residual it produces on a particular right-hand side is about an eighth of it, and that eighth is what the difference between a matrix norm and its action on one vector costs. A quantity that refuses to move across two decades of size and ten of tolerance is the strongest form this collection’s habit takes — the hypothesis was a line of algebra, and eighteen measurements sit on the correct side of it with a slope of one.

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 — and is a choice of units rather than a fact about difficulty. The backward error was typed. So the prediction is a prediction:

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.

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. It is natural to expect that if the original matrix is symmetric the perturbation is not quite symmetric, since the (1, 2) and (2, 1) blocks are compressed independently. The next section measures that, and it is not what happens.

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.

The perturbation keeps the symmetry, and not by luck #

Two blocks compressed independently could discard different tails, so A_H − A could fail to be symmetric by something of the size of the tolerance — which would put a chosen backward error outside the class of the problem posed, and is exactly what a nearby problem of the wrong kind is about. It is the right worry and it is worth measuring rather than conceding.

On a symmetric 128-square, the perturbation’s relative asymmetry — ‖E − Eᵀ‖ over ‖E‖ — rises sharply as the tolerance falls: 6.0·10⁻¹³ at ε = 10⁻², then 7.1·10⁻⁹, then 1.0·10⁻⁴, and 7.4·10⁻³ at ε = 10⁻¹². Read on its own that is the defect arriving, and arriving fast.

Now read the absolute figure, ‖E − Eᵀ‖ over ‖A‖. It is 1.022·10⁻¹⁵, 1.021·10⁻¹⁵, 1.022·10⁻¹⁵ and 1.021·10⁻¹⁵ — constant to three digits across ten decades of tolerance, while the perturbation itself shrinks by nine orders of magnitude. That is not truncation. It is about nine times the unit roundoff, and it is the rounding of the compression arithmetic; the relative figure climbs only because its denominator is disappearing underneath it.

So the asymmetry is at the arithmetic’s floor at every tolerance, and the reason is structural rather than lucky. For a symmetric A the (2, 1) block is the transpose of the (1, 2) block, so the two have identical singular values — a tolerance therefore selects the same rank in both — and a truncated singular value decomposition of Bᵀ is the transpose of the truncated decomposition of B. Compressing the two independently produces transposed answers by construction. The word independently is true of the code and false of the result.

That settles the question this essay inherited from the field before it, and settles it in the strongest available direction. The backward error a caller types is not merely small and not merely of the same kind as the problem posed: it is a symmetric perturbation of a symmetric kernel matrix, to the working precision, at every tolerance. There is no residue of the wrong kind to account for.

It also reprices the repair the paragraph above used to suggest. Compressing one triangle and copying it does buy something — half the compression arithmetic, and an asymmetry of exactly zero rather than of 10⁻¹⁵ — but what it buys is arithmetic, not structure, and it should be argued for on that basis. A reader told that it removes a structural defect would be told something the measurement does not support.

The general habit is worth taking away from this more than the number. A relative measure whose denominator is the quantity being driven to zero will report a rising fraction whether or not anything is rising, and there is nothing wrong with the fraction — it is the right thing to look at when the denominator is fixed. Here it is not, and the two readings of the same six runs disagree by ten orders of magnitude about whether a defect exists. Printing both is what settles it, and printing one is how a phenomenon like this survives.

One consequence for the dial this essay is about. Since the perturbation is symmetric to the working precision, a caller solving a symmetric positive definite problem through a hierarchical representation is entitled to the symmetric backward-error theory rather than the general one — the perturbed matrix is symmetric, so its eigenvalues are real, its own condition number is the ratio the caller expects, and the bound in the section above is being applied to a matrix of the class it was derived for. That is not a small licence: a great deal of the analysis this collection quotes assumes symmetry, and the assumption is usually inherited rather than checked.

What is not established by any of this is positive definiteness. A truncation can move a small eigenvalue across zero, and nothing above says it cannot; the perturbation being symmetric makes the question well posed rather than answering it. At ε = 10⁻² the perturbation is 1.7·10⁻³ of the matrix, so a family whose smallest eigenvalue sits below that fraction of the largest has no guarantee left — which is a statement about κ and ε together, and is the same product the whole essay is about read in the other direction.

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.

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.

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.

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.

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.