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

Question

假设我们采用两个'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

Answer 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 中一样），算法在计算中遵循的不同操作顺序（由于矩阵大小的变化）导致不同的数值误差。数值误差还解释了因环境而略有不同的结果，以及在某些情况下，对于某些环境，数值误差可能不存在。