The error the method already knows

A method that cannot be trusted is usually described as one whose error is unknown. That description is almost always right, and on the exponential’s Taylor series it is wrong. The error of that sum is not unknown. It is computable from two numbers the summation loop is already holding, it costs one division, and no implementation returns it.

The two numbers are the norm of the largest term the series produced and the norm of the sum those terms add up to. Their ratio is the cancellation — how much bigger the intermediates were than the answer they collapsed to — and it is already on the machine because a series that stops when its terms have become negligible against the running sum has to compute both quantities to decide when to stop. The claim of this essay is that this ratio, multiplied by the unit roundoff, tracks the relative error the sum actually achieved, well enough to act on, over the entire range in which the method goes from excellent to worthless.

The matrix is Moler and Van Loan’s 2×2, with eigenvalues −1 and −17, and the reason to use it is that its exponential is available in closed form, so the error being tracked is a measured quantity rather than a second computation’s disagreement with a first. The scalar in front of it is the one knob: the argument here is not about a family of hard problems but about one matrix, scaled, which is the transformation that changes nothing about the problem and everything about the method.

The estimate is a division, and the numbers are already in hand #

At the full matrix the cancellation reads 5.42·10⁶. The unit roundoff of a double is 1.11·10⁻¹⁶, and the product of those two numbers is 6.02·10⁻¹⁰. The relative error the series actually returned, against the closed form, is 5.20·10⁻⁹.

That is a factor of 8.6 out, and the reason it counts as a good reading rather than a poor one is the size of the thing being estimated. An error of 5·10⁻⁹ has to be distinguished from an error of 10⁻¹⁵ and from an error of 10⁻¹, and a quantity that lands within a factor of nine of it has answered that question completely. The estimate is being asked how many digits survived, and it says a little over nine where a little over eight survived.

The reasoning behind it is the oldest one in the subject, and the digits that cancel were never wrong: the subtraction that destroys them is exact, and what is left after it is the rounding the inputs carried in. A term of norm 1.4·10⁷ is stored with an absolute error of about 1.4·10⁷ times the unit roundoff, so the sum inherits that absolute error whatever it cancels down to. Dividing by the answer’s own norm turns it into a relative error, which is exactly the ratio above times the unit roundoff. Nothing in that argument is specific to the exponential, or to matrices, or to this matrix.

What is specific is that the ingredients are free. The loop forms the term norm to compare against the sum, and it keeps the running maximum for no other reason than that it is one comparison. Turning those into an error estimate is one division and it happens once, at the end, against a sum that cost 77 matrix products to build.

Where the series is a good method, the estimate says so too #

An estimate that only fires on failure is a warning label. This one reads correctly where the method is fine, which is what makes it a number rather than an alarm.

At a quarter of the matrix the cancellation is 15.2, the prediction is 1.69·10⁻¹⁵, and the measured error is 1.87·10⁻¹⁵ — a ratio of 1.11. The series at this setting is not merely acceptable, it is the better of the two methods on the page: scaling and squaring on the same matrix returns 1.65·10⁻¹⁴, an order of magnitude worse. The estimate says the sum is good to fifteen digits, and it is good to fifteen digits.

Halving the matrix instead of quartering it moves both quantities together. The cancellation becomes 955, the prediction 1.06·10⁻¹³, the measurement 1.60·10⁻¹³ — a ratio of 1.51, and the ordering against scaling and squaring has flipped by a factor of four. Between these two settings the series stops being the better method, and the estimate crosses the same line at the same place without being told where the line is.

Six settings, and a factor of eight of scatter #

The whole slider, with the prediction beside the measurement:

The measured error moves from 1.87·10⁻¹⁵ to 0.107, which is 5.8·10¹³ — a shade under fourteen orders of magnitude. The ratio of prediction to measurement moves from 1.11 to 8.65, which is a factor of 7.8. That is the whole result: the scatter is a factor of eight against a signal of nearly fourteen decades, so the estimate resolves the answer to within a digit of how many digits are left, everywhere on the range.

It is worth saying what the alternative is. A caller with no reference answer and no estimate has ‖A‖, which moves by a factor of eight across that table while the error moves by 10¹⁴ — so the norm orders the settings correctly and says nothing about the magnitude. The number of terms behaves the same way: 36, 51, 77, 90, 103 and 127, a factor of 3.5 against fourteen decades. Neither of the two numbers a caller would reach for first carries the information, and the ratio that does is formed from the same loop.

The setting between them that refuses to be drawn #

There is a gap in that table between 0.5 and 1, and the obvious way to fill it is at 0.75. The figure declines to draw it, and the refusal is a measurement rather than an inconvenience.

