numpy: broadcast matrix multiply accross array

Question

I have a 3xN array, conceptually an array of N 3-vectors, I want to construct the array which results from matrix multiplying a given 3x3 matrix with each column of the array. Is there a good way to do this in a vectorized manner?

Currently, my problem is 3xN, but I may need to consider 3xNxM (or more) in the future.

Loopy approach

U=numpy.rand( [3,24] )

R=numpy.eye(3) # placeholder

for i in xrange( U.shape[1]):
    U[:,i]=numpy.dot( R, U[:,i] )

@MrE that answer is pretty good, and for my shape (3xN) I avoid the transposes you mentioned there. — Dave
– Dave, Commented Apr 9, 2014 at 17:42

Magellan88 · Accepted Answer · 2014-04-09 21:00:21Z

4

Using np.einsum function you can do it even for the multi dimension problem:

U = np.random.rand(3,24,5) 
R = np.eye(3,3)
result = np.einsum( "ijk,il", U,R )

The notation is a little tricky: the string you give first states the indeces of the dimensions of the arrays; so for U the indeces ijk are the running indeces for each dimension. It follows the einstein summation convention so indeces that have the same letter in the string will be summed over. For details read ScipyDocs. I'm sure in your case the dot is faster, because there is less overhead and it will probably use some blas routine, but as you stated that you want to expand to even more dimensions this might be the way to go.

answered Apr 9, 2014 at 21:00

Magellan88

2,5734 gold badges30 silver badges37 bronze badges

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

1 Comment

Eelco Hoogendoorn Over a year ago

In actual benchmarks, einsum pulls ahead quite often. Especially in cases such as these, involving contractions over small dimensions. I imagine BLAS is invoked by means of the gufunc mechanism; that is, a separate BLAS call on each tiny submatrix. In situations such as these, it is BLAS which will suffer due to an O(N) overhead, which quickly overwhelms the O(1) overhead of einsum.

Community · Accepted Answer · 2017-05-23 12:07:15Z

3

In this case you can simply call np.dot(R, U) and it will work:

import numpy as np

np.random.seed(0)

U = np.random.rand(3,24) 
R = np.random.rand(3,3)

result = np.empty_like(U)

for i in range( U.shape[1]): 
    result[:,i] = np.dot(R, U[:,i])

print result

print np.allclose(result, np.dot(R, U))

For the (3xNxM) case you can reshape to (3x(N.M)), dot and reshape the result back, similar to my answer here

edited May 23, 2017 at 12:07

CommunityBot

11 silver badge

answered Apr 9, 2014 at 17:41

YXD

32.6k15 gold badges79 silver badges117 bronze badges

1 Comment

Magellan88 Over a year ago

I think using np.dot when you can is actually the way to go because it often makes use of really fast libraries such as Blas or Lapack that are higly optimized.

Collectives™ on Stack Overflow

numpy: broadcast matrix multiply accross array

2 Answers 2

1 Comment

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related