Matrix array multiplication in MATLAB

Question

I have a data set consisting of two 1800 x 900 x 3 x 3 arrays, each of which should be interpreted as an 1800x900 array of 3x3 matrices. As part of the analysis, at some point I need to create another such array so that at each point in the 1800x900 array I multiply the corresponding 3x3 matrices together.

There seem to be two ways to do this, which I can think of. The obvious way is

C = zeros(size(A))
for i = 1:900
    for j = 1:1800
        C(i,j,:,:) = A(i,j,:,:)*B(i,j,:,:)
    end
end

But this is a rather long cycle and does not actually use the MATLAB vector. Another way -

C = zeros(size(A))
for i = 1:3
    for j = 1:3
        for k = 1:3
            C(:,:,i,j) = C(:,:,i,j) + A(:,:,i,k).*B(:,:,k,j)
        end
    end
end

where large sizes become vectorized, and I mainly use for loops to implement the Einstein summation convention. It seems really inelegant, though, which makes me think that there should be a better way. There is?

Answer 1

bsxfun permute:

C = permute(sum(bsxfun(@times, A, permute(B, [1 2 5 3 4])), 4), [1 2 3 5 4]);

R2016b bsxfun :

C = permute(sum(A .* permute(B, [1 2 5 3 4]), 4), [1 2 3 5 4]);