title: Write double (triple) sum as inner product?
body: <p>Since my <code>np.dot</code> is accelerated by OpenBlas and Openmpi I am wondering if there was a possibility to write the double sum</p>
  
  <pre><code>for i in range(N):
       for j in range(N):
           B[k,l] += A[i,j,k,l] * X[i,j]
  </code></pre>
  
  <p>as an inner product. Right at the moment I am using</p>
  
  <pre><code>B = np.einsum("ijkl,ij-&gt;kl",A,X)
  </code></pre>
  
  <p>but unfortunately it is quite slow and only uses one processor.
  Any ideas?</p>
  
  <p>Edit:
  I benchmarked the answers given until now with a simple example, seems like they are all in the same order of magnitude:</p>
  
  <pre><code>A = np.random.random([200,200,100,100])
  X = np.random.random([200,200])
  def B1():
      return es("ijkl,ij-&gt;kl",A,X) 
  def B2():
      return np.tensordot(A, X, [[0,1], [0, 1]])
  def B3():
      shp = A.shape
      return np.dot(X.ravel(),A.reshape(shp[0]*shp[1],1)).reshape(shp[2],shp[3])
  
  %timeit B1()
  %timeit B2()
  %timeit B3()
  
  1 loops, best of 3: 300 ms per loop
  10 loops, best of 3: 149 ms per loop
  10 loops, best of 3: 150 ms per loop
  </code></pre>
  
  <p>Concluding from these results I would choose np.einsum, since its syntax is still the most readable and the improvement with the other two are only a factor 2x. I guess the next step would be to externalize the code into C or fortran.</p>