In place matrix multiplication numpy

Question

I have a list of point arrays, and as I iterate through it I want to do in-place matrix multiplication i.e. I want the result stored in the same matrix.

The code is essentially:

for p in p_list:
    # R is a 3x3 matrix
    p[:,:] = np.matmul(R,p)

This code shows no error but the result is incorrect, as if the multiplication is performed in the arrays and substituted as calculated, so it creates a wrong output matrix. Removing [:,:] give a correct multiplication.

1) Why does this happen exactly? 2) The main reason I used [:,:] is to make sure the result is stored back in the list p_list. Is there a correct way to do that (without using an intermediate variable)?

Answer 1

matmul accepts an out parameter

If p_list is an ndarray, with shape N, 3, then you can achieve the entire multiplication in one matmul:

np.matmul(p_list, R.T, out=p_list)

Answer 2

With @sacul's example:

In [59]: R.shape
Out[59]: (3, 3)
In [60]: p_list.shape
Out[60]: (2, 3)

In [58]: np.array([np.matmul(R,p) for p in p_list])
Out[58]: 
array([[1.54190819, 2.86411033, 2.13484841],
       [3.7506842 , 6.75463885, 4.89299192]])

einsum generates the same (new) array without external loop:

In [61]: np.einsum('ij,kj->ki',R,p_list)
Out[61]: 
array([[1.54190819, 2.86411033, 2.13484841],
       [3.7506842 , 6.75463885, 4.89299192]])

And like ufunc, it accepts an out parameter:

In [63]: np.einsum('ij,kj->ki',R,p_list, out=p_list)
Out[63]: 
array([[1.54190819, 2.86411033, 2.13484841],
       [3.7506842 , 6.75463885, 4.89299192]])
In [64]: p_list
Out[64]: 
array([[1.54190819, 2.86411033, 2.13484841],
       [3.7506842 , 6.75463885, 4.89299192]])

I'm sure it uses an intermediate buffer, but should be faster than iterating row by row. Using the out is a bit slower than simply letting it return a new array.

With a tweak in dimensions, matmul can perform the whole calc in one call (the key is to pair the last dim of R with the second to the last of p_list (modified)).

In [84]: (R@p_list[:,:,None])[:,:,0]
Out[84]: 
array([[1.54190819, 2.86411033, 2.13484841],
       [3.7506842 , 6.75463885, 4.89299192]])