Python - a quick way to summarize external products?

Question

I am looking for a quick way to calculate the sum of n external products.

Essentially, I start with two matrices generated from normal distributions - there are n vectors with v elements:

A = np.random.normal(size = (n, v)) B = np.random.normal(size = (n, v)) 

I would like to calculate the external products of each vector of size v in A and B and sum them together.

Note that A * BT does not work - A is of size nxv, while B is of size vx n.

The best I can do is create a loop in which external products are created and then summed up later. I have it like this:

 outers = np.array([A[i] * B[i].T]) 

This creates an nxvxv array (the loop is within the comprehension of the list, which is subsequently converted to an array), which can then be summarized using np.sum(outers, axis = 0) . However, this is pretty slow, and I was wondering if there is a vector function that I could use to speed it up.

If anyone has any advice, I would really appreciate it!

Answer 1

It seems to me that all you need to do is change the order of the transpositions and make AT * B instead of A * BT .

If this is not quite what you need, check out np.einsum , which can make very powerful voodoo. In the above example, you should:

 np.einsum('ij,ik->jk', A, B) 

Answer 2

Also consider np.outer .

 np.array([np.outer(A, B) for i in xrange(n)]).sum(0) 

although np.einsum proposed by @Jamie is a clear winner.

 In [63]: %timeit np.einsum('ij,ik->jk', A, B) 100000 loops, best of 3: 4.61 us per loop In [64]: %timeit np.array([np.outer(A[i], B[i]) for i in xrange(n)]).sum(0) 10000 loops, best of 3: 169 us per loop 

and, of course, their results are identical:

 In [65]: np.testing.assert_allclose(method_outer, method_einsum) 

But, aside, I did not find that AT * B or A * BT successfully transmitted.