<div dir="ltr"><div class="gmail_extra"><br><div class="gmail_quote">On 22 May 2015 at 12:15, Mathieu Blondel <span dir="ltr"><<a href="mailto:mathieu@mblondel.org" target="_blank">mathieu@mblondel.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>Right now I am using np.sum(A * B.T, axis=1) for dense data and I have implemented a Cython routine for sparse data.<br></div><div>I haven't benched np.sum(A * B.T, axis=1) vs. np.einsum("ij,ji->i", A, B) yet since I am mostly interested in the sparse case right now.</div></blockquote></div><br><br></div><div class="gmail_extra">In my system, einsum seems to be faster.<br><br><br>In [3]: N = 256<br><br>In [4]: A = np.random.random((N, N))<br><br>In [5]: B = np.random.random((N, N))<br><br>In [6]: %timeit np.sum(A * B.T, axis=1) <br>1000 loops, best of 3: 260 µs per loop<br><br>In [7]: %timeit  np.einsum("ij,ji->i", A, B)<br>10000 loops, best of 3: 147 µs per loop<br><br><br>In [9]: N = 1023<br><br>In [10]: A = np.random.random((N, N))<br><br>In [11]: B = np.random.random((N, N))<br><br>In [12]: %timeit np.sum(A * B.T, axis=1)<br>100 loops, best of 3: 14 ms per loop<br><br>In [13]: %timeit  np.einsum("ij,ji->i", A, B)<br>100 loops, best of 3: 10.7 ms per loop<br><br><br></div><div class="gmail_extra">I have ATLAS installed from the Fedora repos, so not tuned; but einsum is only using one thread anyway, so probably it is not using it (definitely not computing the full dot, because that already takes 200 ms).<br><br></div><div class="gmail_extra">If B is in FORTRAN order, it is much faster (for N=5000).<br><br>In [25]: Bf = B.copy(order='F')<br><br>In [26]: %timeit  np.einsum("ij,ji->i", A, Bf)<br>10 loops, best of 3: 25.7 ms per loop<br><br>In [27]: %timeit  np.einsum("ij,ji->i", A, B)<br>1 loops, best of 3: 404 ms per loop<br><br>In [29]: %timeit np.sum(A * Bf.T, axis=1)<br>10 loops, best of 3: 118 ms per loop<br><br>In [30]: %timeit np.sum(A * B.T, axis=1)<br>1 loops, best of 3: 517 ms per loop<br><br></div><div class="gmail_extra">But the copy is not worth it:<br><br>In [31]: %timeit Bf = B.copy(order='F')<br>1 loops, best of 3: 463 ms per loop<br><br></div><div class="gmail_extra"><br></div><div class="gmail_extra"><br></div><div class="gmail_extra">/David.<br></div><div class="gmail_extra"><br></div></div>