I've run into a case where these two expressions give substantially different results:

a.sum(axis=0)[0]  # 3223672.5
a[:, 0].sum(axis=0)  # 3223679.2

I imagine I'm running into the limits of floating point precision in some form - the data are float32. But I'm puzzled why they're different - I'd expect these two to sum up exactly the same values and so be affected by finite precision in the same way.

Other details:

• The difference is large enough to fail np.testing.assert_allclose
• As far as I can find, all values are 'normal' numbers - no NaN or infinity
• They don't span a particularly large range of values - less than two orders of magnitude.
• All values are positive (I vaguely recall that floating point error can be magnified if large positive and negative numbers cancel each other out).

I can see a difference with random data, but the data I'm using seems to produce a particularly big difference, so I've made an example in a gist including a .npy array with this data.

git clone https://gist.github.com/07c7922f90749ae30d7d45fa02ab6851.git numpy_sum_behaviour
cd numpy_sum_behaviour/
python3 numpy_slice_sum.py

Numpy/Python version information:

>>> import sys, numpy; print(numpy.__version__, sys.version)
1.16.3 3.7.3 (default, May 11 2019, 00:38:04) 
[GCC 9.1.1 20190503 (Red Hat 9.1.1-1)]

I originally saw this on 1.16.4 on another machine, and cut the example down on my laptop with 1.16.3.