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

Question

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

14

假设我们对两个'float32' 2D数组执行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]

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

案例4：检查除[0]以外的索引；与[0]类似，结果在创建后固定数量的数组扩展后开始稳定。输出

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越大，这个阈值似乎就越晚出现，但它总是存在的。

总之：为什么ND-ND数组的np.dot不精确（且不一致）？相关Git

附加信息:

环境: Win-10 操作系统, Python 3.7.4, Spyder 3.3.6 集成开发环境, Anaconda 3.0 2019/10
CPU: i7-7700HQ 2.8 GHz
Numpy v1.16.5

可能的问题库: Numpy MKL - 还有 BLASS 库; 感谢 Bi Rico 的指出

压力测试代码：如上所述，不一致性随着数组大小的增加而加剧；如果上述情况无法重现，则应尝试下面的情况（如果仍无法，请尝试更大的维度）。我的输出

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]))

问题严重程度：所显示的差异“小”，但在对数十亿个数字进行几秒钟乘法运算以及整个运行时间超过万亿时，情况就不再如此；根据这个帖子，报告的模型准确性相差整整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 Desktop -- 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

备注：

与我的系统观察到的情况相比，链接 Colab笔记本和jupyter环境显示出更小的差异（仅适用于前两行）。此外，案例2从未（至今）显示出不精确。
在这个非常有限的样本中，当前（Docker化的）Jupyter环境比IPython环境更容易受到影响。
np.show_config()太长了无法发布，但总结起来：IPython envs基于BLAS/LAPACK；Colab基于OpenBLAS。在IPython Linux envs中，BLAS库是系统安装的--在Jupyter和Colab中，它们来自/opt/conda/lib

更新：被接受的答案是准确的，但是过于笼统和不完整。对于能够解释代码层面行为的人，问题仍然存在 - 即np.dot使用的精确算法，以及它如何解释上述结果中观察到的“一致的不一致性”（还请参见注释）。以下是一些超出我的解密范围的直接实现：sdot.c -- arraytypes.c.src

- OverLordGoldDragon

评论不适合进行长时间的讨论；此对话已被移至聊天室。 - Samuel Liew

ndarrays的通用算法通常忽略数值精度损失。因为为了简单起见，它们沿着每个轴进行“reduce-sum”，所以操作顺序可能不是最优的...请注意，如果您关心精度误差，最好使用float64。 - Vitor SRG

我明天可能没有时间进行审核，所以现在就授予悬赏奖励。 - Paul

@Paul 无论如何，它都将自动授予得票最高的答案 - 但好吧，谢谢你通知我。 - OverLordGoldDragon

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- mmj · Accepted Answer

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

为了更清晰地说明，我们知道，数学上，结果矩阵的每个元素可以被计算为两个向量（长度相等的数字序列）的点积。但这并不是NumPy计算结果矩阵元素的方式。事实上，还有更有效但更复杂的算法，比如Strassen算法，可以在不直接计算行列点积的情况下获得相同的结果。

在使用这种算法时，即使结果矩阵C=A B的元素Cij在数学上被定义为A的第i行与B的第j列的点积，如果你将具有与A相同第i行的矩阵A2与具有与B相同第j列的矩阵B2相乘，则元素C2ij实际上将按照不同的操作顺序计算（取决于整个A2和B2矩阵），可能会导致不同的数值误差。

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