2

I calculated the sum over an array and over a zero padded version of the same array:

import numpy as np

np.random.seed(3635250408)

n0, n1 = int(2**16.9), 2**17

xx = np.random.randn(n0)
yy = np.zeros(n1)
yy[:n0] = xx

sx, sy = np.sum(xx), np.sum(yy)

print(f"sx = {sx}, sy = {sy}") # -> sx = -508.33773983674155, sy = -508.3377398367416
print(f"sy - sx:", sy - sx)  # -> sy - sx: -5.68434188608e-14
print("np.ptp(yy[:n0] - xx) =", np.ptp(yy[:n0] - xx)) # -> 0

Why don't I get identical results?

Interestingly, I am able to show similar effects in Mathematica. I am using Python 3.6 (Anaconda 5.0 with MKL support) and Numpy 1.13.3. Perhaps, could it be an MKL issue?

Update: @rich-l and @jkim noted that rounding problems might be the cause. I am not convinced, because adding zero should not alter a floating point number (The problem arose, when investigating a data set of that size - where the deviations were significantly larger).

Dietrich
  • 5,241
  • 3
  • 24
  • 36
  • 3
    Possible duplicate of [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) – Rich L Nov 09 '17 at 12:08
  • The same issue can be seen using small `n` values: `n0, n1 = 22, 24`. (And there are probably even smaller examples.) – unutbu Nov 09 '17 at 12:32
  • @unutbu very interesting - I'm able to reporduce that as well. – Dietrich Nov 09 '17 at 13:42
  • Same case with Python 2. Also, `np.sum(yy[:n0+1])` is equal to `sx`, but `np.sum(yy[:n0+2])` is equal to `sy` while `np.sum(yy[:n0+1]) + yy[n0+1]` is equal to `sx`. So I can't say exactly what's going on, but numpy is probably trying to do some optimizations under the hood that cause it run into a floating point error. – Rob Watts Nov 09 '17 at 23:14

2 Answers2

2

You might be running into floating-point precision issues at this point.

By default, numpy uses double precision floats for storing the values, with 16 digits of precision. The first result outputs 17 digits.

I suspect that in the former case the fluctuations in values result in the two values being rounded slightly differently, with the former being resulting in a rounding to a half (5.5e-16), and the latter exceeding the threshold to be rounded to a full number (6.0e-16).

However, this is just a hypothesis - I don't know for sure how numpy does rounding for the least significant digit.

jkm
  • 704
  • 4
  • 7
2

Floating-point arithmetic is not associative:

In [129]: ((0.1+0.2)+0.3) == (0.1+(0.2+0.3))
Out[129]: False

So the order in which the items are added affects the result. numpy.sum usually uses pairwise summation. It reverts to naive summation (from left to right) when the length of the array is less than 8 or when summing over a strided axis.

Since pairwise summation recursively breaks the sequence into two groups, the addition of zero padding affects the midpoint where the sequence gets divided and hence alters the order in which the values are added. And since floating-point arithmetic is not associative, zero padding can affect the result.

For example, consider 

import numpy as np

np.random.seed(3635250408)
n0, n1 = 6, 8
xx = np.random.randn(n0)
# array([ 1.8545852 , -0.30387171, -0.57164897, -0.40679684, -0.8569989 ,
#         0.32546545])

yy = np.zeros(n1)
yy[:n0] = xx
# array([ 1.8545852 , -0.30387171, -0.57164897, -0.40679684, -0.8569989 ,
#         0.32546545,  0.        ,  0.        ])

xx.sum() and yy.sum() are not the same value:

In [138]: xx.sum()
Out[138]: 0.040734223419930771

In [139]: yy.sum()
Out[139]: 0.040734223419930826

In [148]: xx.sum() == yy.sum()
Out[148]: False

Since len(xx) < 8, the values in xx are summed from left to right:

In [151]: xx.sum() == (((((xx[0]+xx[1])+xx[2])+xx[3])+xx[4])+xx[5])
Out[151]: True

Since len(yy) >= 8, pairwise summation is used to compute yy.sum():

In [147]: yy.sum() == (yy[0]+yy[1]+yy[2]+yy[3])+(yy[4]+yy[5]+yy[6]+yy[7])
Out[147]: True

Related NumPy developer discussions:

numpy.sum does not use Kahan nor Shewchuk summation (used by math.fsum). I believe these algorithms would produce a stable result under the zero-padding issue that you've raised but I'm not expert enough to say for sure.

unutbu
  • 842,883
  • 184
  • 1,785
  • 1,677