The generator carries a two-sided assertion, because the slider crosses a boundary: at or below half the matrix it requires the cancellation to stay under 10⁴ and the series to be within a hundred of scaling and squaring, and above half it requires the cancellation to exceed 10⁶ and the series to be at least 10⁴ behind. At 0.75 the cancellation is 6.99·10⁴ and the assertion fails on exactly that number — seven times above the small-end limit and fourteen times below the large-end one.

So the boundary between the two regimes is not a point, and the figure’s own refusal is what says so. Between half the matrix and three-quarters of it the series has lost nearly five digits and is still returning 1.69·10⁻¹¹, which is neither the regime where it competes nor the regime where it fails. The estimate has no such difficulty: it reads 7.76·10⁻¹² against a measured 1.69·10⁻¹¹, a ratio of 2.18, which is its closest reading at any setting above half the matrix. An assertion written to separate two regimes rejects the setting between them; a number written to interpolate does not. That is the difference between a test and an estimate, and it is why a library should return the second rather than only running the first.

Where the estimate stops being a number #

Every quantity on this site has a range, and this one’s is sharp enough to state as a rule.

At twice the matrix the series has no correct digits: the peak term is 2.39·10¹⁴, the answer has norm 1.01, and the measured relative error is 0.107. The prediction is 2.62·10⁻², a ratio of 4.10, so the estimate is still working — but it is working at the edge of what it can say, because a prediction of 0.026 is a prediction that between one and two digits remain, and the arithmetic that produced it has about that many left.

Past that point the estimate saturates, and the mechanism is worth naming because it is not a weakness in the reasoning. The denominator of the cancellation ratio is the computed sum, not the true answer. While the sum still has digits, its norm is the answer’s norm and the ratio means what it says. Once the sum is pure rounding noise, its norm is the noise’s norm, and the ratio stops growing while the error keeps going.

Widening the eigenvalue spread on the same 2×2 shows it happening. At eigenvalues −1 and −45 the computed sum has norm 1.14·10⁴ against a true 2.575, so the denominator is four orders too large; the cancellation reads 1.09·10¹⁵ and the error is 3.98·10³. At −1 and −60 the cancellation is 2.41·10¹⁵ and the error 6.47·10⁹. At −1 and −100 the cancellation has fallen to 1.07·10¹⁵ while the error has risen to 2.47·10²⁷. The ratio is no longer monotone in the difficulty, because the quantity it divides by has stopped being the answer.

The rule that follows is one-sided and is the useful direction. Over 3,040 matrices built as V diag(−a, −b) V⁻¹ — eight integer eigenvector matrices, a from 0.5 to 5, b from 8 to 45 — the estimate reads 10⁻² or worse on 1,081 of them, and on every one of those the measured error is above 2.63·10⁻³, with none below 10⁻³. The estimate never calls an answer worthless when it is good. It can understate how bad a bad answer is, by any amount; it cannot report a bad answer as a good one. That is the reverse of the one-sidedness the estimator every library ships for a condition number has, whose error runs in the flattering direction; this one runs towards caution.

It is not a property of this matrix #

A relation demonstrated on one 2×2 that has been scaled seven ways is a relation demonstrated once. The sweep above is also the check, and it was run for that.

On the 1,950 of those 3,040 matrices where the series cancels at least a thousand to one and the estimate predicts that at least two digits survive — the region where the claim is being made — the ratio of measurement to prediction has a median of 1.26 and lies inside a single decade on 99.4% of them. Its extremes are 0.082 and 17.9. Predictions in that set span 1.1·10⁻¹³ to 10⁻², which is eleven decades of estimate against a factor of 218 of spread in how well it does.

The check was repeated on matrices that are not 2×2 and not built from the same construction. Taking normal matrices at n = 4, 6 and 8 — eigenvector matrices that are exactly orthogonal, so the reference exponential carries no conditioning of its own — with eigenvalues spread from −0.5 down to as far as −40, and keeping the 72 cases where the series cancels at least a thousand to one and the estimate predicts that at least two digits survive, the ratio runs from 0.84 to 3.97 with a median of 1.74. Every one of the 72 is inside a factor of four.

Two things follow. The relation is not an artefact of the 2×2 or of non-normality, which is the obvious objection, since these matrices are normal and the estimate holds tighter on them than on the sweep. And the scatter is not noise around a fitted constant: nothing has been fitted. The predictor is cancellation times the unit roundoff, with no coefficient in front of it, and it is not a re-derivation of anything the method already returned — it is a quantity the method computes, discards, and could have handed back.

The repair has the same property from the other side #

