有人知道这种行为的文档吗?
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 对对称产品很聪明,以节省失败并确保对称输出......
但我处理的是混沌系统,当调试时,这个小差异很快就会变得明显。所以我想知道到底发生了什么。
最佳答案
此行为是为 NumPy 1.11.0 引入的更改的结果,在 pull request #6932 中.来自 release notes对于 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 来说很便宜,因为转置的 2d 数组只是原始数组的 View ,具有倒置的形状和步幅,因此检查数组上的一些元数据就足够了(而不必比较实际的数组数据)。
这里有一个稍微简单的案例来显示差异。请注意,在 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
关于python - Numpy dot 对对称乘法太聪明了,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/43453707/