cross products with einsums
<p>I'm trying to compute the cross-products of many 3x1 vector pairs as fast as possible. This</p>



<pre><code>n = 10000
a = np.random.rand(n, 3)
b = np.random.rand(n, 3)
numpy.cross(a, b)
</code></pre>

<p>gives the correct answer, but motivated by <a href="http://stackoverflow.com/a/20910319/353337">this answer to a similar question</a>, I thought that <code>einsum</code> would get me somewhere. I found that both</p>



<pre><code>eijk = np.zeros((3, 3, 3))
eijk[0, 1, 2] = eijk[1, 2, 0] = eijk[2, 0, 1] = 1
eijk[0, 2, 1] = eijk[2, 1, 0] = eijk[1, 0, 2] = -1

np.einsum('ijk,aj,ak-&gt;ai', eijk, a, b)
np.einsum('iak,ak-&gt;ai', np.einsum('ijk,aj-&gt;iak', eijk, a), b)
</code></pre>

<p>compute the cross product, but their performance is disappointing: Both methods perform much worse than <code>np.cross</code>:</p>

<pre><code>%timeit np.cross(a, b)
1000 loops, best of 3: 628 Âµs per loop
</code></pre>



<pre><code>%timeit np.einsum('ijk,aj,ak-&gt;ai', eijk, a, b)
100 loops, best of 3: 9.02 ms per loop
</code></pre>



<pre><code>%timeit np.einsum('iak,ak-&gt;ai', np.einsum('ijk,aj-&gt;iak', eijk, a), b)
100 loops, best of 3: 10.6 ms per loop
</code></pre>

<p>Any ideas of how to improve the <code>einsum</code>s?</p>

<p>The count of multiply operation of <code>einsum()</code> is more then <code>cross()</code>, and in the newest NumPy version, <code>cross()</code> doesn't create many temporary arrays. So <code>einsum()</code> can't be faster than <code>cross()</code>. </p>

<p>Here is the old code of cross:</p>

<pre><code>x = a[1]*b[2] - a[2]*b[1]
y = a[2]*b[0] - a[0]*b[2]
z = a[0]*b[1] - a[1]*b[0]
</code></pre>

<p>Here is the new code of cross:</p>

<pre><code>multiply(a1, b2, out=cp0)
tmp = array(a2 * b1)
cp0 -= tmp
multiply(a2, b0, out=cp1)
multiply(a0, b2, out=tmp)
cp1 -= tmp
multiply(a0, b1, out=cp2)
multiply(a1, b0, out=tmp)
cp2 -= tmp
</code></pre>

<p>To speedup it, you need cython or numba.</p>

python performance numpy cross-product numpy-einsum