Numpy dot对于对称乘法太聪明了

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

此行为是在拉取请求#6932中为NumPy 1.11.0引入的更改的结果.从1.11.0 的发行说明:

以前,gemm BLAS操作用于所有基质产品.现在,如果矩阵乘积在矩阵和它的转置之间,它将使用syrk BLAS操作来提高性能.此优化已扩展为@,numpy.dot,numpy.inner和numpy.matmul.

在该PR的更改中,可以找到以下注释:

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

因此NumPy正在明确检查矩阵的时间及其转置,并在这种情况下调用不同的底层BLAS函数.正如@hpaulj在评论中指出的那样,这样的检查对于NumPy来说是便宜的,因为转置的2d数组只是原始数组上的视图,具有反转的形状和步幅,因此足以检查数组上的几个元数据(而不是必须比较实际的数组数据).

这是一个稍微简单的案例,显示了差异.请注意,使用.copy其中一个参数dot足以击败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