为啥 np.dot 不精确？ （n 维数组）

Posted 2023-02-23

【中文标题】为啥 np.dot 不精确？ （n 维数组）

【英文标题】：Why is np.dot imprecise? (n-dim arrays)为什么 np.dot 不精确？ （n 维数组）

【发布时间】：2020-03-03 13:45:25

【问题描述】：

假设我们采用两个'float32' 二维数组中的np.dot：

res = np.dot(a, b)   # see CASE 1
print(list(res[0]))  # list shows more digits

[-0.90448684, -1.1708503, 0.907136, 3.5594249, 1.1374011, -1.3826287]

数字。除了，他们可以改变：

---

案例1：切片a

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(6, 6).astype('float32')

for i in range(1, len(a)):
    print(list(np.dot(a[:i], b)[0])) # full shape: (i, 6)

[-0.9044868,  -1.1708502, 0.90713596, 3.5594249, 1.1374012, -1.3826287]
[-0.90448684, -1.1708503, 0.9071359,  3.5594249, 1.1374011, -1.3826288]
[-0.90448684, -1.1708503, 0.9071359,  3.5594249, 1.1374011, -1.3826288]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]
[-0.90448684, -1.1708503, 0.907136,   3.5594249, 1.1374011, -1.3826287]

结果不同，即使打印的切片来自完全相同的数字相乘。

---

CASE 2

：将a 展平，取b 的一维版本，然后 切片a：

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(1, 6).astype('float32')

for i in range(1, len(a)):
    a_flat = np.expand_dims(a[:i].flatten(), -1) # keep 2D
    print(list(np.dot(a_flat, b)[0])) # full shape: (i*6, 6)

[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]
[-0.3393164, 0.9528787, 1.3627989, 1.5124314, 0.46389243, 1.437775]

---

CASE 3：更强的控制；将所有未涉及的整数设置为 零：将 a[1:] = 0 添加到 CASE 1 代码。结果：差异仍然存在。

---

CASE 4：检查[0]以外的索引；就像[0] 一样，结果从它们的创建点开始稳定固定 # 的数组扩大。 Output

np.random.seed(1)
a = np.random.randn(9, 6).astype('float32')
b = np.random.randn(6, 6).astype('float32')

for j in range(len(a) - 2):
    for i in range(1, len(a)):
        res = np.dot(a[:i], b)
        try:    print(list(res[j]))
        except: pass
    print()

---

因此，对于 2D * 2D 情况，结果不同 - 但对于 1D * 1D 是一致的。从我的一些读数来看，这似乎源于使用简单加法的 1D-1D，而 2D-2D 使用“更高级”的性能提升加法，可能不太精确（例如，成对加法正好相反）。尽管如此，我无法理解为什么一旦a 超过设定的“阈值”，情况 1 中的差异就会消失； a 和b 越大，这个阈值似乎越靠后，但它始终存在。

所有人都说：为什么 np.dot 对于 ND-ND 阵列不精确（且不一致）？ Relevant Git

---

其他信息：

环境：Win-10 OS、Python 3.7.4、Spyder 3.3.6 IDE、Anaconda 3.0 2019/10 CPU：i7-7700HQ 2.8 GHz Numpy v1.16.5

可能的罪魁祸首库：Numpy MKL - 也是 BLASS 库；感谢Bi Rico 的关注

---

压力测试代码：如前所述，较大数组的频率差异会加剧；如果上面不能重现，下面应该是（如果不是，请尝试更大的暗淡）。 My output

np.random.seed(1)
a = (0.01*np.random.randn(9, 9999)).astype('float32') # first multiply then type-cast
b = (0.01*np.random.randn(9999, 6)).astype('float32') # *0.01 to bound mults to < 1

for i in range(1, len(a)):
    print(list(np.dot(a[:i], b)[0]))

---

问题严重性：显示的差异“很小”，但在神经网络上运行时不再如此，在几秒钟内乘以数十亿个数字，在整个运行时达到数万亿个数字；根据this thread，报告的模型精度相差整整百分之十。

下面是向模型提供数组的 gif 图像，该模型基本上是 a[0]，w/ len(a)==1 与 len(a)==32：

---

其他平台结果，感谢Paul的测试：

复制案例 1（部分）：

Google Colab 虚拟机 -- Intel Xeon 2.3 G-Hz -- Jupyter -- Python 3.6.8 Win-10 Pro Docker 桌面 -- Intel i7-8700K -- jupyter/scipy-notebook -- Python 3.7.3 Ubuntu 18.04.2 LTS + Docker -- AMD FX-8150 -- jupyter/scipy-notebook -- Python 3.7.3

注意：这些产生的错误比上面显示的要低得多；第一行的两个条目与其他行中的相应条目的最低有效位相差 1。

案例 1 未复制：

