python - Numpy sum(Numpy 1.15.4，与 MKL 链接)

Question

我正在运行来自 Anaconda 的最新版本的 Python/Numpy(1.15.4)，与 MKL 链接。我使用以下 epsilon:

epsilon = 2**(-53)

因此 1.0 + epsilon 等于 1.0。然后我定义了以下 numpy 数组，除了前 8 个元素等于 1 之外，它填充了 epsilon。

import numpy as np
n = 1000000
a = np.full(n, epsilon)
a[0:8] = 1.0

如果您使用从左到右的经典缩减来计算数组的总和，您应该得到恰好 8.0，因为所有 epsilon 都不会随着缩减的进行而改变。但我很惊讶地看到这一点

print(np.sum(a))

返回 8.000000000111008，这不是我所期望的。我试图将那些放在中间和数组的末尾以检查总和是否没有向后完成，但没有任何变化。知道如何求和吗？

PS:我很清楚浮点运算很棘手，而且 + 与浮点没有关联，这使得归约的结果取决于求和的顺序。但是我找不到这里使用的求和顺序。

Answer 1

如评论中所述，numpy 使用成对求和。因此，当调用堆栈开始解析时，在递归成对求和结束时，求和(大致)看起来像这样:

(1+1) + (1+1) + (1+1) + (1+1) + (epsilon + epsilon) + ... + (epsilon + epsilon)
(2+2) + (2+2) + (2*epsilon) + (2*epsilon) + ... + (2*epsilon)
(4+4) + (4*epsilon) + (4*epsilon) + ... + (4*epsilon)
8 + (8*epsilon) + (8*epsilon) + ... + (8*epsilon)
8 + (16*epsilon) + (16*epsilon) + ... + (16*epsilon)
...
8 + (999992*epsilon)

您说 1.0 + epsilon 等于 1.0 是正确的。很容易认为 x + epsilon == x 对所有 x。当 x 为“大”时它确实成立，但是，当 x == epsilon(即 epsilon + epsilon != epsilon).因此 epsilon + epsilon 术语将开始堆积:

In [27]: epsilon = 2**(-53)

In [28]: 1.0 + epsilon == 1.0
Out[28]: True

In [29]: 2.0 + epsilon == 2.0
Out[29]: True

In [30]: epsilon + epsilon == epsilon
Out[30]: False

In [31]: epsilon
Out[31]: 1.1102230246251565e-16

In [32]: epsilon + epsilon
Out[32]: 2.220446049250313e-16

In [33]: 123*epsilon
Out[33]: 1.3655743202889425e-14

我不能完全得到 numpy 的答案，但我们可以非常接近:

In [36]: 8 + (999992*epsilon)
Out[36]: 8.000000000111022

In [62]: def pairwise_sum(arr):
    ...:   if len(arr) <= 2:
    ...:     return sum(arr)
    ...:   midpoint = len(arr)//2
    ...:   first_half = arr[:midpoint]
    ...:   second_half = arr[midpoint:]
    ...:   return pairwise_sum(first_half) + pairwise_sum(second_half)
    ...:

In [63]: pairwise_sum(a)
Out[63]: 8.0000000001110205

所以这里显然有一些其他差异(一些实现细节)，但希望这能让您相信 numpy 实际上是在使用成对求和。

HTH.