2

I have a tensor phi = np.random.rand(n, n, 3) and a matrix D = np.random.rand(3, 3). I want to multiply the matrix D along the last axis of phi so that the output has shape (n, n, 3). I have tried this

np.einsum("klj,ij->kli", phi, D)

But I am not confident in this notation at all. Basically I want to do

res = np.zeros_like(phi)           
for i in range(n):
    for j in range(n):
        res[i, j, :] = D.dot(phi[i, j, :])
Mad Physicist
  • 107,652
  • 25
  • 181
  • 264
dba
  • 325
  • 1
  • 6
  • 16
  • 1
    Which axes should be multiplied together and which axes are to be reduced by summation? From your desired output shape it looks like you want to multiply the last axis in each array and reduce the first axis in ```D```? The notation for this would be ```ijk,lk->ijk``` (this is almost the same as what you wrote). – Kevin Apr 22 '21 at 23:25
  • 1
    What you wrote in your question is multiplying the last axis in both arrays **and** also reducing the last axis in both arrays. – Kevin Apr 22 '21 at 23:26
  • If phi is an rgb-image, the new image must be a (n,n,3)-matrix where each pixel is transformed by the matrix `D` – dba Apr 22 '21 at 23:29
  • 1
    What is a transformation? Think in terms of multiplication and summing specified axes. Each subscript will denote one axis. By repeating a subscript in both operands you are telling einsum to **multiply** theese two axes. By omitting a subscript in the desired output you are telling einsum to **sum** along this axis. – Kevin Apr 22 '21 at 23:33
  • With transformation I mean, using the matrix multiplication, i.e. `D.dot(phi[i,j,:])`. I basically want to do `D.dot(phi[i,j,:])` on each pixel of `phi`. – dba Apr 22 '21 at 23:37
  • Hmm, do you want to broadcast the last dimension of ```phi``` over ```D```? I think ```ijk,kk->ijk``` will do this - no axes are reduced in this case. – Kevin Apr 22 '21 at 23:42
  • I added the intended operation to my post – dba Apr 22 '21 at 23:47
  • 1
    Thanks, but isn't the intended operation exactly what you are doing with your einsum solution - ```np.einsum("klj,ij->kli",phi,D)``` ? Or are you looking for some explanation on why it works? – Kevin Apr 23 '21 at 00:06
  • I was not sure about that, actually. I learned a lot! thanks! – dba Apr 23 '21 at 00:08
  • 1
    For a small `n` you compare both calculations. – hpaulj Apr 23 '21 at 01:08
  • I think this can be closed. It is correct. – dba Apr 23 '21 at 01:24

1 Answers1

1

You are treating phi as an n, n array of vectors, each of which is to be left-multiplied by D. So you want to keep the n, n portion of the shape exactly as-is. The last (only) dimension of the vectors should be multiplied and summed with the last dimension of the matrix (the vectors are implicitly 3x1):

np.einsum('ijk,lk->ijl', phi, D)

OR

np.einsum('ij,klj->kli', D, phi)

It's likely much simpler to use broadcasting with np.matmul (the @ operator):

np.squeeze(D @ phi[..., None])

You can omit the squeeze if you don't mind the extra unit dimension at the end.

Mad Physicist
  • 107,652
  • 25
  • 181
  • 264