title: operation of Einstein sum of 3D matrices
body: <p>The following code indicates that the Einstein sum of two 3D (2x2x2) matrices is a 4D (2x2x2x2) matrix.</p>
  <pre><code>$ c_{ijlm} = \Sigma_k a_{i,j,k}b_{k,l,m} $
  
  $ c_{0,0,0,0} = \Sigma_k a_{0,0,k}b_{k,0,0} = 1x9 + 5x11 = 64 $ 
  </code></pre>
  <p>But, c_{0,0,0,0} = 35 according to the result below:</p>
  <pre><code>&gt;&gt;&gt; a=np.array([[[1,2],[3,4]],[[5,6],[7,8]]])
  &gt;&gt;&gt; b=np.array([[[9,10],[11,12]],[[13,14],[15,16]]])
  &gt;&gt;&gt; c=np.einsum('ijk,klm-&gt;ijlm', a,b)
  &gt;&gt;&gt; c
  array([[[[ 35,  38],
           [ 41,  44]],
  
          [[ 79,  86],
           [ 93, 100]]],
  
  
         [[[123, 134],
           [145, 156]],
  
          [[167, 182],
           [197, 212]]]])
  </code></pre>
  <p>Could someone explain how the operation is carried out?</p>