Compare the accuracy of two implementations#

The problem#

You rewrote a kernel and want to know whether the new version is more or less accurate than the old one, or you are choosing between two summation methods. Comparing both against a float64 result measures their distance from another rounded number, not from the true answer.

Why the usual approach falls short#

Both implementations may be correct, so a pass or fail check says nothing about which is better. A float64 reference has its own rounding error, which blurs small differences, and an error in absolute terms is hard to compare across magnitudes.

The solution#

Measure each implementation's error exactly, against r.exact, in units in the last place (ulps) of the output format. Here a plain sequential float32 loop and NumPy's np.sum, which adds in a tree, both add the same 100,000 values:

Code
import math
import random
from fractions import Fraction
import numpy as np
import exactref as er

# Seeded, so every run uses the same values
rng = random.Random(0)
# 100,000 random float32 values between 0 and 1
values = np.array([rng.uniform(0, 1) for _ in range(100000)], dtype=np.float32)

# Implementation 1: add one value at a time, rounding to float32 each time
sequential = np.float32(0)
for v in values:
    sequential = np.float32(sequential + v)

# Implementation 2: NumPy's sum, which adds in a tree
pairwise = np.sum(values)

# ExactRef's reference for the same sum
r = er.sum(values)
# The exact sum as a Fraction. A sum's enclosure has radius 0, so its
# midpoint is the exact value.
exact = r.exact.value().midpoint.to_fraction()

# How far an output is from the exact sum, in float32 ulps at the answer
def error_in_ulps(output):
    # frexp writes the answer as m * 2**exponent, with 0.5 <= m < 1
    exponent = math.frexp(r.rounded.tolist())[1]
    # float32 keeps 24 significant bits, so one ulp there is 2**(exponent - 24)
    ulp = Fraction(2) ** (exponent - 24)
    # Fraction(float(output)) is the output's exact value, so no rounding here
    return float(abs(Fraction(float(output)) - exact) / ulp)

# Each implementation's error, to 2 decimals
print("sequential:", round(error_in_ulps(sequential), 2))
print("pairwise:", round(error_in_ulps(pairwise), 2))
Output
sequential: 33.36
pairwise: 0.36

The tree sum is almost as accurate as float32 allows; the loop is 33 ulps off. Both are correct, and the tolerance shows how much of its budget each uses:

Code
# The most a correct float32 sum of these values can differ from r.rounded
print("tolerance:", r.tolerance().tolist())

# Could a correct implementation produce each output? covers() answers per
# element, and er.asarray turns NumPy's result into an ExactRef array.
print("sequential covered:", r.covers(er.asarray(sequential)).all())
print("pairwise covered:", r.covers(er.asarray(pairwise)).all())
Output
tolerance: 149.6640625
sequential covered: True
pairwise covered: True

To compare methods before writing any code, the algorithms' error models give the worst case from the shape and format alone:

Code
# The module that holds the algorithms. A sum is a dot product with a vector
# of ones, so take B = 1 in their bounds.
products = er.linalg.products

# Products added one at a time in float32, each rounded before it is added.
# error_model(k) bounds the error of one output element built from k products.
sequential_algorithm = products.matmul_sequential("float32", fma=False)
print("sequential:", sequential_algorithm.error_model(100000))

# The same products added in a balanced tree
pairwise_algorithm = products.matmul_pairwise("float32", fma=False)
print("pairwise:", pairwise_algorithm.error_model(100000))
Output
sequential: 0.006 * sum over k of |A[i,k] * B[k,j]| + 1.41e-40
pairwise: 1.13e-06 * sum over k of |A[i,k] * B[k,j]| + 2.1e-40