Efficiently summing outer product for 1D NumPy arrays
<p>I have a function of the form</p>

<p><a href="https://i.stack.imgur.com/MwLCs.gif" rel="nofollow"><img src="https://i.stack.imgur.com/MwLCs.gif" alt="enter image description here"></a></p>

<p>One way to implement this function in numpy is to assemble a matrix to sum over:</p>

<pre><code>y = a*b - np.sum(np.outer(a*b, b), axis=0)
</code></pre>

<p>Is there a better way to implement this function with numpy, one that doesn't involve creating an NxN array?</p>

<p>You could use <a href="https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html" rel="nofollow"><code>np.einsum</code></a> -</p>

<pre><code>y = a*b - np.einsum('i,i,j-&gt;j',a,b,b)
</code></pre>

<p>We can also perform <code>a*b</code> and feed to <code>einsum</code>  -</p>

<pre><code>y = a*b - np.einsum('i,j-&gt;j',a*b,b)
</code></pre>

<p>On the second approach, we can save some runtime by storing <code>a*b</code> and reusing.</p>

<p>Runtime test -</p>

<pre><code>In [253]: a = np.random.rand(4000)

In [254]: b = np.random.rand(4000)

In [255]: %timeit np.sum(np.outer(a*b, b), axis=0)
10 loops, best of 3: 105 ms per loop

In [256]: %timeit np.einsum('i,i,j-&gt;j',a,b,b)
10 loops, best of 3: 24.2 ms per loop

In [257]: %timeit np.einsum('i,j-&gt;j',a*b,b)
10 loops, best of 3: 21.9 ms per loop
</code></pre>

python numpy optimization