Python tensor product
<p>I have the following problem. For performance reasons I use <code>numpy.tensordot</code> and have thus my values stored in tensors and vectors. 
One of my calculations look like this:</p>

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

<p><code>&lt;w_j&gt;</code> is the expectancy value of <code>w_j</code> and <code>&lt;sigma_i&gt;</code> the expectancy value of <code>sigma_i</code>. (Perhaps I should now have called is sigma, because it has nothing to do with standart deviation) Now for further calculations I also need the variance. To the get Variance I need to calculate:
<a href="https://i.stack.imgur.com/KLmxY.png" rel="nofollow"><img src="https://i.stack.imgur.com/KLmxY.png" alt="enter image description here"></a></p>

<p>Now when I implemented the first formula into python with <code>numpy.tensordot</code> I was really happy when it worked because this is quite abstract and I am not used to tensors. The code does look like this:</p>

<pre><code>erc = numpy.tensordot(numpy.tensordot(re, ewp, axes=1), ewp, axes=1)
</code></pre>

<p>Now this works and my problem is to write down the correct form for the second formula. One of my attempts was:</p>

<pre><code>serc = numpy.tensordot(numpy.tensordot(numpy.tensordot(numpy.tensordot
(numpy.tensordot(re, re, axes=1), ewp, axes=1), ewp, axes=1)
, ewp, axes=1), ewp, axes=1)
</code></pre>

<p>But this does give me a scalar instead of a vector. Another try was:</p>

<pre><code>serc = numpy.einsum('m, m', numpy.einsum('lm, l -&gt; m',
numpy.einsum('klm, k -&gt; lm', numpy.einsum('jklm, j -&gt; klm',
numpy.einsum('ijk, ilm -&gt; jklm', re, re), ewp), ewp), ewp), ewp)
</code></pre>

<p>The vectors have lenght <code>l</code> and the dimension of the tensor is <code>l * l * l</code>. I hope my problem is understandable and thank you in advance!</p>

<p>Edit: The first formula can in python also written down like: <code>erc2 = numpy.einsum('ik, k -&gt; i', numpy.einsum('ijk, k -&gt; ij', re, ewp), ewp)</code></p>

<p>You could do that with a series of reductions, like so -</p>

<pre><code>p1 = np.tensordot(re,ewp,axes=(1,0))
p2 = np.tensordot(p1,ewp,axes=(1,0))
out = p2**2
</code></pre>

<p><strong>Explanation</strong></p>

<p>First off, we could separate it out into two groups of operations :</p>

<pre><code>Group 1: R(i,j,k) , &lt; wj &gt; , &lt; wk &gt; 
Group 2: R(i,l,m) , &lt; wl &gt; , &lt; wl &gt; 
</code></pre>

<p>The operations performed within these two groups are identical. So, one could compute for one group and derive the final output based off it.</p>

<p>Now, to compute <code>R(i,j,k)</code> , &lt; <code>wj</code> >, &lt; <code>wk</code> and end up with <code>(i)</code> , we need to perform element-wise multiplication along the second and third axes of <code>R</code> with <code>w</code> and then perform <code>sum-reduction</code> along those axes. Here, we are doing it in two steps with two <code>tensordots</code> -</p>

<pre><code>[1] R(i,j,k) , &lt; wj &gt; to get p1(i,k)
[2] p1(i,k) , &lt; wk &gt; to get p2(i)
</code></pre>

<p>Thus, we end up with a vector <code>p2</code>. Similarly with the second group, the result would be an identical vector. So, to get to the final output, we just need to square that vector, i.e. <code>p**2</code>.</p>

python performance numpy numpy-einsum