关于对称乘法的小圆点太聪明了 [英] Numpy dot too clever about symmetric multiplications

问题描述

有人知道有关此行为的文档吗?

Anybody know about documentation for this behaviour?

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.

diff.max() is near machine-precision non-zero, e.g. 4.4e-16.

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

This (the discrepancy from 0) is usually fine... in a finite-precision world we should not be surprised.

此外，我猜想numpy对于对称产品很聪明，可以节省翻牌并确保对称输出...

Moreover, I would guess that numpy is being smart about symmetric products, to save flops and ensure symmetric output...

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

But I deal with chaotic systems, and this small discrepancy quickly becomes noticeable when debugging. So I'd like to know exactly what's going on.

推荐答案

此行为是发行说明中获得1.11.0:

This behaviour is the result of a change introduced for NumPy 1.11.0, in pull request #6932. From the release notes for 1.11.0:

以前，gem BLAS操作用于所有基质产品. 现在，如果矩阵乘积在矩阵与其转置之间，则它 将使用syrk BLAS操作来提高性能.这 优化已扩展到@，numpy.dot，numpy.inner和 numpy.matmul.

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的更改中，找到了

In the changes for that PR, one finds 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数组只是原始数组上的视图，具有颠倒的形状和跨度，因此只需检查数组中的一些元数据(而不是必须比较实际的数组数据.

So NumPy is making an explicit check for the case of a matrix times its transpose, and calling a different underlying BLAS function in that case. As @hpaulj notes in a comment, such a check is cheap for NumPy, since a transposed 2d array is simply a view on the original array, with inverted shape and strides, so it suffices to check a few pieces of metadata on the arrays (rather than having to compare the actual array data).

这是一个稍微简单的案例，显示出差异.请注意，在dot的参数之一上使用.copy足以击败NumPy的特殊外壳.

Here's a slightly simpler case that shows the discrepancy. Note that using a .copy on one of the arguments to dot is enough to defeat NumPy's special-casing.

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时得到对称结果，而不是仅根据数值误差而对称的矩阵.作为对此(当然不是很好)的测试，我尝试了:

I guess one advantage of this special-casing, beyond the obvious potential for speed-up, is that you're guaranteed (I'd hope, but in practice it'll depend on the BLAS implementation) to get a perfectly symmetric result when syrk is used, rather than a matrix which is merely symmetric up to numerical error. As an (admittedly not very good) test for this, I tried:

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

这篇关于关于对称乘法的小圆点太聪明了的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！