Conditioning, Tolerances, and the Digits a Reproduction Can Demand

Two runs of the same analysis on two machines produce 0.31415926535 and 0.31415926538. Whether that is a reproduction or a discrepancy is not a matter of taste, and it is not answered by choosing a tolerance after seeing the difference. The tolerable difference is a property of the computation, derivable before either run: it is the width of the interval in which the true answer could lie given the conditioning of the problem and the arithmetic used to solve it.

The interval a single result occupies

Double precision carries 53 bits of significand, so the unit roundoff is $u = 2^{-53} \approx 1.11 \times 10^{-16}$. For a scalar operation that is the whole story: the result is the exact value rounded, so the relative error is at most $u$.

Composition changes that. A sequence of operations accumulates error, and the accumulated error is amplified by the sensitivity of the answer to perturbations in the inputs. The standard framework separates the two causes precisely: a backward error analysis shows that the computed result is the exact result of a slightly perturbed problem, and a condition number relates that perturbation to the answer.

For a linear system $Ax = b$, the normwise condition number is $\kappa(A) = |A| |A^{-1}|$. A backward-stable solver produces $\hat{x}$ satisfying $(A + \Delta A)\hat{x} = b$ with $|\Delta A| \le c, n, u, |A|$ for a modest constant $c$. The forward relative error then satisfies

$$\frac{|\hat{x} - x|}{|x|} \lesssim \kappa(A) \cdot c, n, u.$$

The product $\kappa(A) u$ is the relative width of the interval around the answer. For $\kappa = 10^{12}$ and $u \approx 1.1\times10^{-16}$, that width is about $10^{-4}$: four decimal digits of the answer are determined, and the rest are noise from the arithmetic, not information about the problem. Reporting twelve digits is then not a claim of precision but a display of roundoff.

This is why a reproduction tolerance has to come from the analysis rather than from convention. A common default in test suites is a relative tolerance around $10^{-6}$, which is generous for a well-conditioned direct solve with $\kappa \approx 10$ and far too tight for an ill-conditioned one; applying it uniformly means either failing correct results or accepting wrong ones.

Tolerance relative to what quantity

The choice of denominator matters as much as the magnitude. Three conventions appear in practice and mean different things.

Relative to the result. $|\hat{y} - y| / |y|$. Fails when $y$ is near zero, where the ratio blows up while the absolute difference is tiny. A computed value of zero and a true value of $10^{-18}$ give a relative difference of one.

Relative to the scale of the data. $|\hat{y} - y| / \max_i |y_i|$ over a table or vector. This is stable and is the natural choice for a residual or a set of coefficients, but it hides a badly computed small component inside a well-scaled total.

Mixed, with an absolute floor. $|\hat{y} - y| \le a + r|y|$ for an absolute tolerance $a$ and relative tolerance $r$. The form used by many numerical test harnesses, because it degrades gracefully near zero while still scaling for large values.

For vectors and matrices, a normwise comparison can pass while individual entries are badly wrong: two vectors differing by $10^{-3}$ in every component but $10^{-1}$ in one may have a small normwise difference if the vectors are large. Componentwise comparisons with a per-entry scale, as in $\max_i |\hat{y}_i - y_i| / (\epsilon + |y_i|)$, catch that at the cost of failing when an entry is near zero through cancellation rather than error.

The distinction has a concrete consequence for verification. If a step’s output is a residual vector that is expected to be near zero in exact arithmetic, a relative tolerance on the norm is meaningless; the right check is whether the residual norm is below a threshold derived from the conditioning and the arithmetic precision, which is what the standard termination criteria for iterative solvers express.

Where the tolerance budget comes from

A tolerance on a final output has to accommodate everything on the path, and the contributions do not simply add.

Input uncertainty. Measured data carries uncertainty from the instrument, and that uncertainty propagates through the model. If the condition number of the map from inputs to output is $\kappa$, a relative input uncertainty of $\delta$ produces a relative output uncertainty of order $\kappa \delta$. For a well-conditioned model and data accurate to four digits, the output cannot be more accurate than four digits, and no amount of extra precision in the arithmetic helps.

