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?
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(' ');