Непонятное поведение numpy.sum при добавлении нулей

Question

Я понимаю, как математически эквивалентные арифметические операции могут приводить к различным результатам из-за численных ошибок (например, суммирования поплавков в разных порядках).

Однако меня удивляет, что добавление нулей в sum может изменить результат. Я думал, что это всегда относится к поплавкам, независимо от того, что: x + 0. == x.

Вот пример. Я ожидал, что все линии будут равны нулю. Может кто-нибудь объяснить, почему это происходит?

M = 4 # number of random values 
Z = 4 # number of additional zeros 
for i in range(20): 
    a = np.random.rand(M) 
    b = np.zeros(M+Z) 
    b[:M] = a 
    print a.sum() - b.sum() 

-4.4408920985e-16 
0.0 
0.0 
0.0 
4.4408920985e-16 
0.0 
-4.4408920985e-16 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
2.22044604925e-16 
0.0 
4.4408920985e-16 
4.4408920985e-16 
0.0 

Это, кажется, не произойдет при меньших значениях M и Z.

Я также позаботился a.dtype==b.dtype.

Вот еще один пример, который также демонстрирует встроенные sum ведет себя питона, как и ожидалось:

a = np.array([0.1,  1.0/3,  1.0/7,  1.0/13, 1.0/23]) 
b = np.array([0.1, 0.0, 1.0/3, 0.0, 1.0/7, 0.0, 1.0/13, 1.0/23]) 
print a.sum() - b.sum() 
=> -1.11022302463e-16 
print sum(a) - sum(b) 
=> 0.0 

Я использую Numpy V1.9.2.

Answer 1

Короткий ответ: Вы видите разницу между

a + b + c + d 

и

(a + b) + (c + d) 

, который из-за точечных неточностей плавающих не то же самое.

Длинный ответ: Numpy реализует парное суммирование как оптимизацию как скорости (это позволяет упростить векторизации), так и ошибки округления.

Реализация суммы numpy может быть найдена here (функция pairwise_sum_@TYPE@). Он по существу выполняет следующие действия:

1. Если длина массива меньше 8, выполняется регулярное суммирование по петле. Вот почему странный результат не наблюдается, если в вашем случае W < 4 - в обоих случаях будет использоваться одно и то же суммирование по петле.
2. Если длина составляет от 8 до 128, она аккумулирует суммы в 8 ячейках r[0]-r[7], а затем суммирует их на ((r[0] + r[1]) + (r[2] + r[3])) + ((r[4] + r[5]) + (r[6] + r[7])).
3. В противном случае он рекурсивно суммирует две половины массива.

Таким образом, в первом случае вы получаете a.sum() = a[0] + a[1] + a[2] + a[3] и во втором случае b.sum() = (a[0] + a[1]) + (a[2] + a[3]), что приводит к a.sum() - b.sum() != 0.