Iterate over vectors in a multidimensional numpy array

Question

I have a 3xNxM numpy array a, and I would like to iterate over the last two axes: a[:,x,y]. The inelegant approach is:

import numpy as np
a = np.arange(60).reshape((3,4,5))
M = np. array([[1,0,0],
               [0,0,0],
               [0,0,-1]])

for x in arange(a.shape[1]):
    for y in arange(a.shape[2]):
        a[:,x,y] = M.dot(a[:,x,y])

Can this be done with nditer? The goal of this is to perform a matrix multiplication on each entry, e.g. a[:,x,y] = M[:,:,x,y].dot(a[:,x,y]). An alternative, MATLAB-style approach would be to reshape a as (3,N*M) and M as (3,3*N*M) and take a dot product, but that tends to eat up a lot of memory.

Jaime · Accepted Answer · 2013-02-22 16:23:16Z

5

While fooling around with the shapes may make what you are trying to accomplish much more clear, the easiest way of handling this type of problems without thinking too much is with np.einsum:

In [5]: np.einsum('ij, jkl', M, a)
Out[5]: 
array([[[  0,   1,   2,   3,   4],
        [  5,   6,   7,   8,   9],
        [ 10,  11,  12,  13,  14],
        [ 15,  16,  17,  18,  19]],

       [[  0,   0,   0,   0,   0],
        [  0,   0,   0,   0,   0],
        [  0,   0,   0,   0,   0],
        [  0,   0,   0,   0,   0]],

       [[-40, -41, -42, -43, -44],
        [-45, -46, -47, -48, -49],
        [-50, -51, -52, -53, -54],
        [-55, -56, -57, -58, -59]]])

Plus it often comes with a performance bonus:

In [17]: a = np.random.randint(256, size=(3, 1000, 2000))

In [18]: %timeit np.dot(M, a.swapaxes(0,1))
10 loops, best of 3: 116 ms per loop

In [19]: %timeit np.einsum('ij, jkl', M, a)
10 loops, best of 3: 60.7 ms per loop

EDIT einsum is very powerful voodoo. You can also do what the OP asks in the comment below as follows:

>>> a = np.arange(60).reshape((3,4,5))
>>> M = np.array([[1,0,0], [0,0,0], [0,0,-1]])
>>> M = M.reshape((3,3,1,1)).repeat(4,axis=2).repeat(5,axis=3)
>>> np.einsum('ijkl,jkl->ikl', M, b)
array([[[  0,   1,   2,   3,   4],
        [  5,   6,   7,   8,   9],
        [ 10,  11,  12,  13,  14],
        [ 15,  16,  17,  18,  19]],

       [[  0,   0,   0,   0,   0],
        [  0,   0,   0,   0,   0],
        [  0,   0,   0,   0,   0],
        [  0,   0,   0,   0,   0]],

       [[-40, -41, -42, -43, -44],
        [-45, -46, -47, -48, -49],
        [-50, -51, -52, -53, -54],
        [-55, -56, -57, -58, -59]]])

edited Feb 22, 2013 at 16:23

answered Feb 22, 2013 at 2:41

Jaime

67.7k19 gold badges128 silver badges164 bronze badges

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

2 Comments

emarti Over a year ago

Thanks! Is there a way to include element-wise multiplication in einsum? For instance, I have an array of vectors (as above) a of size 3xMxN, and I want to perform a matrix rotation on each vector with a tensor M of size 3x3xMxN. a = arange(60).reshape((3,4,5)) M = array([[1,0,0], [0,0,0], [0,0,-1]]) M = M.reshape((3,3,1,1)).repeat(4,axis=2).repeat(5,axis=3) for x in arange(a.shape[1]): for y in arange(a.shape[2]): a[:,x,y] = M[:,:,x,y].dot(a[:,x,y]) What's the right einsum approach? I've got: b = einsum('ijkl,jxy', M, a) c = b[:,0,0,:,:]

Jaime Over a year ago

@emarti I've edited my answer to cover your question in the comment.

Community · Accepted Answer · 2020-06-20 09:12:55Z

for x in np.arange(a.shape[1]):
    for y in np.arange(a.shape[2]):
        a[:,x,y] = M.dot(a[:,x,y])

is equivalent to

a = np.dot(M,a.swapaxes(0,1))

In [73]: np.dot(M,a.swapaxes(0,1))
Out[73]: 
array([[[  0,   1,   2,   3,   4],
        [  5,   6,   7,   8,   9],
        [ 10,  11,  12,  13,  14],
        [ 15,  16,  17,  18,  19]],

       [[  0,   0,   0,   0,   0],
        [  0,   0,   0,   0,   0],
        [  0,   0,   0,   0,   0],
        [  0,   0,   0,   0,   0]],

       [[-40, -41, -42, -43, -44],
        [-45, -46, -47, -48, -49],
        [-50, -51, -52, -53, -54],
        [-55, -56, -57, -58, -59]]])

Explanation:

For multidimensional arrays, np.dot(M,a) performs a sum product over the last axis of M and the second-to-last axis of a.

a has shape (3,4,5), but we want to sum over the axis with shape 3. Since the second-to-last axis is going to be summed over, we need a.swapaxis(0,1) -- which has shape (4,3,5) -- to move the 3 into the second-to-last axis.

M has shape (3,3), a.swapaxis(0,1) has shape (4,3,5). Removing the last axis of M and the second-to-last axis of a.swapaxis(0,1) leaves you with (3,) and (4,5), so the result returned by np.dot is an array of shape (3,4,5) -- just what we want.

Collectives™ on Stack Overflow

Iterate over vectors in a multidimensional numpy array

2 Answers 2

2 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

2 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related