python - 为什么 np.dot 不精确？ (n-dim 阵列)

Question

假设我们取 np.dot两个 'float32'二维数组:

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]

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

案例 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]

案例 3 :更强的控制力；将所有非参与元素设置为零:添加 a[1:] = 0到案例 1 代码。结果:差异持续存在。

案例 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 使用“更高级”的、提高性能的加法，这可能不太精确(例如，成对加法正好相反)。尽管如此，我无法理解为什么在案例 1 中差异消失一次 a切过设定的“阈值”；较大的 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，报告的模型准确度相差整整 10 个百分点.

下面是通过向模型馈送而产生的数组的 gif 基本上是 a[0] , 带 len(a)==1对比 len(a)==32 :



其他平台结果，根据并感谢 Paul的测试:

案例1转载(部分) :

Google Colab VM -- 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 -- 英特尔 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 notebook 和 jupyter 环境显示的差异(并且仅针对前两行)比在我的系统上观察到的要小得多。此外，案例 2 从未(还)表现出不精确。

在这个非常有限的样本中，当前(Dockerized)Jupyter 环境比 IPython 环境更容易受到影响。

np.show_config()发布时间太长，但总而言之:IPython 环境基于 BLAS/LAPACK； Colab 基于 OpenBLAS。在 IPython Linux 环境中，BLAS 库是系统安装的——在 Jupyter 和 Colab 中，它们来自/opt/conda/lib

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

Answer 1

这看起来像是不可避免的数字不精确性。如解释 here , NumPy 使用高度优化、精心调整的 BLAS 方法进行矩阵乘法。这意味着当矩阵的大小发生变化时，乘以 2 个矩阵所遵循的操作序列(和和乘积)可能会发生变化。

为了更清楚，我们知道，从数学上讲，结果矩阵的每个元素都可以计算为两个 vector (等长数字序列)的点积。但这是不是 NumPy 如何计算结果矩阵的元素。事实上，还有更高效但更复杂的算法，比如 Strassen algorithm ，在不直接计算行列点积的情况下获得相同的结果。

使用此类算法时，即使元素 C 结果矩阵的 ij C = A B 在数学上定义为 第 i 行的点积第 j 列 乙 ，如果你乘以矩阵 A2 与 具有相同的第 i 行带矩阵 B2 与 具有相同的第 j 列乙 , 元素 C2 ij 将按照不同的操作序列实际计算(这取决于整个 A2 和 B2 矩阵)，可能会导致不同的数值错误。

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