Python/Numpy - Matrix Multiply a 2D Array and Each Row of another 2D Array

Question

What is the best way to do this?

a = 3x3 array
b = 20x3 array
c = 20x3 array = some_dot_function(a, b) where:
c[0] = np.dot(a, b[0])
c[1] = np.dot(a, b[1])
c[2] = np.dot(a, b[2])
...etc...

I know this can be done with a simple python loop or using numpy's apply_along_axis, but I'm wondering if there is any good way to do this entirely within the underlying C code of numpy. I looked at tensordot and some other functions, but didn't have any luck. I also tried the following:

c = np.dot(a, b[:, :, np.newaxis]
#c.shape = (3, 59, 1)

This actually ran and gave results that looked approximately right, except that the resulting array is not 20x3. I may be able to find a way to reshape it into the array I want, but I figured that there must be an easier/cleaner/clearer built-in method that I'm missing?

senderle · Accepted Answer · 2014-01-31 00:54:45Z

7

This gives (what looks to me like) the correct result:

numpy.dot(b, a.T)

Here's some example output:

>>> a = numpy.arange(9).reshape(3, 3)
>>> b = numpy.arange(60).reshape(20, 3)
>>> numpy.dot(b, a.T)
array([[   5,   14,   23],
       [  14,   50,   86],
       [  23,   86,  149],
       [  32,  122,  212],
       ....

edited Jan 31, 2014 at 0:54

answered Jan 31, 2014 at 0:48

senderle

152k36 gold badges218 silver badges244 bronze badges

Sign up to request clarification or add additional context in comments.

Comments

usethedeathstar · Accepted Answer · 2014-02-03 09:28:26Z

import numpy
a = numpy.arange(9).reshape(3,3)
b = numpy.arange(60).reshape(20,3)
c1 = numpy.dot(b, a.T) # as in the answer of senderle
c2 = numpy.einsum('ji,ki->kj',a,b)

and the resulting c1 and c2 are both the same as what you wish (verified with your c[i] = np.dot(a, b[i]) )

the advantage of numpy.einsum is that this trick 'ji,ki->kj' telling what has to be done on what dimension also works for larger dimensions.

more explanation on einsum

for example, if you want to do the following operation:

a = numpy.arange(60.).reshape(3,4,5)
b = numpy.arange(24.).reshape(4,3,2)
d1 = numpy.zeros((5,2))

for i in range(5):
    for j in range(2):
        for k in range(3):
            for n in range(4):
                d1[i,j] += a[k,n,i] * b[n,k,j]

you can do the same thing much faster by doing:

d2 = numpy.einsum('kni,nkj->ij', a, b) 
# the 'kni,nkj->ij' is what you otherwise do with the indices as in 
# d1[i,j] += a[k,n,i] * b[n,k,j]

or if you do not like this way of specifying what has to happen, you can also use numpy.tensordot instead of numpy.einsum, and specify the axes as follows:

d3 = numpy.tensordot(a,b, axes=([1,0],[0,1]))

so this einsum method is very general and can be used to shortcut for-loops (that are otherwise slow, if you do them in python), and are very interesting for funky tensor-stuff

for more info, see http://docs.scipy.org/doc/numpy/reference/generated/numpy.tensordot.html and http://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html

I saw this method, but I can't seem to wrap my head around how it works. It does seem very powerful and I just spend some time reading about it. Even the documentation says the best way to understand it is with examples. Maybe I need to play around with it more, but I don't understand off the bat how the summations are working here.
@ScottB, think of the string as specifying the form that the assignment statement would take if you were writing it in pure python using a series of nested loops. i is the outermost index, j, the next, and so on, and the format is {a indices},{b indices}->{output indices}. So in other words, output[i,j] += a[i,j] * b[j,k] becomes 'ij,jk->ij'. See how that works? IMO, using dot and T is more legible here, but for unusual summations, einsum is absolutely the way to go.
@ScottB i added a few lines after the d2=... to explain a bit more, what you do with those indices, it is quite simple actually once you get your head around it, (just write as in a quadruple for-loop like d1[i,j] = ..., and than write what you do with indices, but that in one line as in d2 = ...

Collectives™ on Stack Overflow

Python/Numpy - Matrix Multiply a 2D Array and Each Row of another 2D Array

2 Answers 2

Comments

3 Comments

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

3 Comments

Related