title: Numpy: How to perform the outer product of four vectors and then summing them?
body: <p><code>A = [a1, a2]</code>, <code>B = [b1, b2]</code>, <code>C = [c1, c2]</code>, <code>D = [d1, d2]</code>.</p>
  
  <p>The outer product of <code>A, B, C, D</code> is just:</p>
  
  <pre><code>[
      [
          [[a1*b1*c1*d1, a1*b1*c1*d2], [a1*b1*c1*d1, a1*b1*c2*d2]], 
          [[a1*b2*c2*d1, a1*b2*c1*d2], [a1*b2*c2*d1, a1*b2*c2*d2]]
      ],
      [
          [[a2*b1*c1*d1, a2*b1*c1*d2], [a2*b1*c2*d1, a2*b1*c2*d2]], 
          [[a2*b2*c1*d1, a2*b2*c1*d2], [a2*b2*c2*d1, a2*b2*c2*d2]]
      ]
  ]
  </code></pre>
  
  <p>In numpy, the outer product is thus an array of shape <code>(2, 2, 2, 2)</code>. I know how to do this. The problem is, I have 100 this kind of A, B, C, D. They are taken from A_k, B_k, C_k, D_k by fixing k. So they are actually array of shape (2, 100) or (100, 2). Actually may be (20, 100) or (100, 20). I use 2 to simplify the writing. What I want to do is: do this kind of outer product for 100 times and then sum them up. That is: to sum 100 arrays of shape (2, 2, 2, 2). How can I do this efficiently with numpy? If better, is there any GPU acceleration way of doing this?</p>
  
  <p>If A, B, C, D are just 1d array, I know how to do this. The following one-line code works fine:</p>
  
  <pre><code>a[:, None, None, None] * b[None, :, None, None] * c[None, None, :, None] * d[None, None, None, :]
  </code></pre>