How to operate with multidimensional arrays (Einstein summation)?

Let's imagine that I have an array of dimension 4 (A_ijkl) that I want to multiply with an array of dimension 2 (B_kl) to obtain an array of dimension 2 (C_ij; i.e., summing over indices k and l). What could be the most efficient (and clean) manner of doing so in Matlab?

And let's go one step further, what if we want to do C_ijxy=A_ijklxy * B_klxy ?

Thank you

Answer 1

[Matt J]

You would reshape() to transform it to a 2D matrix multiplication. In the first case, this would be

   A=rand(3,4,2,2); 
   B=rand(2,2); 

    C=reshape(A,12,4)*B(:); 
   C=reshape(  C   ,3,4);

For the latter case, you would have to permute() and reshape()

    A=rand(3,4,2,2,6,5); B=rand(2,2,7,8);

    Ap=permute(A,[1,2,5,6,3,4]); 
    Apr=reshape(Ap,[],4); Br=reshape(B,4,[]);
    C=reshape(Arp*Br,[3,4,6,5,7,8]);

[LuisGarcia]

Thank you very much for your answer Matt J. I get the approach but I do not see much gain in efficiency when compared to traditional loops (how can that be? in-built functions should be much faster).

Also, I suspect that there is a typo in the last part, as Arp has not been defined (but Apr has). In any case, C is of dimension 6, while it should be of dimension 4. Can you correct those issues so that I can accept the answer?

[Edit: I've also realized that the second example does not fulfil C_ijxy=A_ijklxy * B_klxy; B should be defined as (7,8,2,2)]

[Matt J]

I get the approach but I do not see much gain in efficiency when compared to traditional loops (how can that be? in-built functions should be much faster).

I would have to see how you've coded the loops and the data size you are testing.

In any case, C is of dimension 6, while it should be of dimension 4. Can you correct those issues

I'm no longer sure I understand the definition of C. Should this really be

   C_ijxy= sum_{kl} A_ijklxy * B_klxy

If so, I would download MTIMESX and modify the 2nd example as follows:

     A=rand(3,4,5,6,7,8); B=rand(5,6,7,8);

       Ar=reshape(A,  3*4, 5*6 ,7,8);
       Br=reshape(B,  5*6,  1  ,7,8);

     C=reshape(mtimesx(Ar,Br),[3,4,7,8]);

[LuisGarcia]

Thank you, Matt. You got the definition right, it is just Einstein notation (summation over repeated indices).

However, I aim at doing this without using any external packages. Isn't it possible? In Python, it is extremely straightforward. Thank you!

[Matt J]

You got the definition right, it is just Einstein notation (summation over repeated indices)

Yes, but in your 2nd case the indices x,y are repeated, yet no summation is done over them.

I am at doing this without using any external packages

If you have the Parallel Computing Toolbox (I don't know if you consider this an external package) and a compatible graphics card, then you can do it on the GPU as follows

       A=rand(3,4,5,6,7,8); B=rand(5,6,7,8);

       Ar=gpuArray(reshape(A,  3*4, 5*6 ,7,8));
       Br=gpuArray(reshape(B,  5*6,  1  ,7,8));

     C= gather( reshape(pagefun(@mtimes,Ar,Br),[3,4,7,8]) );

Ideally, you would try to keep your data on the GPU as continuously as possible, and thereby avoid the need to do GPU/CPU transfers every time you want to do the above operation.

There is no other option for the 2nd case that I can see. That's why mtimesx is a popular download.

[LuisGarcia]

This is part of a larger code and using the GPU will be a mess. I'm amazed that there is no way to deal with such operations (other than using loops). But many thanks anyway!

[Matt J]

Sorry, one other way just occurred to me:

     A=rand(3,4,5,6,7,8); B=rand(5,6,7,8);

       Ar=reshape(A,  3*4, 5*6 ,7,8);
       Br=reshape(B,  1,  5*6  ,7,8);

        p=Ar.*Br; %assumes at least R2016b, otherwise use p=bsxfun(@times, Ar,Br)

     C=reshape(sum(p,2),[3,4,7,8]);