The method that replaces the series scales the matrix down, evaluates a Padé approximant where nothing cancels, and squares the result back s times. The essay that priced both of its costs established that the parameter has a floor rather than a direction: too few squarings leaves the approximant outside its range, too many multiply its rounding by roughly 2ˢ. What that essay measured was the shape. What is measurable beyond it is the value, and the value on the right-hand arm is set by s alone.

At ‖A‖ = 1.2 the error at s = 0 is 4.01·10⁻¹⁵ and at ‖A‖ = 12.2 it is 5.55·10⁻¹¹ — the left wall moves by four orders of magnitude across a tenfold range of the norm, which is what a wall being set by the matrix looks like. On the right arm the same twelvefold range moves nothing. At fifteen squarings the relative error is 3.15887·10⁻¹² at μ = 1, 2 and 4, 3.15882·10⁻¹² at μ = 8 and 3.15912·10⁻¹² at μ = 12: five significant figures of agreement over a matrix whose norm has grown by a factor of ten.

The spread across those five matrices falls from a factor of 13,850 at s = 0 to 4.40 at s = 2, 1.07 at s = 5, 1.004 at s = 11 and 1.0001 at s = 15. Changing the dimension does not disturb it either: at s = 15 the error reads 3.158871·10⁻¹² at n = 4 and 3.158865·10⁻¹² at n = 12, six significant figures across a threefold change in size. Past the optimum the matrix has stopped being a variable in the problem.

And the number the arm converges to is the one the amplification argument predicts. Dividing the measured error by 2ˢ times the unit roundoff, over 250 combinations of dimension, norm and squaring count from s = 6 upward, gives a quotient between 0.100 and 3.115 with a median of 0.733 — every one of the 250 inside a single decade of one. For s ≥ 11 it tightens to between 0.470 and 1.126. The right-hand arm is 2ˢ·u, with an order-one constant and no fitting, and s is a number the caller chose.

So the second half of the standard method carries an error estimate as well, and it is cheaper than the first: it needs no division, because both of its inputs are constants of the run. A code that takes three squarings more than it needs is spending a factor of about eight, and it can say so before the matrix is looked at. That is a different statement from the one the amplification argument already makes — the error grows like 2ˢ is a rate, and the error is 2ˢ·u to within a factor of three is a value.

What follows for a caller #

Report the cancellation. It is one division at the end of a sum that cost seventy-seven matrix products, and it converts an answer into an answer with a bar on it. This is this collection’s standing habit applied to a method rather than to a parameter: print the quantity that decided something beside the thing it decided.

A convergent series is still not a method, and the estimate does not rescue it. At the full matrix the series returns 5.2·10⁻⁹ where scaling and squaring returns 8.5·10⁻¹⁴, and knowing that the first number is 5·10⁻⁹ does not make it 10⁻¹⁴. What the estimate buys is that the failure is reported at the time rather than discovered later, which is the difference between a wrong answer and a wrong answer that is labelled.

Report s and the scaled norm, and the arm they land on. A caller who is told the squaring count can price the right-hand cost without a reference answer, because that cost is 2ˢ·u and does not depend on the matrix. Reporting the parameter is worth as much here as reporting the estimate, which is not the usual relationship between the two.

Do not read either number as a bound. The series estimate is out by up to a factor of eighteen on the sweep above, and the squaring estimate by up to a factor of three; both are quantities to act on rather than guarantees, in the same way that a backward error chosen in advance is a promise and a forward estimate is not.

And the estimate exists because the intermediates were formed. A method that never builds a large intermediate has nothing to divide — the route that computes eᴬb without ever forming the matrix has no cancellation ratio to report, and needs a different stopping test for exactly that reason. The estimate here is a consequence of the defect it measures, which is why it is available on the method that needs it and not on the method that does not.

The refusal published at the top of this page is the second half stated as something that must fail. Squaring is exact in the algebra, so more squarings can only take the approximant further into its range, and a parameter whose costs sit at both ends has no such monotonicity in floating point. The sweep is run, the best value is taken, and the claim that the last value is no worse than the best is fed to an assertion that has to reject it: 6.88·10⁻¹⁶ at s = 2 against 3.16·10⁻¹² at s = 15, a factor of 4,593 bought with thirteen extra matrix products. It is the same object as a rounding applied repeatedly to its own output, and the reason it is published rather than asserted in a sentence is that a claim about what a machine does with an exact identity is exactly the kind of claim that has to be run.

And it is worth ending on which half of the subject this belongs to. The definition of a matrix function is not a method, and neither is a series that converges. What separates the two methods on this page from each other is not that one has an error and the other does not; both have errors, both errors are large in some regime, and both are computable from what the run already holds. The difference is that one of them reports its number and neither of them was written to.