Sums of random number reorderings combine to recurring values

Question

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,

for MATLAB 1. try resetting the 35 long state vector of the generator used by rand with the following : rand(‘state’,sum(100*clock)) this resets generator to different state each time. 2. rand('state') returns the 35 long vector that is the state of the uniform generator used by command rand. 3. To change the state of the generator: rand(‘state’,s) resets state to s 4. rand(‘state’,0) resets generator to initial state 5. rand(‘state’,k) resets generator to k-th state — John BG, Aug 14 '23 at 14:47
Check out: https://stackoverflow.com/questions/588004/is-floating-point-math-broken — JonSG, Aug 14 '23 at 14:47
Sometimes it is 0, but most of the time you end up adding two floats along the way that don't perfectly equal a third float. It manifests because we change the order of the elements that are added. — JonSG, Aug 14 '23 at 14:50
@JonSG Yes, the question is why they keep adding to the exact same values despite the numbers being completely different — OverLordGoldDragon, Aug 14 '23 at 14:51
@JohnBG Thanks, though surely the numbers are completely different each time... it's just 100 numbers. The repetition happens only after a few trials, if not one after another (though that's rare) — OverLordGoldDragon, Aug 14 '23 at 14:54
You are running into the default precision of floats. Here is an example output where I also print `sum(g0)` and `sum(g1)` --> `2.220446049250313e-16 1.0386152428024817 1.0386152428024815` take those last two numbers and subtract them `1.0386152428024817 - 1.0386152428024815` now play around with altering the final digit of the second number and see what you get. — JonSG, Aug 14 '23 at 15:01
@JonSG Ah, individual bits - multiples of float epsilon (or something of sort)... Feel free to post as answer. — OverLordGoldDragon, Aug 14 '23 at 15:12

Cris Luengo · Accepted Answer · 2023-08-14T15:58:59.357

2

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.

edited Aug 14 '23 at 15:58

answered Aug 14 '23 at 15:13

Cris Luengo

55,762
10
62
120

1

`1e100` is a valid floating-point number, and `1e-100` is too. But if you add them together, the second one is way below what can be represented at the magnitude of the first, so the result would be equal to the first. Imagine you have 3 decimal digits of precision (with an arbitrary exponent). You can exactly represent 10.0, and you can represent 0.001. Adding them together you get 10.001, but you only have 3 digits of precision, so you must round it to 10.0. `eps` is `eps(1)`, the smallest value you can add to 1 to get to the next representable floating-point number. [continued] – Cris Luengo Aug 14 '23 at 15:35
[continued] Summing 100 `randn` values is expected to produce a relatively small value, maybe 4, maybe 16. Sometimes it'll be larger or smaller of course. Anyway, the next value up from 4 is `eps(4) == 4*eps`, the next value up from 16 is `eps(16) == 16*eps`. So this is the precision of the operations you are performing here. [Note I'm using powers of 2 on purpose here, `eps(5) == eps(4)`, 5 and 4 have the same magnitude in binary.] – Cris Luengo Aug 14 '23 at 15:40
`eps(x)` returns a unique value for sufficiently different `x`, but not for any different `x`. `numel(unique(eps(randn(1, 100)))) == ` 7 through 11, avg 8.65 over 1e6 trials. I guess it's within realm of reason, though still strikes me as unlikely that I'd get repetitions that often. I guess it's likelier since the different roundoff realizations can cancel. Also (2). -- (1) You said "`1e-100` is a valid floating-point", then "cannot be represented as floating-point"? -- (2) There's within-sum rounding, hence my pre-(1) - also see [sum vs np.sum](https://stackoverflow.com/a/71056283/10133797) – OverLordGoldDragon Aug 15 '23 at 13:47
[Comments snapshot](https://i.stack.imgur.com/cWbh1.png) -- removed mine (but yours nicely extend the answer). (Also, @ not needed if only a single user replied to your Q/A) – OverLordGoldDragon Aug 15 '23 at 13:52
1

@OverLordGoldDragon "is a valid floating-point number that cannot be represented exactly" You can represent it, but you'll have a small rounding error. `1e-600` is invalid, it cannot be represented at all, and is rounded to 0. `eps(x)` depends on the (binary) order of magnitude of `x`, not on its exact value. It's the value of the lowest bit in the number. When you increase the exponent by 1, the lowest bit's value is doubled. So all numbers between 1 and 2 have the same `eps`, all numbers between 2 and 4 have double that, all numbers between 4 and 8 have double that, etc. – Cris Luengo Aug 15 '23 at 14:57
1

@OverLordGoldDragon Here's some more reading about floating-point. It's a very complex topic with lots of details and gotchas. [This is **the** reference](https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html), in case you like the hard maths. [This is the simplified version of that paper](https://floating-point-gui.de/) (though I haven't read it, I don't know if it goes into sufficient detail to help you). – Cris Luengo Aug 15 '23 at 15:00
Yep, I've seen that link in like the first week I started actively using Stack Overflow. Somebody gotta make an SO version of that - just the results, not why's and how's. (and looks like your second URL tries that - I'll check it out) Thanks for all the explanations, Cris. – OverLordGoldDragon Aug 15 '23 at 15:10

Sums of random number reorderings combine to recurring values

1 Answers1