Numpy dot在对称乘法方面太聪明了
有人知道关于这种行为的文档吗？
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 - Sym2diff.max()接近机器精度的非零，例如4.4e-16。
这(与0的差异)通常很好...在一个有限精度的世界里，我们不应该感到惊讶。
此外，我猜测numpy在对称产品方面很聪明，以节省flops并确保对称输出……
但我处理的是混沌系统，在调试时，这种微小的差异很快就会变得明显起来。所以我想知道到底发生了什么。
回答 2
Stack Overflow用户
发布于 2017-04-17 23:28:48
此行为是pull request #6932中引入的NumPy 1.11.0更改的结果。从1.11.0的release notes：
以前，所有矩阵产品都使用gemm BLAS操作。现在，如果矩阵乘积位于矩阵及其转置之间，它将使用syrk BLAS操作来提高性能。这种优化已经扩展到@、numpy.dot、numpy.inner和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来说是很便宜的，因为转置的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())除了明显的加速潜力之外，我猜这种特殊外壳的一个优点是，当使用syrk时，您可以保证(我希望，但实际上它将取决于BLAS实现)获得一个完美对称的结果，而不是一个仅仅是对称到数值误差的矩阵。作为对此的测试(诚然不是很好)，我尝试了：
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: FalseStack Overflow用户
发布于 2017-04-17 23:11:22
我怀疑这与将中间浮点寄存器提升到80位精度有关。从某种程度上证实了这个假设，如果我们使用更少的浮点数，我们的结果将始终为0。
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)https://stackoverflow.com/questions/43453707
复制相似问题