How tolerance works#
Finite-precision arithmetic rounds, so a correct implementation rarely returns the exact answer. The tolerance tells you how far from the exact answer a correct implementation can land, for the inputs you actually have. ExactRef works it out from the rounding steps any implementation of the operation has to take, so it is not a fixed number per dtype.
The examples on this page use plain Python floats (float64) unless they say otherwise.
One rounding leaves no room#
import exactref as er
# ExactRef's result for 0.1 * 3.0, with 0.1 in a float64 array
r = er.multiply(er.asarray([0.1]), 3.0)
# The exact product, rounded once to float64
print("rounded:", r.rounded.tolist())
# How far a correct implementation's output can be from rounded
print("tolerance:", r.tolerance().tolist())
# Python's own float multiplication
print("Python:", 0.1 * 3.0)
rounded: [0.30000000000000004]
tolerance: [0.0]
Python: 0.30000000000000004
A product of two floats is rounded once. Every correct implementation
rounds the same exact product the same way, so the tolerance is 0 and
Python's * gives exactly the reference.
The order of additions changes the answer#
# ExactRef's result for the float64 sum of 1e16, 1.0 and -1e16
r = er.sum([1e16, 1.0, -1e16])
# The exact sum, rounded once
print("rounded:", r.rounded.tolist())
# What a loop that adds in index order returns
print("normal:", r.normal().tolist())
# How far a correct implementation's output can be from rounded
print("tolerance:", r.tolerance().tolist())
# Could a correct implementation return 0.0? And 8.0?
print("0.0 covered:", r.covers(0.0).all())
print("8.0 covered:", r.covers(8.0).all())
rounded: 1.0
normal: 0.0
tolerance: 1.1102230246251574
0.0 covered: True
8.0 covered: False
The exact sum is 1.0. Adding in index order, 1e16 + 1.0 rounds back to
1e16, because float64 values that large are 2 apart, and the loop ends
at 0.0. An implementation that adds -1e16 first keeps the 1.0. Both are
correct, so the tolerance has to allow both: 0.0 passes and 8.0 does not.
The tolerance comes from the inputs#
# The tolerance of each sum, worked out from its own inputs
print("ten 0.1s:", er.sum([0.1] * 10).tolerance().tolist())
print("1e16, 1.0, -1e16:", er.sum([1e16, 1.0, -1e16]).tolerance().tolist())
ten 0.1s: 8.881784197001252e-16
1e16, 1.0, -1e16: 1.1102230246251574
Both sums are 1.0, but the second one's partial sums pass through 1e16, and each addition there can be off by about 1. A fixed tolerance per dtype would be far too loose for the first sum or would fail correct implementations of the second.
normal() is what an ordinary loop computes#
# Add ten 0.1s in index order, rounding to float64 after each addition
total = 0.0
for x in [0.1] * 10:
total += x
print("loop:", total)
# ExactRef's emulation of that loop for the same values
print("normal:", er.sum([0.1] * 10).normal().tolist())
loop: 0.9999999999999999
normal: 0.9999999999999999
normal() emulates, exactly, an implementation that adds in index order in
the compute format. Use it to see what an ordinary implementation gets, next
to the reference.
The compute format matters#
# ExactRef's result for the sum of 300 bfloat16 ones
r = er.sum(er.asarray([1.0] * 300, "bfloat16"))
# The exact sum, rounded once to bfloat16
print("rounded:", r.rounded.tolist())
# The tolerance when the additions happen in float32, the default here
print("tolerance:", r.tolerance().tolist())
# The tolerance when they happen in bfloat16
print("bfloat16 tolerance:", r.tolerance("bfloat16").tolist())
rounded: 300.0
tolerance: 0.0
bfloat16 tolerance: 580.0
For bfloat16 inputs, ExactRef assumes by default that the additions happen
in float32. float32 holds every partial sum here exactly, so the tolerance
is 0. An implementation that keeps its running sum in
bfloat16 gets stuck at 256, because bfloat16 values from 256 to 512 are 2
apart and 256 + 1 rounds back down. Judged as such an implementation,
with tolerance("bfloat16"), 256 is a correct answer:
# Is 256 correct for an implementation that adds in float32, the default?
print("float32 covered:", r.covers(er.asarray(256.0, "bfloat16")).all())
# And for one that adds in bfloat16, named by the second argument?
print("bfloat16 covered:", r.covers(er.asarray(256.0, "bfloat16"), "bfloat16").all())
float32 covered: False
bfloat16 covered: True
The roundings in the compute format are to nearest, whatever rounding an
operation is given: rounding sets the reference's rounding and the final
rounding into the output format. A named algorithm sets its
own.
Normal and any-order: three verdicts#
tolerance() covers normal accumulation, the way ordinary kernels add: in
index order, or in 2 or 4 interleaved running sums as SIMD code keeps them.
tolerance(standard="any_order") covers every order of additions that does
not overflow. Between them, an output gets one of three verdicts.
# 1e16 and -1e16, padded with zeros, then 96 ones
xs = [1e16, 0.0, 0.0, 0.0, -1e16, 0.0, 0.0, 0.0] + [1.0] * 96
# ExactRef's result for their float64 sum
r = er.sum(xs)
# The exact sum, rounded once
print("rounded:", r.rounded.tolist())
# The tolerance for normal accumulation, the default standard
print("normal tolerance:", r.tolerance().tolist())
# The tolerance for any order of additions that does not overflow
print("any-order tolerance:", r.tolerance(standard="any_order").tolist())
rounded: 96.0
normal tolerance: 11.102230246252248
any-order tolerance: 228.70594307278597
Three implementations of this sum:
# Adds the values one at a time, in index order
def index_order(xs):
total = 0.0
for x in xs:
total += x
return total
# Adds element i to running sum i % 8, then adds up the 8 sums
def eight_lanes(xs):
lanes = [0.0] * 8
for i, x in enumerate(xs):
lanes[i % 8] += x
return index_order(lanes)
# Buggy: counts xs[0] twice
def first_twice(xs):
return index_order([xs[0]] + xs)
for f in (index_order, eight_lanes, first_twice):
# The implementation's output
out = f(xs)
# Whether it passes the normal tolerance, and the any-order one
normal = r.covers(out).all()
any_order = r.covers(out, standard="any_order").all()
print(f"{f.__name__}: {out}, normal: {normal}, any_order: {any_order}")
index_order: 96.0, normal: True, any_order: True
eight_lanes: 72.0, normal: False, any_order: True
first_twice: 1e+16, normal: False, any_order: False
index_order is exact here: no worse than normal accumulation can be.
eight_lanes keeps 8 running sums, as 8-wide SIMD code does. Lanes 0 and 4 start at
1e16and-1e16and swallow every 1.0 added to them, so it is off by 24. That is correct, because some order of additions gives it, but worse than normal accumulation can be.first_twice adds the first element twice: a bug, and no order of additions explains it.
What a tolerance cannot tell you#
# A buggy sum that drops the last element. Does the tolerance pass it?
print("covered:", r.covers(index_order(xs[:-1])).all())
covered: True
Dropping the last 1.0 gives 95.0, which is inside the tolerance. A tolerance says whether rounding can explain a difference. A bug whose effect is smaller than the rounding allowed on those inputs passes, so test with inputs where the answer is clear-cut too.
Where to go next#
Reading complicated tolerances: whole tensors, overflow, infinities and NaN, several outputs, complex outputs, and kernels that round more than the definition.
Results: everything a result holds and every option of
tolerance().