MATLAB: How to vector-multiply two arrays of matrices?

Question

I have two 3-dimensional arrays, the first two dimensions of which represent matrices and the last one counts through a parameterspace, as a simple example take

A = repmat([1,2; 3,4], [1 1 4]);

(but assume A(:,:,j) is different for each j). How can one easily perform a per-j matrix multiplication of two such matrix-arrays A and B?

C = A; % pre-allocate, nan(size(A,1), size(B,2)) would be better but slower
for jj = 1:size(A, 3)
  C(:,:,jj) = A(:,:,jj) * B(:,:,jj);
end

certainly does the job, but if the third dimension is more like 1e3 elements this is very slow since it doesn't use MATLAB's vectorization. So, is there a faster way?

Answer 1

I did some timing tests now, the fastest way for 2x2xN turns out to be calculating the matrix elements:

C = A;
C(1,1,:) = A(1,1,:).*B(1,1,:) + A(1,2,:).*B(2,1,:);
C(1,2,:) = A(1,1,:).*B(1,2,:) + A(1,2,:).*B(2,2,:);
C(2,1,:) = A(2,1,:).*B(1,1,:) + A(2,2,:).*B(2,1,:);
C(2,2,:) = A(2,1,:).*B(1,2,:) + A(2,2,:).*B(2,2,:);

---

In the general case it turns out the for loop is actually the fastest (don't forget to pre-allocate C though!).

Should one already have the result as cell-array of matrices though, using cellfun is the fastest choice, it is also faster than looping over the cell elements:

C = cellfun(@mtimes, A, B, 'UniformOutput', false);

However, having to call num2cell first (Ac = num2cell(A, [1 2])) and cell2mat for the 3d-array case wastes too much time.

---

Here's some timing I did for a random set of 2 x 2 x 1e4:

 array-for: 0.057112
 arrayfun : 0.14206
 num2cell : 0.079468
 cell-for : 0.033173
 cellfun  : 0.025223
 cell2mat : 0.010213
 explicit : 0.0021338

Explicit refers to using direct calculation of the 2 x 2 matrix elements, see bellow. The result is similar for new random arrays, cellfun is the fastest if no num2cell is required before and there is no restriction to 2x2xN. For general 3d-arrays looping over the third dimension is indeed the fastest choice already. Here's the timing code:

n = 2;
m = 2;
l = 1e4;

A = rand(n,m,l);
B = rand(m,n,l);

% naive for-loop:
tic
%Cf = nan(n,n,l);
Cf = A;
for jl = 1:l
    Cf(:,:,jl) = A(:,:,jl) * B(:,:,jl);
end;
disp([' array-for: ' num2str(toc)]);

% using arrayfun:
tic
Ca = arrayfun(@(k) A(:,:,k)*B(:,:,k), 1:size(A,3), 'UniformOutput',false);
Ca = cat(3,Ca{:});
disp([' arrayfun : ' num2str(toc)]);

tic
Ac = num2cell(A, [1 2]);
Bc = num2cell(B, [1 2]);
disp([' num2cell : ' num2str(toc)]);

% cell for-loop:
tic
Cfc = Ac;
for jl = 1:l
    Cfc{jl} = Ac{jl} * Bc{jl};
end;
disp([' cell-for : ' num2str(toc)]);

% using cellfun:
tic
Cc = cellfun(@mtimes, Ac, Bc, 'UniformOutput', false);
disp([' cellfun  : ' num2str(toc)]);

tic
Cc = cell2mat(Cc);
disp([' cell2mat : ' num2str(toc)]);

tic
Cm = A;
Cm(1,1,:) = A(1,1,:).*B(1,1,:) + A(1,2,:).*B(2,1,:);
Cm(1,2,:) = A(1,1,:).*B(1,2,:) + A(1,2,:).*B(2,2,:);
Cm(2,1,:) = A(2,1,:).*B(1,1,:) + A(2,2,:).*B(2,1,:);
Cm(2,2,:) = A(2,1,:).*B(1,2,:) + A(2,2,:).*B(2,2,:);
disp([' explicit : ' num2str(toc)]);

disp(' ');