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] )

See my answer here: http://stackoverflow.com/a/22081723/553404 — YXD, Apr 09 '14 at 17:38
@MrE that answer is pretty good, and for my shape (3xN) I avoid the transposes you mentioned there. — Dave, Apr 09 '14 at 17:42

score 4 · Accepted Answer · answered Apr 09 '14 at 21:00

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.

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. — Eelco Hoogendoorn, Apr 09 '14 at 22:00

score 3 · Answer 2 · edited May 23 '17 at 12:07

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 '17 at 12:07

Community

1
1

answered Apr 09 '14 at 17:41

YXD

31,741
15
75
115

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. – Magellan88 Apr 09 '14 at 21:17

numpy: broadcast matrix multiply accross array

2 Answers2

Linked