sum along axis gives different answer for C v F ordered array

[dvincentwest]

when doing a sum along an axis, I get different answers for C versus Fortran ordered arrays.

In the example below I have a large 2D array where all values are zero except for the first row. Slicing the first row and summing gives the same answer for both arrays. However, summing along the second axis and taking the first element of the result gives a different answer for the F ordered array while taking the total sum of the array gives a third answer, while for the C-ordered array, all three methods yield the expected same answer

def test_array_sum_function(arr):
    idx=0
    val1 = arr[idx, :].sum()
    val2 = arr.sum(axis=(1))[idx]
    print('axis sums:', val1)
    print('          ', val2)
    print('    equal:', val1 == val2)
    print('total sum:', arr.sum())

n = 2_000_000
np.random.seed(42)
rnd = np.random.random(n)

print('Fortran order:')
arrF = np.zeros((2, n), order='F')
arrF[0, :] = rnd
test_array_sum_function(arrF)

print('\nC order:')
arrC = np.zeros((2, n), order='C')
arrC[0, :] = rnd
test_array_sum_function(arrC)

prints:

Fortran order:
axis sums: 999813.1414744433
           999813.1414744079
    equal: False
total sum: 999813.1414744424

C order:
axis sums: 999813.1414744433
           999813.1414744433
    equal: True
total sum: 999813.1414744433

Numpy/Python version information:

1.15.4
3.7.1 (default, Dec 14 2018, 19:28:38)
[GCC 7.3.0]

[seberg]

This is expected due to limited numerical precision causing different round-off errors when you do the sum in a different order (which numpy does to optimize speed).

Note that there is one extra caveat here, in that the numerical precision is actually improved over the naive implementation for the "fast" sum (C-order in your case).

[dvincentwest]

Is it expected that this could lead to differences as large as 1-10 percent? I will post an example when I get back to my computer in a few hours

[seberg]

@dvincentwest yes of course. The errors can get arbitrarily large, because large_number + tiny_number == large_number so that at some point things do not change any more, no matter how many small numbers you add. Exactly this type of failure will occur with a lower probability (or at least lower impact) for summing the faster axis (e.g. C-order with axis=-1).

Now, I did not try to estimate the expected errors, but numpy will not make more errors then what you get if you just just do sum += element in C.

I will close this, but please keep discussing. Also if you have a specific suggestion as to what we can improve, we can reopen (I think there probably is a documentation issue open somewhere about the difference between the memory order).

[seberg]

If you are curious, you can check out, e.g.:

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

Numpy does the pairwise (but only chunked) and not always, which is why I said precision can differ.

[dvincentwest]

The code below shows nature of the errors I am seeing. The accumulation of errors is quite systematic and extremely large considering that the number of elements seems to me to be relatively small. Can you comment further on this and what leads to this extremely large discrepancy? Is it really just accumulated round-off errors?

%pylab inline

def random_theta(n, th_lower, th_upper):
    """uniform sampling in spherical coordiantes"""
    costh_lower = cos(th_lower)
    costh_upper = cos(th_upper)
    rng = costh_lower - costh_upper
    costh = np.random.random(n) * rng + costh_lower
    return arccos(costh)

sum_contiguous = list()
sum_discontiguous = list()
n = 5_000_000
for j in range(90):
    dtype = np.float32
#     dtype = np.float64
    arr = np.zeros(shape=(4, n), dtype=dtype, order='F')
    theta = random_theta(n, radians(j), radians(j+1)).astype(dtype)
    arr[0, :] = cos(theta)
    sum_contiguous.append(arr[0, :].sum())
    sum_discontiguous.append(arr.sum(axis=1)[0])
    if j % 30 == 0: print(j)

now visualizing after renaming the lists for the two different dtype options

fig, (ax1, ax2) = subplots(1, 2, figsize=(9, 4), dpi=90, facecolor='white')
ax1.plot(range(90), sum_contiguous_32, label='continguous sum')
ax1.plot(range(90), sum_discontiguous_32, label='discontinguous sum')
ax2.plot(range(90), sum_contiguous_64, label='continguous sum')
ax2.plot(range(90), sum_discontiguous_64, label='discontinguous sum')
ax1.legend()
ax1.set_title('32-bit float')
ax2.legend()
ax2.set_title('64-bit float')
fig.tight_layout()

[seberg]

Is it really just accumulated round-off errors?

Yes, definitely. Note that 32 bits only has about 8 Decimal digits of precision, and 64bit only has about 16 (I think rather less). So 5 Million numbers 5e6 numbers is not small compared to the numerical precision.

Note that the error for your first example is not that huge compared to the precision of floating point.

[dvincentwest]

@seberg thanks very much for the commentary. To convince myself I wrote the below function in numba and saw the same result as the 32-bit discontiguous sum above. I guess you learn something new everyday.

Can you tell me what is the reason for not using the more accurate summation algorithm in a case like the one above? I wonder if a warning from the sum function might be useful to alert the user that they are in danger of significant round-off errors given the precision / number of elements / algorithm being used to sum the numbers. I don't know, just a thought. Thanks again very much

example numba code to confirm the round-off error:

# this code gives the same round-off error seen in the graphs above
from numba import njit
@njit
def nsum(arr):
    rows, cols = arr.shape
    s = np.float32(0)
    for i in range(cols):
        s += arr[0, i]
    return s

[seberg]

Simple: the reason is speed, you should notice that your numba version is around the same speed as the numpy version for the normal sum, but significantly slower for the F-order sum.

EDIT: Well if you make the numba version actually iterate also the first axis (as an outer loop).

[dvincentwest]

@seberg thank you again for taking the time to walk me through all of this, especially since I now realize the underpinnings of my error are quite an old CS problem that I am just now learning. Since you mentioned suggestions, I saw there were other github issues discussing the numpy.sum function and whether or not it was sufficient. I suppose my only suggestion would be to repeat one that I already saw, which is some kind of option to willfully choose a slower routine for the sake of accuracy, either with a specific function call, or a keyword argument to sum. I couldn't find any such option in the numpy documentation. Of course, I have many options open to me right now, summing contiguous arrays, using doubles, but I do have a gut sense that an explicit option for summation algorithms in the function call would be more robust than managing the context of the arrays being summed. At the least, I think that an explicit declaration in the sum documentation that different algorithms are used in different contexts would be a good idea.