blocks|key|2855051|text|我怀疑这与将中间浮点寄存器提升到80位精度有关。从某种程度上证实了这个假设，如果我们使用更少的浮点数，我们的结果将始终为0。|type|unstyled|depth|inlineStyleRanges|entityRanges|data|2855052|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)|code-block|syntax|javascript|2855053|entityMap^0|0|0^^$0|@$1|2|3|4|5|6|7|I|8|@]|9|@]|A|$]]|$1|B|3|C|5|D|7|J|8|@]|9|@]|A|$E|F]]|$1|G|3|-4|5|6|7|K|8|@]|9|@]|A|$]]]|H|$]]

I suspect this is to do with promotion of intermediate floating point registers to 80 bit precision. Somewhat confirming this hypothesis is that if we use fewer floats we consistently get 0 in our results, ala

<pre><code>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)
</code></pre>

blocks|key|2643523|text|此行为是pull+request+#6932中引入的NumPy+1.11.0更改的结果。从1.11.0的release+notes：|type|unstyled|depth|inlineStyleRanges|entityRanges|offset|length|data|2643524|2643525|以前，所有矩阵产品都使用gemm+BLAS操作。现在，如果矩阵乘积位于矩阵及其转置之间，它将使用syrk+BLAS操作来提高性能。这种优化已经扩展到@、numpy.dot、numpy.inner和numpy.matmul。|blockquote|2643526|2643527|2643528|在该PR的更改中，可以找到this+comment|2643529|/*
+*+Use+syrk+if+we+have+a+case+of+a+matrix+times+its+transpose.
+*+Otherwise,+use+gemm+for+all+other+cases.
+*/|code-block|syntax|javascript|2643530|因此，NumPy对矩阵乘以转置的情况进行显式检查，并在这种情况下调用不同的底层BLAS函数。正如@hpaulj在评论中指出的那样，这样的检查对于NumPy来说是很便宜的，因为转置的2d数组只是原始数组上的一个视图，形状和步长都颠倒了，所以它足以检查数组上的一些元数据(而不是必须比较实际的数组数据)。|2643531|这里有一个稍微简单一点的案例，它显示了这种差异。请注意，在dot的一个参数上使用.copy就足以击败NumPy的特殊大小写。|style|CODE|2643532|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())|2643533|除了明显的加速潜力之外，我猜这种特殊外壳的一个优点是，当使用syrk时，您可以保证(我希望，但实际上它将取决于BLAS实现)获得一个完美对称的结果，而不是一个仅仅是对称到数值误差的矩阵。作为对此的测试(诚然不是很好)，我尝试了：|2643534|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())|2643535|我机器上的结果：|2643536|Sym1+symmetric:++True
Sym2+symmetric:++False|2643537|entityMap|0|LINK|mutability|MUTABLE|url|https://github.com/numpy/numpy/pull/6932|1|https://docs.scipy.org/doc/numpy/release.html#id11|2|https://github.com/numpy/numpy/blob/v1.12.1/numpy/core/src/multiarray/cblasfuncs.c#L715-L718^0|4|I|0|1G|D|1|0|0|0|0|0|D|C|2|0|0|0|T|3|14|5|0|0|U|4|0|0|0|0^^$0|@$1|2|3|4|5|6|7|1I|8|@]|9|@$A|1J|B|1K|1|1L]|$A|1M|B|1N|1|1O]]|C|$]]|$1|D|3|-4|5|6|7|1P|8|@]|9|@]|C|$]]|$1|E|3|F|5|G|7|1Q|8|@]|9|@]|C|$]]|$1|H|3|-4|5|6|7|1R|8|@]|9|@]|C|$]]|$1|I|3|-4|5|6|7|1S|8|@]|9|@]|C|$]]|$1|J|3|K|5|6|7|1T|8|@]|9|@$A|1U|B|1V|1|1W]]|C|$]]|$1|L|3|M|5|N|7|1X|8|@]|9|@]|C|$O|P]]|$1|Q|3|R|5|6|7|1Y|8|@]|9|@]|C|$]]|$1|S|3|T|5|6|7|1Z|8|@$A|20|B|21|U|V]|$A|22|B|23|U|V]]|9|@]|C|$]]|$1|W|3|X|5|N|7|24|8|@]|9|@]|C|$O|P]]|$1|Y|3|Z|5|6|7|25|8|@$A|26|B|27|U|V]]|9|@]|C|$]]|$1|10|3|11|5|N|7|28|8|@]|9|@]|C|$O|P]]|$1|12|3|13|5|6|7|29|8|@]|9|@]|C|$]]|$1|14|3|15|5|N|7|2A|8|@]|9|@]|C|$O|P]]|$1|16|3|-4|5|6|7|2B|8|@]|9|@]|C|$]]]|17|$18|$5|19|1A|1B|C|$1C|1D]]|1E|$5|19|1A|1B|C|$1C|1F]]|1G|$5|19|1A|1B|C|$1C|1H]]]]

This behaviour is the result of a change introduced for NumPy 1.11.0, in <a href="https://github.com/numpy/numpy/pull/6932" rel="noreferrer">pull request #6932</a>. From the <a href="https://docs.scipy.org/doc/numpy/release.html#id11" rel="noreferrer">release notes</a> for 1.11.0:

<blockquote>
 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.
</blockquote>

In the changes for that PR, one finds <a href="https://github.com/numpy/numpy/blob/v1.12.1/numpy/core/src/multiarray/cblasfuncs.c#L715-L718" rel="noreferrer">this comment</a>:

<pre><code>/*
 * Use syrk if we have a case of a matrix times its transpose.
 * Otherwise, use gemm for all other cases.
 */
</code></pre>

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).

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

<pre><code>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())
</code></pre>

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 <code>syrk</code> 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:

<pre><code>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())
</code></pre>

Results on my machine:

<pre><code>Sym1 symmetric: True
Sym2 symmetric: False
</code></pre>

Anybody know about documentation for this behaviour?

<pre><code>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
</code></pre>

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

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

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.

Numpy dot too clever about symmetric multiplications

有人知道关于这种行为的文档吗？import numpy as npA  = np.random.uniform(0,1,(10,5))w  = np.ones(5)Aw = A*wSym1 = Aw.dot(Aw.T)Sym2 = (A*w).dot((A*w).T)diff = Sym1 - Sym2diff.max...

问Numpy dot在对称乘法方面太聪明了
EN

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

问Numpy dot在对称乘法方面太聪明了EN

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

问Numpy dot在对称乘法方面太聪明了
EN