Random order for a sum

Parallelization usually means a summation is done with a random order. That may lead to different values if the computation is made many times even though the result should be the same. This example compares summation of random permutation of the same array of values.

Setup

from tqdm import tqdm
import numpy as np
import pandas

unique_values = np.array(
    [2.1102535724639893, 0.5986238718032837, -0.49545818567276], dtype=np.float32
)
random_index = np.random.randint(0, 3, 2000)
assert set(random_index) == {0, 1, 2}
values = unique_values[random_index]

s0 = values.sum()
s1 = np.array(0, dtype=np.float32)
for n in values:
    s1 += n

delta = s1 - s0
print(f"reduced sum={s0}, iterative sum={s1}, delta={delta}")

reduced sum=1555.234130859375, iterative sum=1555.218505859375, delta=-0.015625

There are discrepancies.

Random order

Let’s go further and check the sum of random permutation of the same set. Let’s compare the result with the same sum done with a higher precision (double).

def check_orders(values, n=200, bias=0):
    double_sums = []
    sums = []
    reduced_sums = []
    reduced_dsums = []
    for _i in tqdm(range(n)):
        permuted_values = np.random.permutation(values)
        s = np.array(bias, dtype=np.float32)
        sd = np.array(bias, dtype=np.float64)
        for n in permuted_values:
            s += n
            sd += n
        sums.append(s)
        double_sums.append(sd)
        reduced_sums.append(permuted_values.sum() + bias)
        reduced_dsums.append(permuted_values.astype(np.float64).sum() + bias)

    data = []
    mi, ma = min(sums), max(sums)
    data.append(dict(name="seq_fp32", min=mi, max=ma, bias=bias))
    print(f"min={mi} max={ma} delta={ma-mi}")
    mi, ma = min(double_sums), max(double_sums)
    data.append(dict(name="seq_fp64", min=mi, max=ma, bias=bias))
    print(f"min={mi} max={ma} delta={ma-mi} (double)")
    mi, ma = min(reduced_sums), max(reduced_sums)
    data.append(dict(name="red_f32", min=mi, max=ma, bias=bias))
    print(f"min={mi} max={ma} delta={ma-mi} (reduced)")
    mi, ma = min(reduced_dsums), max(reduced_dsums)
    data.append(dict(name="red_f64", min=mi, max=ma, bias=bias))
    print(f"min={mi} max={ma} delta={ma-mi} (reduced)")
    return data


data1 = check_orders(values)

  0%|          | 0/200 [00:00<?, ?it/s]
  8%|▊         | 16/200 [00:00<00:01, 154.41it/s]
 16%|█▌        | 32/200 [00:00<00:01, 124.73it/s]
 23%|██▎       | 46/200 [00:00<00:01, 128.97it/s]
 31%|███       | 62/200 [00:00<00:00, 138.15it/s]
 38%|███▊      | 77/200 [00:00<00:00, 139.32it/s]
 46%|████▌     | 92/200 [00:00<00:00, 128.20it/s]
 54%|█████▎    | 107/200 [00:00<00:00, 134.00it/s]
 61%|██████    | 122/200 [00:00<00:00, 138.12it/s]
 68%|██████▊   | 136/200 [00:01<00:00, 137.92it/s]
 75%|███████▌  | 150/200 [00:01<00:00, 133.65it/s]
 83%|████████▎ | 166/200 [00:01<00:00, 140.33it/s]
 90%|█████████ | 181/200 [00:01<00:00, 141.66it/s]
 98%|█████████▊| 196/200 [00:01<00:00, 141.59it/s]
100%|██████████| 200/200 [00:01<00:00, 136.69it/s]
min=1555.2161865234375 max=1555.2200927734375 delta=0.00390625
min=1555.2341522574425 max=1555.2341522574425 delta=0.0 (double)
min=1555.2340087890625 max=1555.234130859375 delta=0.0001220703125 (reduced)
min=1555.2341522574425 max=1555.2341522574425 delta=0.0 (reduced)

This example clearly shows the order has an impact. It is usually unavoidable but it could reduced if the sum it close to zero. In that case, the sum would be of the same order of magnitude of the add values.

Removing the average

Computing the average of the values requires to compute the sum. However if we have an estimator of this average, not necessarily the exact value, we would help the summation to keep the same order of magnitude than the values it adds.

mean = unique_values.mean()
values -= mean
data2 = check_orders(values, bias=len(values) * mean)

  0%|          | 0/200 [00:00<?, ?it/s]
  6%|▌         | 12/200 [00:00<00:01, 113.70it/s]
 14%|█▍        | 28/200 [00:00<00:01, 138.98it/s]
 22%|██▏       | 44/200 [00:00<00:01, 144.79it/s]
 30%|██▉       | 59/200 [00:00<00:00, 142.84it/s]
 37%|███▋      | 74/200 [00:00<00:00, 142.51it/s]
 44%|████▍     | 89/200 [00:00<00:00, 143.34it/s]
 52%|█████▏    | 104/200 [00:00<00:00, 142.12it/s]
 60%|██████    | 120/200 [00:00<00:00, 145.70it/s]
 68%|██████▊   | 136/200 [00:00<00:00, 149.34it/s]
 76%|███████▌  | 152/200 [00:01<00:00, 151.56it/s]
 84%|████████▍ | 168/200 [00:01<00:00, 111.97it/s]
 91%|█████████ | 182/200 [00:01<00:00, 118.13it/s]
 98%|█████████▊| 196/200 [00:01<00:00, 106.88it/s]
100%|██████████| 200/200 [00:01<00:00, 127.20it/s]
min=1555.2335205078125 max=1555.2335205078125 delta=0.0
min=1555.234149158001 max=1555.234149158001 delta=0.0 (double)
min=1555.234130859375 max=1555.234130859375 delta=0.0 (reduced)
min=1555.234149158001 max=1555.234149158001 delta=0.0 (reduced)

The differences are clearly lower.

df = pandas.DataFrame(data1 + data2)
df["delta"] = df["max"] - df["min"]
piv = df.pivot(index="name", columns="bias", values="delta")
print(piv)

bias     0.000000    1475.613037
name
red_f32     0.000122         0.0
red_f64          0.0         0.0
seq_fp32    0.003906         0.0
seq_fp64         0.0         0.0

Plots.

ax = piv.plot.barh()
ax.set_title("max(sum) - min(sum) over random orders")
ax.get_figure().tight_layout()
ax.get_figure().savefig("plot_check_random_order.png")