Fastest way to multiply arrays of matrices in Python (numpy)

I have two arrays of 2-by-2 complex matrices, and I was wondering what would be the fastest method of multiplying them. (I want to do matrix multiplication on the elements of the matrix arrays.) At present, I have

numpy.array(map(lambda i: numpy.dot(m1[i], m2[i]), range开发者_开发百科(l)))

But can one do better than this?

Thanks,

v923z

---

numpy.einsum is the optimal solution for this problem, and it is mentioned way down toward the bottom of DaveP's reference. The code is clean, very easy to understand, and an order of magnitude faster than looping through the array and doing the multiplication one by one. Here is some sample code:

import numpy
l = 100

m1 = rand(l,2,2)
m2 = rand(l,2,2)

m3 = numpy.array(map(lambda i: numpy.dot(m1[i], m2[i]), range(l)))
m3e = numpy.einsum('lij,ljk->lik', m1, m2)

%timeit numpy.array(map(lambda i: numpy.dot(m1[i], m2[i]), range(l)))
%timeit numpy.einsum('lij,ljk->lik', m1, m2)

print np.all(m3==m3e)

Here are the return values when run in an ipython notebook:
1000 loops, best of 3: 479 µs per loop
10000 loops, best of 3: 48.9 µs per loop
True

---

I think the answer you are looking for is here. Unfortunately it is a rather messy solution involving reshaping.

---

If m1 and m2 are 1-dimensional arrays of 2x2 complex matrices, then they essentially have shape (l,2,2). So matrix multiplication on the last two axes is equivalent to summing the product of the last axis of m1 with the second-to-last axis of m2. That's exactly what np.dot does:

np.dot(m1,m2)

Or, since you have complex matrices, perhaps you want to take the complex conjugate of m1 first. In that case, use np.vdot.

PS. If m1 is a list of 2x2 complex matrices, then perhaps see if you can rearrange your code to make m1 an array of shape (l,2,2) from the outset.

If that is not possible, a list comprehension

[np.dot(m1[i],m2[i]) for i in range(l)]

will be faster than using map with lambda, but performing l np.dots is going to be slower than doing one np.dot on two arrays of shape (l,2,2) as suggested above.

---

If m1 and m2 are 1-dimensional arrays of 2x2 complex matrices, then they essentially have shape (l,2,2). So matrix multiplication on the last two axes is equivalent to summing the product of the last axis of m1 with the second-to-last axis of m2. That's exactly what np.dot does:

But that is not what np.dot does.

 a = numpy.array([numpy.diag([1, 2]), numpy.diag([2, 3]), numpy.diag([3, 4])])

produces a (3,2,2) array of 2-by-2 matrices. However, numpy.dot(a,a) creates 6 matrices, and the result's shape is (3, 2, 3, 2). That is not what I need. What I need is an array holding numpy.dot(a[0],a[0]), numpy.dot(a[1],a[1]), numpy.dot(a[2],a[2]) ...

[np.dot(m1[i],m2[i]) for i in range(l)]

should work, but I haven't yet checked, whether it is faster that the mapping of the lambda expression.

Cheers,

v923z

EDIT: the for loop and the map runs at about the same speed. It is the casting to numpy.array that consumes a lot of time, but that would have to be done for both methods, so there is no gain here.

---

May be it is too old question but i was still searching for an answer.

I tried this code

a=np.asarray(range(1048576),dtype='complex');b=np.reshape(a//1024,(1024,1024));b=b+1J*b
%timeit c=np.dot(b,b)
%timeit d=np.einsum('ij, ki -> jk', b,b).T

The results are : for 'dot'

10 loops, best of 3: 174 ms per loop

for 'einsum'

1 loops, best of 3: 4.51 s per loop

I have checked that c and d are same

(c==d).all()
True

still 'dot' is the winner, I am still searching for a better method but no success