Accept nondeterministic reductions#
The problem#
A GPU reduction that uses atomic additions, a split-K matrix product, or a parallel sum gives a slightly different answer each run, because the order of additions changes with thread scheduling. A test that compares against a fixed expected value fails at random, and the usual fix is to loosen it until it never fails.
Why the usual approach falls short#
Every one of those answers may be correct, or some may come from a real race or a lost update. A loosened tolerance cannot tell them apart, and a fixed expected value is just one of the correct answers.
The solution#
Check each run against the any-order tolerance: it covers every order of additions that does not overflow. Here five shuffled orders of the same 10,000 float32 values stand in for five runs:
import random
import numpy as np
import exactref as er
# Seeded, so every run uses the same values
rng = random.Random(1)
# 10,000 random float32 values: normal, with mean 0 and standard deviation 1
values = np.array([rng.gauss(0, 1) for _ in range(10000)], dtype=np.float32)
# ExactRef's reference for their sum
r = er.sum(values)
# Five runs, each adding the values in its own order
runs = []
for seed in range(5):
# This run's order: the indexes, shuffled with the run's own seed
order = list(range(len(values)))
random.Random(seed).shuffle(order)
# Add one value at a time, rounding to float32 each time
total = np.float32(0)
for i in order:
total = np.float32(total + values[i])
runs.append(total)
# How many different answers the five runs gave
print("different answers:", len(set(runs)))
# Could a correct implementation that adds in any order produce each run's
# answer? er.asarray turns NumPy's result into an ExactRef array.
covered = [bool(r.covers(er.asarray(run), standard="any_order").all())
for run in runs]
print("covered:", covered)
different answers: 4
covered: [True, True, True, True, True]
All four answers are correct. The any-order tolerance here is about 4.7, compared with about 0.013 for normal accumulation:
# The most a correct float32 sum can differ from r.rounded: by normal
# accumulation, and by any order of additions
print("normal tolerance:", r.tolerance().tolist())
print("any-order tolerance:", r.tolerance(standard="any_order").tolist())
normal tolerance: 0.012960433959960938
any-order tolerance: 4.724847793579102
The normal tolerance describes the orders ordinary sequential and SIMD code uses. An arbitrary order can be worse than that and still correct. A lost update (a value added twice or not at all) is usually far outside both.