Ubuntu 18.04.3 LTS -- Intel i7-8700K -- IPython 5.5.0 -- Python 2.7.15+ 和 3.6.8（2 次测试） Ubuntu 18.04.3 LTS -- Intel i5-3320M -- IPython 5.5.0 -- Python 2.7.15+ Ubuntu 18.04.2 LTS -- AMD FX-8150 -- IPython 5.5.0 -- Python 2.7.15rc1

注意事项：

linked Colab 笔记本和 jupyter 环境显示的差异（仅在前两行）比在我的系统上观察到的要小得多。此外，案例 2 从未（还）表现出不精确性。 在这个非常有限的示例中，当前（Docker 化的）Jupyter 环境比 IPython 环境更容易受到影响。 np.show_config() 太长无法发布，但总而言之：IPython 环境是基于 BLAS/LAPACK 的； Colab 基于 OpenBLAS。在 IPython Linux 环境中，BLAS 库是系统安装的——在 Jupyter 和 Colab 中，它们来自 /opt/conda/lib

---

更新：接受的答案是准确的，但广泛且不完整。对于任何能够在代码级别解释行为的人来说，这个问题仍然悬而未决——即np.dot 使用的精确算法，以及它如何解释在上述结果中观察到的“一致的不一致”（另见 cmets ）。下面是一些超出我理解的直接实现：sdot.c -- arraytypes.c.src

【问题讨论】：

评论不用于扩展讨论；这个对话是moved to chat。

ndarrays 的通用算法通常忽略数值精度损失。因为为简单起见，他们沿每个轴reduce-sum，操作的顺序可能不是最佳顺序...请注意，如果您介意精度错误，不妨使用float64

我明天可能没有时间复习，所以现在奖励赏金。

@Paul 无论如何，它都会自动授予投票最高的答案 - 但是好的，感谢您的通知

【参考方案1】：

这看起来像是不可避免的数字不精确。正如here 解释的那样，NumPy 使用高度优化、仔细调整的 BLAS 方法进行矩阵乘法。这意味着当矩阵的大小发生变化时，乘以 2 个矩阵的操作序列（求和和乘积）可能会发生变化。

为了更清楚一点，我们知道，数学上，结果矩阵的每个元素都可以计算为两个向量的点积（等长序列数）。但这不是 NumPy 计算结果矩阵的元素的方式。事实上，还有更高效但更复杂的算法，例如Strassen algorithm，无需直接计算行列点积即可获得相同的结果。

当使用这样的算法时，即使结果矩阵 C = AB 的元素 C ij 是数学上定义为 A 的 i-th 行与 B 的 j-th 列的点积，如果您将具有与 A 相同的 i-th 行的矩阵 A2 与具有与 B 相同的 j-th 列，元素 C2 ij 将按照不同的操作顺序实际计算（这取决于整个 A2 和 B2 矩阵），可能会导致不同的数值错误。

这就是为什么，即使在数学上 C ij = C2 ij （就像在您的案例 1 中一样），算法在计算中遵循的不同操作顺序（由于矩阵大小的变化）导致不同的数值误差。数值误差还解释了因环境而略有不同的结果，以及在某些情况下，对于某些环境，数值误差可能不存在。

【讨论】：

感谢您的链接，它似乎包含相关信息 - 但是，您的答案可能更详细，因为到目前为止它是问题下方 cmets 的释义。例如，链接的 SO 直接显示 C 代码，并提供算法级别的解释，因此它正朝着正确的方向前进。

这也不是“不可避免的”，如问题底部所示 - 不精确程度因环境而异，仍然无法解释

@OverLordGoldDragon: (1) 一个简单的加法示例：取数字n，取数字k，使其低于k的最后一个尾数位的精度。对于 Python 的原生浮点数，n = 1.0 和 k = 1e-16 有效。现在让ks = [k] * 100。看到sum([n] + ks) == n，而sum(ks + [n]) > n，即求和顺序很重要。 (2) 现代 CPU 有多个单元可以并行执行浮点 (FP) 运算，并且没有定义在 CPU 上计算 a + b + c + d 的顺序，即使命令 a + b 出现在 c + d 之前也是如此机器码。

@OverLordGoldDragon 您应该知道，您要求程序处理的大多数数字不能用浮点数精确表示。试试format(0.01, '.30f')。如果即使像0.01 这样的简单数字也不能用 NumPy 浮点数精确表示，则无需了解 NumPy 矩阵乘法算法的深层细节即可理解我的回答要点；即不同的起始矩阵导致不同的操作序列，因此数学上相等的结果可能由于数值误差而略有不同。

@OverLordGoldDragon 回复：黑魔法。有一篇论文需要在几个 CS MOOC 上阅读。我的回忆不是那么好，但我认为就是这样：itu.dk/~sestoft/bachelor/IEEE754_article.pdf