Numpy 点对对称乘法太聪明了

Question

有人知道有关此行为的文档吗？

import numpy as np
A  = np.random.uniform(0,1,(10,5))
w  = np.ones(5)
Aw = A*w
Sym1 = Aw.dot(Aw.T)
Sym2 = (A*w).dot((A*w).T)
diff = Sym1 - Sym2

diff.max() 接近机器精度非零，例如4.4e-16.

这（与 0 的差异）通常很好...在有限精度的世界中，我们不应该感到惊讶。

此外，我猜想 numpy 在对称产品方面很聪明，可以节省失败并确保对称输出...

但我处理的是混乱的系统，当调试时，这个小差异很快就会变得明显。所以我想知道到底发生了什么。

Answer 1

我怀疑这与将中间浮点寄存器提升到 80 位精度有关。在某种程度上证实了这个假设的是，如果我们使用更少的浮点数，我们的结果始终为 0，ala

A  = np.random.uniform(0,1,(4,2))
w  = np.ones(2)
Aw = A*w
Sym1 = Aw.dot(Aw.T)
Sym2 = (A*w).dot((A*w).T)
diff = Sym1 - Sym2
# diff is all 0's (ymmv)

Answer 2

此行为是为 NumPy 1.11.0 引入的更改的结果，在 pull request #6932. From the release notes for 1.11.0:

Previously, gemm BLAS operations were used for all matrix products. Now, if the matrix product is between a matrix and its transpose, it will use syrk BLAS operations for a performance boost. This optimization has been extended to @, numpy.dot, numpy.inner, and numpy.matmul.

在该 PR 的更改中，可以找到 this comment:

/*
 * Use syrk if we have a case of a matrix times its transpose.
 * Otherwise, use gemm for all other cases.
 */

所以 NumPy 正在对矩阵乘以其转置的情况进行显式检查，并在这种情况下调用不同的底层 BLAS 函数。正如 @hpaulj 在评论中指出的那样，这种检查对于 NumPy 来说很便宜，因为转置二维数组只是原始数组的一个视图，具有倒置的形状和步幅，所以检查数组上的一些元数据就足够了（而不是必须比较实际的数组数据）。

这里有一个稍微简单的例子来说明差异。请注意，在 dot 的参数之一上使用 .copy 足以击败 NumPy 的特殊外壳。

import numpy as np
random = np.random.RandomState(12345)
A = random.uniform(size=(10, 5))
Sym1 = A.dot(A.T)
Sym2 = A.dot(A.T.copy())
print(abs(Sym1 - Sym2).max())

我想除了明显的加速潜力之外，这种特殊外壳的一个优点是你可以保证（我希望，但实际上这将取决于 BLAS 实现）得到使用 syrk 时的完美对称结果，而不是仅对称到数值误差的矩阵。作为对此的（诚然不是很好）测试，我尝试了：

import numpy as np
random = np.random.RandomState(12345)
A = random.uniform(size=(100, 50))
Sym1 = A.dot(A.T)
Sym2 = A.dot(A.T.copy())
print("Sym1 symmetric: ", (Sym1 == Sym1.T).all())
print("Sym2 symmetric: ", (Sym2 == Sym2.T).all())

我机器上的结果：

Sym1 symmetric:  True
Sym2 symmetric:  False

Numpy 点对对称乘法太聪明了

Numpy dot too clever about symmetric multiplications

python

floating-point

numpy

linear-algebra

floating-accuracy