title: Python: Element-wise multiplication of 3d with 3d arrays
body: <p>I am having some problems implementing the following equation in a performant way using Python:</p>
  <p><a href="https://i.stack.imgur.com/hVsO1.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/hVsO1.gif" alt="Equation" /></a></p>
  <p>beta and gamma are cartesian coordinates {x,y} and b,m are some index value which can be quite big n=10000. I have a working version of the code which is shown below for the simple case of l=2 and m,b = 4 (l and m always have the same length). I checked the code using timeit and the the bottleneck is the element-wise multiplication with an array of size (3,3) and the reshaping of the resulting array into shape (3m,3m).
  Does anybody has an idea how to increase the performance? (I also noticed that my current version suffers a quite big overhead for large values of l....)</p>
  <pre><code>import numpy as np 
  
  g_l3 = np.array([[1, 4, 5],[2, 6, 7]])
  g_l33  = g_l3.reshape(-1, 3, 1) * g_l3.reshape(-1, 1, 3)
  
  A_lm = np.arange(1, 9, 1).reshape(2, 4)
  B_lb = np.arange(7, 15, 1).reshape(2, 4)
  
  AB_lmb = A_lm.reshape(-1, 4, 1) * B_lb.reshape(-1, 1, 4)
  
  D_lmb33 = np.sum(g_l33.reshape(-1, 1, 1, 3, 3) * AB_lmb.reshape(-1, 4, 4, 1, 1), axis=0)
  D = np.concatenate(np.concatenate(D_lmb33, axis=2), axis=0)
  </code></pre>