浮点非确定性的原因？包括 NumPy？

Posted 2023-02-16

【中文标题】浮点非确定性的原因？包括 NumPy？

【英文标题】：Causes of floating point non-determinism? Including NumPy?

【发布时间】：2018-12-05 06:55:52

【问题描述】：

IEEE 浮点运算是确定性的，但请参阅 How can floating point calculations be made deterministic? 了解整体浮点计算可能是非确定性的一种方式：

... 并行计算在执行浮点计算的顺序方面是不确定的，这可能导致跨运行的结果不精确。

两部分问题：

整体浮点计算怎么可能是不确定的，产生不完全相等的结果？

考虑一个调用 NumPy、CVXOPT 和 SciPy 子例程（例如 scipy.optimize.fsolve()）的单线程 Python 程序，这些子例程又调用本地库（例如 MINPACK 和 GLPK）以及优化的线性代数子例程（例如 BLAS、ATLAS 和 MKL）。 “If your numpy/scipy is compiled using one of these, then dot() will be computed in parallel (if this is faster) without you doing anything.”

这些原生库是否曾经以引入非确定性结果的方式进行并行化？

假设：

同样的软件，同样的输入，同样的硬件。多次运行的输出应该相等。 如果可行，则非常希望测试执行代码重构后的输出是否相等。 （是的，操作顺序的某些更改可能会使某些输出不相等。） 程序中的所有随机数都是伪随机数，在所有运行中以一致的方式使用相同的种子。

没有未初始化的值。 Python 通常以这种方式是安全的，但 numpy.empty() 返回一个新数组而不初始化条目。目前尚不清楚它在实践中的速度要快得多。所以要小心！

@PaulPanzer 的测试表明numpy.empty() 确实返回了一个未初始化的数组，它可以轻松快速地回收最近的数组：

import numpy as np np.arange(100); np.empty(100, int); np.empty(100, int) np.arange(100, 200.0); np.empty(100, float); np.empty(100, float)

为这些例程获取有用的计时测量是很棘手的！在timeit 循环中，numpy.empty() 可以继续重新分配相同的一两个内存节点。时间与数组大小无关。为了防止回收：

from timeit import timeit timeit('l.append(numpy.empty(100000))', 'import numpy; l = []') timeit('l.append(numpy.zeros(100000))', 'import numpy; l = []')

但是将数组大小减小到 numpy.zeros(10000) 需要 15 倍的时间；将其减少到 numpy.zeros(1000) 需要 1.3 倍的时间（在我的 MBP 上）。令人费解。

另请参阅： Hash values are salted in Python 3 and each dict preserves insertion order。这可能会随着运行的不同而改变操作的顺序。 [我在 Python 2.7.15 中解决这个问题。]

【问题讨论】：

您可能对github.com/numpy/numpy/issues/12394和github.com/numpy/numpy/issues/12465感兴趣

“在我的测量中，empty() 所花费的时间大约与 zeros() [...]”一样长。 这真的取决于。在频繁的dealloc/realloc（例如简单的timeit）上下文中，empty 可以比zeros 快几个数量级。 empty 也有有效的用法，排列的反转只是一个例子。

只是一个简短的轶事。我最近遇到了这个确切的问题，其中一些代码在 C 中进行确定性最小化。在我认为应该只是表面上的改变之后，我得到了稍微不同的结果。事实证明，问题在于我将a *= b/c 之类的行更改为a = a*b/c，这给出了very 略有不同的结果。由于这导致了稍微不同的结果，它导致最小化器尝试不同的下一个点并且效果滚雪球，因此我在两次提交之间得到了明显不同的结果。

@PaulPanzer 原来numpy.empty() 正在清零我 MBP 上的数组。这解释了时间相似性。也许这是 NumPy 调用 Apple 的 Accelerate 框架的结果，这是默认的，而不是 openblas。我会用 MKL 重新测试。

我的猜测是 zeros 调用 calloc 和 calloc 是决定是否回收或要求新内存的人确实知道它是否必须自己进行归零或可以假设它已经完成了。请注意，我不是受过训练的计算机科学家。

【参考方案1】：

我发现我遇到的大多数（不是全部）非确定性问题似乎已在 OpenBLAS 0.3.5 的代码中得到解决。

OpenBLAS 早期版本中的一系列线程问题是fixed in release 0.3.4，但该版本有一个 macOS 兼容性错误，该错误已在 0.3.5 版本的代码中修复。 Apple 的 Accelerate 框架版本 1.1 和 Intel 的 MKL mkl==2019.0 也会出现这些错误。

见how to install OpenBLAS and compile NumPy and SciPy on it。

也许我遇到的其余问题是由于链接到 Accelerate 的其他库引起的？

注意：我仍然愿意接受这个问题的更多答案。