BUG: ``np.sum`` gives inconsistent results with ``np.float32`` reshaped arrays

[mlazzarin]

Describe the issue:

If I create an array of dtype np.float32 and of shape (2 ** n, ), then I reshape it, the sum over all elements of the array depends on the shape.

In the code below, I create an array of shape (2 ** n,) that sums to 1.
If I reshape it to (2 ** (n-t), 2 ** t), I sum over the first axis, then over the second one, I obtain 0.5.
If I sum over the second axis, then over the first one, I obtain 1.

More precisely, with n = 25 + k, the np.sum(np.sum(x, axis=0)) gives 2 ** (t-k-1) instead of 1, unless t-k>=1.

Finally, if I replace the array of ones with a random array and manipulate it a bit (commented line in the code below), I can obtain inconsistent results also with smaller values of n e.g. n = 20, even if the difference is less striking than in the previous case.

Reproduce the code example:

import numpy as np

n = 26
t = 1

x = np.ones(2 ** n) / 2 ** n
# x = np.random.random(2 ** n) ** 2 / (2 ** n)
x = x.astype(np.float32)
x = np.reshape(x, (2 ** (n - t), 2 ** t))

print(np.sum(x))
print(np.sum(np.sum(x, axis=0)))
print(np.sum(np.sum(x, axis=1)))

Error message:

No response

NumPy/Python version information:

1.21.4 3.9.7 (default, Sep 16 2021, 13:09:58)
[GCC 7.5.0]

[bashtage]

Possibly these are both within numerical tolerance of a float32.

np.finfo(np.float32).eps*x.shape[0]

4.0

Also:

y.sum(0)

array([0.25, 0.25], dtype=float32)

[mlazzarin]

Thank you for your reply. Now I understand that it may be a precision issue (in fact, the problem disappears if we use np.float64).
I'm still wondering why the following approach (from my code snippet above) has a precision issue:

print(np.sum(np.sum(x, axis=0)))

while the other two give a perfect result:

print(np.sum(x))
print(np.sum(np.sum(x, axis=1)))

Shouldn't they be affected as well?

[seberg]

See the end of the notes of the np.sum documentation: https://numpy.org/doc/stable/reference/generated/numpy.sum.html

For floating point numbers 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 (partial 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 math.fsum function uses a slower but more precise approach to summation. Especially when summing a large number of lower precision floating point numbers, such as float32, numerical errors can become significant. In such cases it can be advisable to use dtype=”float64” to use a higher precision for the output.

[mlazzarin]

Ok, thank you very much, now it's all clear.