Search code examples
pythonnumpymatlabrandomprecision

Sums of random number reorderings combine to recurring values


g0 = randn(1, 100);
g1 = g0;
g1(2:end) = flip(g1(2:end));
sprintf("%.15e", sum(g0) - sum(g1))
g0 = np.random.randn(100)
g1 = g0.copy()
g1[1:] = g1[1:][::-1]
print(sum(g0) - sum(g1))

In both Python and MATLAB, rerunning these commands enough times will repeat the following values (or their negatives; incomplete list):

8.881784197001252e-15
3.552713678800501e-15
2.6645352591003757e-15
4.440892098500626e-16
1.7763568394002505e-15

In fact, the first and second time I ran them - mat -> py -> mat -> py - they were exactly same, making me think they share the RNG on system level with a delay... (but ignore this happened for the question).

I'll sooner fall through the floor than this coinciding, plus across different languages.

What's happening?


Windows 11, Python 3.11.4, numpy 1.24.4, MATLAB 9.14.0.2286388 (R2023a) Update 3,


Solution

  • Your list of values:

    8.881784197001252e-15
    3.552713678800501e-15
    2.6645352591003757e-15
    4.440892098500626e-16
    1.7763568394002505e-15
    

    matches

    >> [40,16,12,2,8].' * eps
    
    ans =
    
       1.0e-14 *
    
       0.888178419700125
       0.355271367880050
       0.266453525910038
       0.044408920985006
       0.177635683940025
    

    It is totally expected that you'd get rounding errors in this range. Getting the exact same two values in two different systems is not all that big of a coincidence. It happened by chance.

    eps is the machine epsilon, the smallest value you can add to 1 to get to the next number. Results of adding 100 (unit normal) random values mostly happen to be in the range 1-32, with smaller values more likely. We also expect the rounding error to be small w.r.t. the precision of the result. So we should be able to write these numbers as <small integer> * eps(<binary magnitude of result>):

    8.881784197001252e-15  == 5 * eps(8)
    3.552713678800501e-15  == 1 * eps(16)
    2.6645352591003757e-15 == 3 * eps(4)
    4.440892098500626e-16  == 1 * eps(2)
    1.7763568394002505e-15 == 1 * eps(8)
    

    Note also that MATLAB recently changed its implementation of sum, explicitly to reduce rounding errors. And that NumPy uses a similar strategy to sum values in an array.