1
$\begingroup$

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

enter image description here

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))
$\endgroup$
12
  • $\begingroup$ The proposal looks good to me. I'm not sure what more there is to add that isn't in the thread. It's also "just" a default value, anyone can of course pick what's right for them. I don't think any value of rtol or rtol can ever be bad. Absolute tolerance does exactly what it is supposed to. How else would you check if a value is close to zero? $\endgroup$ Aug 11, 2023 at 13:17
  • $\begingroup$ @njuffa "based on a maximum ulp difference" meaning what, exactly? What's the formula? $\endgroup$ Aug 12, 2023 at 8:26
  • $\begingroup$ @MikaelÖhman rtol is closely equivalent to measuring error in significant figures. atol, particularly 1e-8 for float64, is troubling since it's nowhere near zero, and gives false positives (as shown in (1)). I favor math.isclose's formulation in that it seems to closely resemble a conditional check - if max(abs(a), abs(b)) < eps: use_atol(). $\endgroup$ Aug 12, 2023 at 8:28
  • $\begingroup$ @njuffa Thanks, though I don't know about bound independence, compute errors grow faster at lower precision. I also don't follow your concern on relative error. It's indeed very useful mathematically (the analogue is relative Euclidean distance, often invoked in theorems) and anytime there's expected significant figures, and I don't see what contexts it wouldn't be useful in, aside near float epsilon. $\endgroup$ Aug 12, 2023 at 10:49
  • $\begingroup$ @njuffa It's a normalized error measure, essentially in percents. Percents are easily interpretable and configurable, if not essentially universal. I'm quite surprised you take this much an issue with it - I'll open a separate question if you're willing to elaborate. Note, it's simply abs(a-b)/abs(a) > rtol (simplified & rewritten), so distance between a and b that's adjusted for a's norm. Wavelet scattering theorems use relative Euclidean distance, sqrt(sum(abs(a - b)**2) / sum(abs(a)**2)). I'm not aware specifically on isclose. $\endgroup$ Aug 12, 2023 at 11:11

1 Answer 1

4
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Thanks. Could you elaborate on what "problem size" means (input size)? And how the $p$'s are determined. Concerning float32 and float16, there's certainly been significant work done in both, since PyTorch and TensorFlow default to 32 and support 16, and much work was put in by NVIDIA and other researchers in making even 16 work, which of course wouldn't be possible without tolerances. Deep neural networks involve a ton of computations, but supporting 16 requires also handling fewer compute. So I guess I should look there. I just expected someone in this network to have such references handy. $\endgroup$ Aug 13, 2023 at 10:34
  • $\begingroup$ (Off-topic) I'm actually not at all happy with NumPy's default atol=1e-8 with assuming "on order of 1". I've done a lot of signal processing work and have a strong sense of what qualifies as "zero" for data on order of 1 - it's at most 1e-15 for float64. 1e-8 is good if you're trying to test if an individual array - a set of measurements - is close to zero, rather than for comparing two arrays which are supposed to be float-equal. The current heuristic I'm thinking of is atol = median(abs(x))) / 1e15 - and just settle on one of a or b since if they're "equal" it won't matter $\endgroup$ Aug 13, 2023 at 11:20
  • $\begingroup$ Ongoing discussion for those interested $\endgroup$ Aug 13, 2023 at 12:55
  • $\begingroup$ @OverLordGoldDragon For summation, $n$ is the number of elements you wish to compute the sum of. For linear algebra, the matrix size. If you want to study more, Higham's Accuracy and stability of numerical algorithms is a good encyclopedic book on the classic algorithms. I know there has been a lot of work on single, half and even quarter precision recently; I agree that it is an area with a lot of potential, but still I think it is impossible to come up with a one-size-fits-all tolerance and careful analysis is needed. $\endgroup$ Aug 13, 2023 at 16:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.