title: Combining element-wise and matrix multiplication with multi-dimensional arrays in NumPy
body: <p>I have two multidimensional NumPy arrays, <code>A</code> and <code>B</code>, with <code>A.shape = (K, d, N)</code> and <code>B.shape = (K, N, d)</code>. I would like to perform an element-wise operation over axis 0 (<code>K</code>), with that operation being matrix multiplication over axes 1 and 2 (<code>d, N</code> and <code>N, d</code>). So the result should be a multidimensional array <code>C</code> with <code>C.shape = (K, d, d)</code>, so that <code>C[k] = np.dot(A[k], B[k])</code>. A naive implementation would look like this:</p>
  
  <pre><code>C = np.vstack([np.dot(A[k], B[k])[np.newaxis, :, :] for k in xrange(K)])
  </code></pre>
  
  <p>but this implementation is <em>slow</em>. A slightly faster approach looks like this:</p>
  
  <pre><code>C = np.dot(A, B)[:, :, 0, :]
  </code></pre>
  
  <p>which uses the default behaviour of <code>np.dot</code> on multidimensional arrays, giving me an array with shape <code>(K, d, K, d)</code>. However, this approach computes the required answer <code>K</code> times (each of the entries along axis 2 are the same). Asymptotically it will be slower than the first approach, but the overhead is much less. I am also aware of the following approach:</p>
  
  <pre><code>from numpy.core.umath_tests import matrix_multiply
  C = matrix_multiply(A, B)
  </code></pre>
  
  <p>but I am not guaranteed that this function will be available. My question is thus, does NumPy provide a standard way of doing this efficiently? An answer which applies to multidimensional arrays in general would be perfect, but an answer specific to only this case would be great too.</p>
  
  <p><strong>Edit:</strong> As pointed out by @Juh_, the second approach is incorrect. The correct version is:</p>
  
  <pre><code>C = np.dot(A, B).diagonal(axis1=0, axis2=2).transpose(2, 0, 1)
  </code></pre>
  
  <p>but the overhead added makes it <em>slower</em> than the first approach, even for small matrices. The last approach is winning by a long shot on all my timing tests, for small and large matrices. I'm now strongly considering using this if no better solution crops up, even if that would mean copying the <code>numpy.core.umath_tests</code> library (written in C) into my project.</p>