Float equality tolerance for single and half precision

Question

Suppose the metric is

abs(a-b) <= rtol * max(abs(a), abs(b))

i.e. math.isclose with atol=0. The default for float64 is rtol=1e-9. numpy.allclose uses rtol=1e-5.

Both 1e-9 and 1e-5 make sense to me for float64, depending on application - former for stricter equality, latter for engineeresque "good enough". What I seek to know is, what are their float32 and float16 equivalents? How about 10 ** (floor(log10(2 ** (-(mantissa_bits + 1) / 2))) - 1) (see (3) below)?

Seeking references, or simply reasoning also works. Below discussion self-admits to lacking a thorough look into this matter, yet numeric computation is an entire field, so I figure there must be something of sort for the most rudimentary task.

---

math.isclose rationale

PEP 485 carefully explains it, including (1) the default rtol=1e-9; and, acknowledging that (2) no sensible default exists for atol, (3) ULP may be more appropriate depending on use case. It also describes alternate proposals to deal with atol issues. Further tolerance selection rationale is described in the Boost project.

For sake of this question, ignore atol, assume a != 0 && b != 0.

A discussion

Highlights from numpy.isclose vs math.isclose:

(1) atol=1e-8 is bad unless all input values are on order of unity. np.allclose(1e-9, 2e-9) is True.

(2) math.isclose author:

Note that the relative tolerance of 1e-8 was not only about half the precision of a float64, but also was large enough that there would be no difference in the results for various methods of comparison -- i.e. it really didn't matter which approach was taken -- again, fewer surprises.

(3) Proposed other-precision tolerances

Clarification

Relative error is by definition the goal. We seek a metric that is (1) robust, (2) order-of-magnitude independent, (3) readily configurable for tasks involving real-world data (known noise threshold). "Robust" defined as error variance matching data variance (so variability not induced by float limitations). rtol fails (1) near float epsilon, hence the need for atol.

Concerning defaults, it's easy to agree that generally, 1% is too weak, and $10^{-20}$% is too strong. A general purpose default can be derived via numeric simulations that sweep many common operations on different input sizes - repeated fft, squaring, square rooting, adding, multiplying. Example (full code at bottom):

x0, x1 = randn(size), randn(size)  # float64
out0 = cast(op(x0, x1), 'float16')
out1 = op(cast(x0, 'float16'), cast(x1, 'float16'))
err = abs(1 - out1/out0)

size = 100
float32 3.03e-08
float16 1.17e-04

size = 10000
float32 2.91e-08
float16 1.92e-05

This'd be done separately for GPUs per significant multithreading-induced error, hence different defaults. A user setting their own defaults is far more error prone. Those knowledgeable enough won't need defaults in the first place - but they also benefit per not having to set new thresholds for every new task, also having a reference (if new th is much greater or lower, think twice).

We'd also have to graph error vs number of operations, and it must be sufficiently stable up to a "reasonable" # of ops. I've just done some brief testing with my aforementioned operations, and it's indeed such. What's "reasonable" can be determined by surveying many real-world applications.

Example's code

from numpy.random import seed, randn
from numpy import sqrt, sign

for N in (100, 10000):
    print("size =", N)
    for prec in ('float32', 'float16'):
        # seed for reproducibility; generate data
        seed(0)
        x0, x1 = randn(N), randn(N)
        x0p, x1p = x0.astype(prec), x1.astype(prec)
        
        # float64 ground truth, cast to `prec` after computation
        x0 = x0**2 * sign(x0)
        x1 = sqrt(abs(x1))
        o = sum(x0 * x1).astype(prec)
        
        # fully lower precision computation
        x0p = x0p**2 * sign(x0p)
        x1p = sqrt(abs(x1p))
        op = sum(x0p * x1p)
        
        print(prec, "%.2e" % abs(1 - o/op))

Answer 1

Those proposed tolerances look fine, but in my (opinionated) view this is really a problem with no satisfying solution, as the comments also argue.

For most algorithms, the error bounds one gets look like $O(n^p\mathsf{u})$, where $\mathsf{u}$ is the machine precision (ulp, eps(1.0)) and $n$ is the problem size. Often, the worst-case $p$ is larger than the average-case $p$, for instance for summation $p_{worst}=1$ and $p_{average}=1/2$. For linear algebra algorithms, the worst-case exponents might be larger, for instance $p=2$, but typically the average-case ones are still linear or less.

In particular, that $n$ means that the expected error level is very problem-dependent, and it is tricky to decide on a sensible default without knowing $n$.

The idea behind $10^{-8}$ is that it is a quantity that figures to be (1) larger than any $O(n^p \mathsf{u})$ you might reasonably encounter, but (2) small enough to avoid that different numbers erroneously compare as equal.

The good thing about double precision is that you have enough digits to satisfy (1) and (2) at the same time. With fewer digits available in your number type, those requirements are in conflict: there simply is no rtol that satisfies both (1) for a wide variety of problems and (2).

This issue seems fundamental to me: double precision is a sensible default; if you want to go lower, you will start to feel cramped. If you want to save on computational costs by skimping on precision you need to turn your brain on and be careful about numerical errors, so that your computation does not produce garbage. If you don't understand the errors in your algorithm well enough to even decide a sensible default tolerance, you're going to have much bigger problems in the rest of your work.