Numpy浮点数舍入误差

Question

在查找一些NumPy的内容时，我看到了一个讨论NumPy.dot()舍入精度的问题：

Numpy：dot(a，b)和(a*b).sum()之间的区别

因为我碰巧有两台搭载Haswell-CPU的不同计算机放在我的桌子上，这应该提供了FMA和其他一切，所以我想我会测试Ophion在第一个答案中给出的示例，我得到了一个令人惊讶的结果：

在更新/安装/修复lapack / blas / atlas / numpy之后，我在两台计算机上都得到了以下结果：

>>> a = np.ones(1000, dtype=np.float128)+1e-14
>>> (a*a).sum()
1000.0000000000199999
>>> np.dot(a,a)
1000.0000000000199948

>>> a = np.ones(1000, dtype=np.float64)+1e-14
>>> (a*a).sum()
1000.0000000000198
>>> np.dot(a,a)
1000.0000000000176

所以标准的乘法+sum()比np.dot()更精确。然而，timeit确认了，在float64和float128上，.dot()版本更快（但不是很快）。

有人可以解释一下吗？

编辑：我意外删除了关于numpy版本的信息：在Python 3.4.0和3.4.1上，1.9.0和1.9.3的结果相同。

- Johannes Hinrichs

---

有趣的是，我只在NumPy 1.9.2上发现了这个差异，而不是NumPy 1.8.2。两者都使用blas + lapack（而不是atlas）。在NumPy 1.8.2上，使用dot和sum产生的结果相同，表明存在相同的舍入事件，在NumPy 1.9.2上，乘法+ sum()更精确。 - Alex Huszagh

请参见http://docs.scipy.org/doc/numpy/release.html#better-numerical-stability-for-sum-in-some-cases。 - Mark Dickinson

此外，https://github.com/numpy/numpy/pull/3685 是实现sum更改的地方。 - Mark Dickinson

有人能确认一下ATLAS是否使用FMA吗？我知道最新的MKL使用了，这应该是一个有趣的比较。 - Daniel

Answer 1

看起来他们最近为了提高数值稳定性，在ndarray.sum中添加了一个特殊的Pairwise Summation。

根据PR 3685，这将影响到：

all add.reduce calls that go over float_add with IS_BINARY_REDUCE true
so this also improves mean/std/var and anything else that uses sum.

请查看这里的代码更改。

- Waylon Flinn

---

一种“分而治之”的算法。是的，那很有道理，谢谢。 - Johannes Hinrichs

这是否意味着我们实际上最好使用替代策略来计算矩阵，而不是集成函数？这将是具有讽刺意味的，因为NumPy的目的是简化科学工作。 - Ando Jurai