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((接近机器precision non-Zero ，例如。4.4e-16。

这个(从0的差异(通常很好...在有限精确的世界中，我们不应该感到惊讶。

此外，我猜想Numpy对对称产品很聪明，以节省拖鞋并确保对称输出...

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

Answer 1

此行为是引入Numpy 1.11.0的更改的结果，在拉请求＃6932中。从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阵列只是对原始数组的视图，并具有倒置的形状和步伐，因此需要检查阵列上的一些元数据(而不是必须比较实际数组数据(。

这是一个稍微简单的情况，显示了差异。请注意，在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

Answer 2

我怀疑这是关于促进中间浮点寄存器至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)