Reading complicated tolerances#
How tolerance works covers single numbers. Real checks run over whole tensors and meet overflow, infinities, NaN, several outputs, complex numbers and kernels that round in their own way. This page shows how to read the result in each case.
A whole tensor at once#
covers() returns a bool array with the output's shape, so you can count
the verdicts and find the elements behind them. Here are 26 rows, each the
sum from the three verdicts
with 0 to 200 ones, computed by the 8-lane implementation:
import exactref as er
# 26 rows of 208 floats: 1e16 and -1e16 padded with zeros, then n ones
# (n = 0, 8, ..., 200), then zeros
rows = [[1e16, 0.0, 0.0, 0.0, -1e16, 0.0, 0.0, 0.0] + [1.0] * n + [0.0] * (200 - n)
for n in range(0, 201, 8)]
# Adds element i to running sum i % 8, then adds up the 8 sums in order
def eight_lanes(xs):
lanes = [0.0] * 8
for i, x in enumerate(xs):
lanes[i % 8] += x
total = 0.0
for lane in lanes:
total += lane
return total
# Its output for every row, as one ExactRef array
out = er.asarray([eight_lanes(row) for row in rows])
# ExactRef's result for the sum of each row
r = er.sum(rows, axis=1)
# Per row, whether the output passes the normal tolerance, as Python bools
normal = r.covers(out).scalars()
# Per row, whether it passes the any-order tolerance
any_order = r.covers(out, standard="any_order").scalars()
# Per row, True if the output passes the any-order tolerance only
worse = [a and not n for n, a in zip(normal, any_order)]
# Count the rows with each verdict
print("no worse than normal:", sum(normal))
print("correct but worse:", sum(worse))
print("wrong:", sum(not a for a in any_order))
no worse than normal: 6
correct but worse: 20
wrong: 0
6 rows are no worse than normal accumulation, 20 are correct but worse, and none is wrong. To see why a row landed where it did, read its numbers:
# Row 12, the one with 96 ones
i = 12
# Its output from eight_lanes
print("output:", out.scalars()[i])
# Its exact sum, rounded once
print("rounded:", r.rounded.scalars()[i])
# How far a correct output can be from rounded under normal accumulation
print("tolerance:", r.tolerance().scalars()[i])
output: 72.0
rounded: 96.0
tolerance: 11.102230246253484
Overflow, infinities and NaN#
The tolerance bounds finite outputs only. Three more methods handle the rest, so a check needs no special cases:
finite_range(): the lowest and highest finite output a correct implementation can produce.may_be("+inf"),may_be("-inf"),may_be("nan"): where a correct implementation may produce that value instead.covers(output): uses both.
Near the top of the range, the tolerance stays finite and the possible infinity goes in a flag:
# ExactRef's result for a float32 sum near the top of float32's range
top = er.sum(er.asarray([3.4028234e38, 2e31, -1e31], "float32"))
# The exact sum, rounded once to float32
print("rounded:", top.rounded.tolist())
# How far a correct finite output can be from rounded
print("tolerance:", top.tolerance().tolist())
# The lowest and highest finite output, (low, high), as plain floats
print("finite range:", [end.tolist() for end in top.finite_range()])
# Whether a correct implementation may return +inf instead
print("may be +inf:", top.may_be("+inf").tolist())
# Does +inf pass the check?
print("+inf covered:", top.covers(er.asarray(float("inf"), "float32")).all())
rounded: 3.4028234663852886e+38
tolerance: 4.056481920730334e+31
finite range: [3.4028230607370965e+38, 3.4028234663852886e+38]
may be +inf: True
+inf covered: True
The reference is the largest float32 value. An implementation that adds
2e31 first overflows to infinity, and that is correct too.
A NaN tolerance means the reference itself is not finite, so there is no distance to measure:
# ExactRef's result for a float32 sum of about 6e38, beyond float32's range
inf = er.sum(er.asarray([3e38, 3e38], "float32"))
# The reference is +inf, so the tolerance is NaN
print("rounded:", inf.rounded.tolist())
print("tolerance:", inf.tolerance().tolist())
# The lowest and highest finite output
print("finite range:", [end.tolist() for end in inf.finite_range()])
# Whether a correct implementation may return +inf
print("may be +inf:", inf.may_be("+inf").tolist())
rounded: inf
tolerance: nan
finite range: [nan, nan]
may be +inf: True
Here no finite output can occur, so the range is empty (both ends NaN) and only +inf passes. But a reference can be infinite while correct implementations still return finite values:
# ExactRef's result for a float32 sum whose exact value rounds to +inf
near = er.sum(er.asarray([3.4028234e38, 6e30, 6e30], "float32"))
# The reference is +inf, so the tolerance is NaN
print("rounded:", near.rounded.tolist())
print("tolerance:", near.tolerance().tolist())
# What a loop that adds in index order returns
print("normal:", near.normal().tolist())
# The lowest and highest finite output that can still occur
print("finite range:", [end.tolist() for end in near.finite_range()])
# Does the largest finite float32 value pass the check?
print("largest finite covered:",
near.covers(er.asarray(3.4028234663852886e+38, "float32")).all())
rounded: inf
tolerance: nan
normal: 3.4028234663852886e+38
finite range: [3.4028230607370965e+38, 3.4028234663852886e+38]
largest finite covered: True
The exact sum is just past the point where float32 rounds to infinity. Added
in index order, each 6e30 is less than half the spacing of float32 values
up there, so both are lost and the loop returns the largest finite value.
finite_range() keeps those finite outputs, so covers() accepts them.
A final rounding that cannot overflow (toward zero, a directed rounding on
its inward side, or overflow="saturate") gives the largest finite value.
The reference is then that value, and the tolerance is measured from it:
# Three float16 values whose exact sum, 70001, is past float16's largest, 65504
x16 = er.asarray([40000.0, 30000.0, 1.0], "float16")
# ExactRef's result for their sum, with the final rounding toward zero
capped = er.sum(x16, rounding=er.Rounding.TOWARD_ZERO)
# The reference and its tolerance
print("rounded:", capped.rounded.tolist())
print("tolerance:", capped.tolerance().tolist())
rounded: 65504.0
tolerance: 0.0
Infinite and NaN inputs follow IEEE rules on the exact sum:
# ExactRef's result for the float32 sum of +inf and -inf
special = er.sum(er.asarray([float("inf"), float("-inf")], "float32"))
# NaN, as IEEE arithmetic gives for inf + (-inf)
print("rounded:", special.rounded.tolist())
# Whether a correct implementation may return NaN, and +inf
print("may be NaN:", special.may_be("nan").tolist())
print("may be +inf:", special.may_be("+inf").tolist())
rounded: nan
may be NaN: True
may be +inf: False
An infinite tolerance means no bound at all: every finite value is in the
range and every flag is set. Single calls on arrays do not produce it.
Chains through exp can. If an earlier step may have stored
+inf, exp can turn that back into a finite value:
# The sum from above, as the one row of a 1x3 array
first = er.sum(er.asarray([[3.4028234e38, 2e31, -1e31]], "float32"), axis=1)
# A correct implementation of it may store +inf
print("first may be +inf:", first.may_be("+inf").tolist())
# The 2x2 identity matrix, as a 1x1x2x2 array
eye = er.asarray([[[[1.0, 0.0], [0.0, 1.0]]]], "float32")
# Keys: eye times first. Passing the result itself, reshaped to 1x1x1x1,
# carries its possible +inf into the keys.
key = er.multiply(eye, first.reshape((1, 1, 1, 1)))
# One query and two values
query = er.asarray([[[[1.0, 2.0]]]], "float32")
value = er.asarray([[[[1.0, 1.0], [2.0, 2.0]]]], "float32")
# ExactRef's result for attention, which applies exp to the scores
a = er.scaled_dot_product_attention(query, key, value)
# The reference and tolerance of attention's first output, named output
print("output:", a.rounded.output.tolist())
print("tolerance:", a.tolerance().output.tolist())
first may be +inf: [True]
output: [[[[2.0, 2.0]]]]
tolerance: [[[[inf, inf]]]]
Several outputs#
Attention returns two outputs, named as in PyTorch: output and
logsumexp. rounded, tolerance() and the other methods hold both, by
name, and r["output"] is a result for just one of them:
# Two queries, three keys and three values, each with 2 features
q = er.asarray([[[1.0, 0.5], [0.25, -1.0]]], "float32")
k = er.asarray([[[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]], "float32")
v = er.asarray([[[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]]], "float32")
# Per query, the keys it keeps (True): keys 0 and 2 for the first query
mask = er.asarray([[[True, False, True], [False, False, False]]], "bool")
# ExactRef's result for attention with that mask
a = er.scaled_dot_product_attention(q, k, v, attn_mask=mask)
# The names of its outputs
print("outputs:", list(a.rounded))
# Each output, read by name
print("output:", a.rounded.output.tolist())
print("logsumexp:", a.rounded.logsumexp.tolist())
print("logsumexp tolerance:", a.tolerance().logsumexp.tolist())
outputs: ['output', 'logsumexp']
output: [[[3.3499159812927246, 4.349915981292725], [0.0, 0.0]]]
logsumexp: [[1.592574954032898, -inf]]
logsumexp tolerance: [[8.344650268554688e-07, nan]]
The mask keeps no key for the second query, so that row's output is 0 and its
logsumexp is -inf (the log of an empty sum), which has a NaN tolerance like
every non-finite reference. logsumexp is in float32 for outputs narrower
than float32, as fused kernels return it.
covers() checks one output at a time, so call it on a["output"] or
a["logsumexp"]. On a itself it raises TypeError.
# An implementation's output: here the reference values printed above
expected = [[[3.3499159812927246, 4.349915981292725], [0.0, 0.0]]]
# a["output"] is a result for that output alone, which covers() can check
print("covered:", a["output"].covers(er.asarray(expected, "float32")).all())
covered: True
Gradients work the same way, with names such as grad_A and grad_query;
Backward checks them against PyTorch's autograd.
Complex outputs#
Each part of a complex element is bounded on its own. tolerance() returns
complex128: its real part bounds the real part of the difference, its
imaginary part the imaginary part. may_be() and covers() add a last axis
of 2, real then imaginary, the layout of torch.view_as_real:
# ExactRef's result for a complex128 sum whose real part is infinite
z = er.sum(er.asarray([complex("inf+1j"), 1 + 1j], "complex128"))
# The exact sum, rounded once
print("rounded:", z.rounded.tolist())
# One bound per part, as one complex number
print("tolerance:", z.tolerance().tolist())
# Per part, real then imaginary: may a correct implementation return +inf?
print("may be +inf:", z.may_be("+inf").tolist())
# Per part: is inf+2.5j correct?
print("inf+2.5j covered:",
z.covers(er.asarray(complex("inf+2.5j"), "complex128")).tolist())
rounded: (inf+2j)
tolerance: (nan+4.440892098500626e-16j)
may be +inf: [True, False]
inf+2.5j covered: [True, False]
The real part is infinite, so it has a NaN tolerance, but the imaginary part is still checked: 2.5 is too far from 2.
Kernels that round more than the definition#
The tolerance of an operation covers implementations that round only where the operation requires. Some kernels round in more places. FlashAttention, for example, rounds the softmax weights to bfloat16 before multiplying them by V. Here is that kernel's arithmetic emulated in PyTorch, on reproducible inputs:
import math
import random
import torch
# Seeded, so every run uses the same values
rng = random.Random(0)
# A float32 tensor of this shape, filled with normally distributed values
def randn(*shape):
values = [rng.gauss(0, 1) for _ in range(math.prod(shape))]
return torch.tensor(values).reshape(shape)
# Queries, keys and values: batch 1, 2 heads, 8 positions, 16 features
q, k, v = randn(1, 2, 8, 16), randn(1, 2, 8, 16), randn(1, 2, 8, 16)
# The scores q times k transposed, scaled by 1/sqrt(16)
scores = (q @ k.transpose(-2, -1)) / 4.0
# True above the diagonal: the later keys each query must not see
causal = torch.ones(8, 8, dtype=torch.bool).triu(1)
# Remove those keys: a score of -inf gets weight 0 in the softmax
scores = scores.masked_fill(causal, float("-inf"))
# The softmax weights, rounded to bfloat16 and back to float32, times v
out = torch.softmax(scores, dim=-1).bfloat16().float() @ v
# ExactRef's result for causal attention on the same inputs
r = er.scaled_dot_product_attention(q, k, v, is_causal=True)
# Does the definition's tolerance cover the kernel's output?
print("covered:", r["output"].covers(out).all())
covered: False
The definition's tolerance does not budget the extra rounding. The matching algorithm does:
# An algorithm that also rounds the weights to bfloat16 before multiplying by v
flash = er.nn.attention.AttentionAlgorithm(weights="bfloat16")
# Does its tolerance cover the kernel's output? covers() takes outputs by name.
print("flash covered:",
flash.covers({"output": out}, q, k, v, is_causal=True).output.all())
# The largest tolerance of any output element: the definition's, then flash's
print("largest tolerance:", max(r.tolerance().output.scalars()))
print("largest flash tolerance:",
max(flash.tolerance(q, k, v, is_causal=True).output.scalars()))
flash covered: True
largest tolerance: 7.62939453125e-06
largest flash tolerance: 0.009089946746826172
Finer than the output format#
r.exact holds a certified enclosure of each exact value: a midpoint and a
radius, with the true value within the radius of the midpoint. Use it when
the implementation you check is more precise than the output format, or to
measure an error exactly.
from fractions import Fraction
# 0.1, 0.2 and 0.3, each rounded to float32
x = er.asarray([0.1, 0.2, 0.3], "float32")
# ExactRef's result for the dot product of x with [1, 3, 7]
d = er.matmul(x, er.asarray([1.0, 3.0, 7.0], "float32"))
# The enclosure of the exact dot product
e = d.exact.value()
# Its midpoint and radius, as exact Fractions
print("midpoint:", repr(e.midpoint.to_fraction()))
print("radius:", repr(e.radius))
# The exact value, rounded once to float32
print("rounded:", d.rounded.tolist())
midpoint: Fraction(375809651, 134217728)
radius: Fraction(0, 1)
rounded: 2.8000001907348633
The radius is 0, so the midpoint is the exact value. A float64 implementation working on these float32 inputs gets it exactly:
# A float64 dot product of x's float32 values, read as Python floats
total = 0.0
for a, b in zip(x.scalars(), [1.0, 3.0, 7.0]):
total += a * b
# Its error: the exact value of total minus the exact dot product
print("error:", repr(Fraction(total) - e.midpoint.to_fraction()))
error: Fraction(0, 1)
Flushing subnormals to zero#
Subnormals are the tiny values below a format's smallest normal number.
Some hardware modes, and torch.set_flush_denormal(True), flush them to
zero. Say so with flush_to_zero=True:
# ExactRef's result for a float32 product of about 1e-40, a float32 subnormal
t = er.multiply(er.asarray([1e-30], "float32"), er.asarray([1e-10], "float32"))
# The reference and its tolerance
print("rounded:", t.rounded.tolist())
print("tolerance:", t.tolerance().tolist())
# The output of an implementation that flushes the product to zero
zero = er.asarray([0.0], "float32")
# Is it correct by default? And with flush_to_zero=True?
print("zero covered:", t.covers(zero).all())
print("zero covered, flushing:", t.covers(zero, flush_to_zero=True).all())
# The tolerance with flushing: 2**-126, the smallest normal float32
print("flushing tolerance:", t.tolerance(flush_to_zero=True).tolist())
rounded: [9.99994610111476e-41]
tolerance: [0.0]
zero covered: False
zero covered, flushing: True
flushing tolerance: [1.1754943508222875e-38]
Matrix exponential#
matrix_exp has no rounding step that every algorithm shares: Taylor and
Pade methods of any degree, with any number of squarings, all compute it.
So its tolerance is a backward error instead: a correct result is
exp(A + dA) for a small dA, rounded to the compute format and then to
the output format.
# ExactRef's result for the matrix exponential of a 2x2 float64 matrix
m = er.matrix_exp(er.asarray([[0.0, 1.0], [-1.0, 0.0]], "float64"))
# The exact exponential, rounded once to float64
print("rounded:")
for row in m.rounded.tolist():
print(row)
# The tolerance of each element, from the backward error
print("tolerance:")
for row in m.tolerance().tolist():
print(row)
# Asking for the normal answer raises an error
try:
m.normal()
except NotImplementedError as error:
print(f"{type(error).__name__}: {error}")
rounded:
[0.5403023058681398, 0.8414709848078965]
[-0.8414709848078965, 0.5403023058681398]
tolerance:
[4.218847493575595e-15, 4.107825191113079e-15]
[4.107825191113079e-15, 4.218847493575595e-15]
NotImplementedError: matrix_exp has no normal answer
There is no single ordinary algorithm to emulate, so there is no normal answer. Guarantees and limits lists the implementations known to exceed this tolerance.