圆点对对称乘法太聪明

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的发行说明中：

以前，gem 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在其中一个参数上使用adot足以击败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