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:

Code
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)
Output
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:

Code
# 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())
Output
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.