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

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

Однако меня удивляет, что добавление нулей в 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.

Вот еще один пример, который также демонстрирует, что встроенный python 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.

Ответ 1

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

a + b + c + d

(a + b) + (c + d)

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

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

Здесь вы можете найти сумму-реализацию numpy здесь (функция [email protected]@). Он по существу делает следующее:

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

Следовательно, в первом случае вы получаете 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.