Discretization error. A numerical method that approximates a continuous problem has a truncation error that shrinks with resolution at a rate set by the method’s order. That error is not roundoff and does not follow the same scaling; it usually dominates. A finite-difference scheme of order $h^2$ on a grid with spacing $10^{-2}$ carries an error around $10^{-4}$ relative to the solution scale, which makes a tolerance derived from unit roundoff an order of magnitude too strict.

Iteration error. An iterative method stopped at a residual threshold leaves an error bounded by the condition number times the residual. The threshold is a choice, and it is recorded as part of the computation’s identity; two runs with different thresholds are different computations even if everything else matches.

Arithmetic error. The roundoff contribution, of order $\kappa u$ as above.

Implementation divergence. Different summation orders, vectorized kernels, fused multiply-add availability, and library versions change the rounding, producing differences at the level of the accumulated roundoff. This is the component most often observed when comparing two machines, and it is bounded by the arithmetic term times a small factor for the difference in order.

The usable tolerance is bounded below by the arithmetic and implementation terms and above by nothing in particular; the sensible value is the sum of the terms that cannot be eliminated, which is usually dominated by discretization or input uncertainty. Setting a tolerance tighter than the arithmetic noise produces a test that fails randomly; setting one looser than the model error produces a test that passes for wrong answers.

Casting tolerance as a statistical statement

Monte Carlo results need a different treatment, because the run-to-run variation is part of the estimate rather than a defect. A reported mean with a standard error $s$ occupies an interval of width proportional to $s$. Two independent runs of the same size produce means that differ by roughly $\sqrt{2} s$ on average, so a reproduction check that demands agreement within a fraction of $s$ will fail most of the time even when both runs are correct.

The defensible check is a two-sample comparison: the difference of means divided by the pooled standard error, compared against a critical value. For two runs with $n$ samples each, the standard error of the difference is $\sqrt{s_1^2 + s_2^2}$, and a difference within a couple of multiples of that is consistent with the same underlying expectation. With few replicates, the $t$ distribution rather than the normal is the right reference, and its critical values are noticeably wider.

The same logic applies to any output that is a statistic rather than a deterministic function. Reporting a tolerance in units of standard error, rather than as a relative difference, is the difference between a check that works and one that has a false failure rate nobody has measured.

import math

def within(a, b, atol, rtol):
    """Mixed tolerance: symmetric, scale-aware, safe near zero."""
    return abs(a - b) <= atol + rtol * max(abs(a), abs(b))

# For an ill-conditioned solve, derive rtol from the conditioning rather than guessing.
kappa, n = 1e12, 2000
rtol = kappa * n * 2.0 ** -53
print(f"rtol ~ {rtol:.3e}")        # ~1e-4: the interval the arithmetic can support
print(within(0.31415926535, 0.31415926538, 0.0, rtol))

Using max of the two magnitudes rather than abs(a) makes the predicate symmetric, so that within(a, b) and within(b, a) agree; asymmetry in a comparison routine produces disagreements between a forward and a reverse check that look like data problems.

The failure mode of tolerance theatre

A tolerance chosen to make a test pass is worse than no test, because it converts a visible disagreement into a certified match.

The pattern is recognizable: a comparison fails at $10^{-8}$, the constant is relaxed to $10^{-6}$, the suite goes green. Nothing has been learned about whether $10^{-8}$ was a real defect. The disciplined version of the same response is to ask which term in the error budget predicts $10^{-8}$, and if none does, to investigate the difference rather than the threshold.

The opposite failure is equally common in published work: reporting a result to a fixed number of digits regardless of conditioning, then treating any deviation in those digits as a failure to reproduce. A result whose condition number supports four digits is reproduced by a value agreeing to four digits, and a comparison that demands eight is measuring the other machine’s rounding mode rather than the science. Agreement to the supported precision, stated with the tolerance and its derivation, is the claim worth making, and it is the one a reader can audit. Content-Addressed Provenance Graphs for Derived Outputs records the tolerance with the node so that the acceptance criterion travels alongside the numbers it qualifies.