Yes, summing along the fast axis gives additional numerical precision because numpy uses pairwise summation. But if you sum along the slow axis (in a multi dimensional array), doing the pairwise summation would be much slower, so it doesn't happen.

It is a strange beast, but since you basically get bonus precision, in the end we have it...

There should be an open issue about it somewhere, that it should be documented better (although I have a doubt it would be easy to discover even if documented).

Thanks for the explanation - I didn't know about pairwise summation.

If it helps with the documentation, I went looking for notes relating to this in the docs for numpy.sum.

It looks like #9393 is the issue to document this. Feel free to close this issue. Now that I have a search term, I see that #11331 is also likely very similar to this.

Hmmm, I tried to write a bit:

    The numerical precision of sum (and ``np.add.reduce``) is in general
    limited by directly adding each number individually to the result
    causing rounding errors in every step.
    However, often numpy will use a  numerically better approach
    (pairwise summation) leading to improved precision in many use cases.
    This improved precision is always provided when no ``axis`` is given.
    When ``axis`` is given, it will depend on which axis is summed.
    Technically, to provide the best speed possible, the improved precision
    is only used when the summation is along the fast axis in memory.
    Note that the exact precision may vary depending on other parameters.
    In contrast to NumPy, Python's ``sum`` function uses a slower but more
    precise approach to summation.

but not sure it isn't too technical... Will probably open a PR with it in any case.

Um... what approach does Python's builtin sum() use? Is it a language feature or an implementation detail?

@bbbbbbbbba I have no idea, I guess imlementation detail, see also:

https://en.wikipedia.org/wiki/Kahan_summation_algorithm

The standard library of the Python computer language specifies an fsum function for exactly rounded summation, using the Shewchuk algorithm[9] to track multiple partial sums.

Maybe I should link that site...

Which means its incorrect though, fsum uses Shewchuck, sum does not (which makes sense, it would have to recognize floats).

Actually, isn't it just that builtin sum is starting from 0 by default?

>>> sum(np.ones((10 ** 8,), dtype = np.float32))
100000000.0
>>> type(0 + np.float32(1))
<class 'numpy.float64'>
>>>
>>> sum(np.ones((10 ** 8,), dtype = np.float32), np.float32(0))
16777216.0
>>> type(np.float32(0) + np.float32(1))
<class 'numpy.float32'>
>>>
>>> sum(np.ones((10 ** 8, 1), dtype = np.float32))
array([16777216.], dtype=float32)
>>> type(0 + np.float32([1]))
<class 'numpy.ndarray'>
>>> (0 + np.float32([1])).dtype
dtype('float32')

@bbbbbbbbba I am not sure I parse... Just try to add lots of float64, and you will see differences, especially comparing to fsum.

Oh, I think we are actually on the same page now. Apparently you were thinking of fsum when you talked about the "slower but more precise approach to summation". I noticed that sum sometimes also seemed to be more precise than the naive approach, but it turned out that it was only because it used float64.

Closed by #13737