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

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:ticCa = arrayfun(@(k) A(:,:,k)*B(:,:,k), 1:size(A,3), 'UniformOutput',false);Ca = cat(3,Ca{:});disp([' arrayfun : ' num2str(toc)]);ticAc = num2cell(A, [1 2]);Bc = num2cell(B, [1 2]);disp([' num2cell : ' num2str(toc)]);% cell for-loop:ticCfc = Ac;for jl = 1:l    Cfc{jl} = Ac{jl} * Bc{jl};end;disp([' cell-for : ' num2str(toc)]);% using cellfun:ticCc = cellfun(@mtimes, Ac, Bc, 'UniformOutput', false);disp([' cellfun  : ' num2str(toc)]);ticCc = cell2mat(Cc);disp([' cell2mat : ' num2str(toc)]);ticCm = 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(' ');

Here is my benchmark test comparing the methods mentioned in @TobiasKienzler answer. I am using the TIMEIT function to get more accurate timings.

function [t,v] = matrixMultTest()    n = 2; m = 2; p = 1e5;    A = rand(n,m,p);    B = rand(m,n,p);    %# time functions    t = zeros(5,1);    t(1) = timeit( @() func1(A,B,n,m,p) );    t(2) = timeit( @() func2(A,B,n,m,p) );    t(3) = timeit( @() func3(A,B,n,m,p) );    t(4) = timeit( @() func4(A,B,n,m,p) );    t(5) = timeit( @() func5(A,B,n,m,p) );    %# check the results    v = cell(5,1);    v{1} = func1(A,B,n,m,p);    v{2} = func2(A,B,n,m,p);    v{3} = func3(A,B,n,m,p);    v{4} = func4(A,B,n,m,p);    v{5} = func5(A,B,n,m,p);    assert( isequal(v{:}) )end%# simple FOR-loopfunction C = func1(A,B,n,m,p)    C = zeros(n,n,p);    for k=1:p        C(:,:,k) = A(:,:,k) * B(:,:,k);    endend%# ARRAYFUNfunction C = func2(A,B,n,m,p)    C = arrayfun(@(k) A(:,:,k)*B(:,:,k), 1:p, 'UniformOutput',false);    C = cat(3, C{:});end%# NUM2CELL/FOR-loop/CELL2MATfunction C = func3(A,B,n,m,p)    Ac = num2cell(A, [1 2]);    Bc = num2cell(B, [1 2]);    C = cell(1,1,p);    for k=1:p        C{k} = Ac{k} * Bc{k};    end;    C = cell2mat(C);end%# NUM2CELL/CELLFUN/CELL2MATfunction C = func4(A,B,n,m,p)    Ac = num2cell(A, [1 2]);    Bc = num2cell(B, [1 2]);    C = cellfun(@mtimes, Ac, Bc, 'UniformOutput', false);    C = cell2mat(C);end%# Loop Unrollingfunction C = func5(A,B,n,m,p)    C = zeros(n,n,p);    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,:);end

The results:

>> [t,v] = matrixMultTest();>> tt =      0.63633      # FOR-loop      1.5902       # ARRAYFUN      1.1257       # NUM2CELL/FOR-loop/CELL2MAT      1.0759       # NUM2CELL/CELLFUN/CELL2MAT      0.05712      # Loop Unrolling

As I explained in the comments, a simple FOR-loop is the best solution (short of loop unwinding in the last case, which is only feasible for these small 2-by-2 matrices).

I highly recommend you use the MMX toolbox of matlab. It can multiply n-dimensional matrices as fast as possible.

The advantages of MMX are:

1. It is easy to use.
2. Multiply n-dimensional matrices (actually it can multiply arrays of 2-D matrices)
3. It performs other matrix operations (transpose, Quadratic Multiply, Chol decomposition and more)
4. It uses C compiler and multi-thread computation for speed up.

For this problem, you just need to write this command:

C=mmx('mul',A,B);

I added the following function to @Amro's answer

%# mmx toolboxfunction C=func6(A,B,n,m,p)    C=mmx('mul',A,B);end

I got this result for n=2,m=2,p=1e5:

    1.6571 # FOR-loop    4.3110 # ARRAYFUN    3.3731 # NUM2CELL/FOR-loop/CELL2MAT    2.9820 # NUM2CELL/CELLFUN/CELL2MAT    0.0244 # Loop Unrolling    0.0221 # MMX toolbox  <===================

I used @Amro's code to run the benchmark.