How to multiply N matrices without a FOR loop? (Slices of 3D array)

I have a 3D matrix 2x2xN which, for my purposes, are essentially N 2x2 matrices and I want to do matrix multiplication with all of them so that I would get the following result:

N = 14;
M = rand(2,2,N);
Z = M(:,:,1)*M(:,:,2)* ... *M(:,:,N);
size(Z) == [2 2]

I can do it with a for loop, but I am looking for a single line approach, something like:

prod(M,3);

but probably with mtimes that would do matrix multiplication along the 3rd dimension (not the element-wise product).

I also converted matrix M into a Nx1 cell array of 2x2 matrices, but this approach did not work either to do the multiplication.

[KL]

Check this one: https://www.mathworks.com/matlabcentral/answers/62382-matrix-multiply-slices-of-3d-matricies

[Edgars Stegenburgs]

In that problem, you are processing two 3D matrices and the resulting matrix is 3D too. I have a subset of N matrices in an array (of 3D) or cell array and want to get 2x2 matrix at the end. I do not want to include or depend on additional libraries as it would make the whole thing bulkier too.

[Stephen23]

@Edgars Stegenburgs: what is actually wrong with using a loop?

[Edgars Stegenburgs]

While my whole code is currently vectorized except this one for loop, I just wonder if the code could run faster due to the fact that in the for loop an array is being accessed on each iteration, and a slice is taken of it.

[Stephen23]

@Edgars Stegenburgs: I doubt that it would make much difference. The slices would have to be taken either inside that operator, or explicitly in the loop as you do now: JIT code acceleration could possibly make it run almost as fast as that hypothetical internal loop (of course this means you need to write functions (which you should be doing anyway!)).

If efficiency really is your main aim then:

1. note that vectorized code is not always faster than loops: e.g. when intermediate arrays are large then it can be slower than looped code.
2. optimizing like this is often a waste of time. Timing of processes can change between MATLAB versions, so creatively tweaking your code to get 1% speedup out of some operator could easily turn out to be counter-productive by using the next MATLAB version or changing hardware.

I would suggest that you focus on writing code that is readable, clear, and easy to debug. This will save you more time in the end. The loop is fine.

[Andrei Bobrov]

@Stephen: +1

[Edgars Stegenburgs]

Thanks for the comment, I will definitely keep this in mind. I may experiment with the sizes of my matrices to see when it is faster by using for loops, but so far it was much slower.

Of course, in either case, you need to take the slices. If there was a MATLAB function that would do what I wanted, maybe it could have some ingenious way how it is handling the work, just like MATLAB is very fast in processing arrays.

Because execution time is fluctuating slightly, I always benchmark self-repeated code 10x,100x,..., depending on the initial required time for execution, and see if there is a noticeable difference, and because this as well is not constant, I do it several times.

In the end, I am maintaining my code and keeping it in the simplest form, quite readable, and keeping it short as well, that helps for debugging purposes, seeing more at a glance, which in turn somewhat improves understanding the code. Elegant...

[Jan]

Stephen's comment is very good.

For the estimation of the effects of optimizing the code, the usual sizes of the inputs matter: Is it really a [2 x 2 x N] array and what sizes of N do you have? For larger rows and columns, the main is done by mtimes, while the loop does not matter much. mtimes calls optimized BLAS or ATLAS functions, such that there is no room for further improvements. But I do not know, if these library function handle tiny 2x2 matrices with unrolled loops. So perhaps a C-Mex function could be more efficient.

Answer 1

[Jan]

If you really have 2x2 sub matrices to accumulate, try a C-Mex function:

#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
 const mwSize *size;
 mwSize N;
 double *p, *q, q11, q12, q21, q22, t11, t21;

   p    = mxGetPr(prhs[0]);
   size = mxGetDimensions(prhs[0]);
   if (size[0] != 2 || size[1] != 2) {
      mexErrMsgIdAndTxt("JSimon:CumMProd2x2:BadInput1",
              "1st input must be a [2 x 2 x N] array.");
   }
   N = size[2];

   q11 = p[0];
   q21 = p[1];
   q12 = p[2];
   q22 = p[3];
   while (--N) {  // Unrolled 2x2 matrix multiplication
      p  += 4;
      t11 = q11 * p[0] + q12 * p[1];
      t21 = q21 * p[0] + q22 * p[1];
      q12 = q11 * p[2] + q12 * p[3];
      q22 = q21 * p[2] + q22 * p[3];
      q11 = t11;
      q21 = t21;
   }

   plhs[0] = mxCreateDoubleMatrix(2, 2, mxREAL);
   q       = mxGetPr(plhs[0]);
   q[0]    = q11;
   q[1]    = q21;
   q[2]    = q12;
   q[3]    = q22;

   return;
  }

[EDITED] This is tested now. The speed is very interesting:

function speed

x = rand(2, 2, 1000);
tic; for k = 1:1000, y = CumMProd2x2(x); end; toc
tic; for k = 1:1000, y = CumMProd2x2_AB(x); end; toc
tic
for k = 1:1000  % Jos (10584)
   iif    = @(varargin) varargin{2*find([varargin{1:2:end}], 1, 'first')}() ;
   mprodf = @(F,M,n) iif (n < 2, M(:,:,1), true, @() F(F,M,n-1) * M(:,:,n)) ;
   out    = mprodf(mprodf, x, size(x, 3));
end
toc

end

function out = CumMProd2x2_AB(M)   % Andrei Bobrov
s   = size(M, 3);
out = M(:,:,1);
for ii = 2:s
   out = out * M(:,:,ii);
end
end

R2016b/64/Win7:

 Elapsed time is  0.011403 seconds.   C-mex
 Elapsed time is  3.884977 seconds.   Loop
 Elapsed time is 96.038754 seconds.   Recursive anonymous function

I was surprised, that Andrei's loop is such slow, although it is clearly the nicest and cleaned solution. Let's try to unroll the loops like in the C-Code:

function out = CumMProd2x2_unroll(M)
q11 = M(1);
q21 = M(2);
q12 = M(3);
q22 = M(4);
c   = 1;
for ii = 2:size(M, 3)
   c   = c + 4;
   t11 = q11 * M(c)   + q12 * M(c+1);
   t21 = q21 * M(c)   + q22 * M(c+1);
   q12 = q11 * M(c+2) + q12 * M(c+3);
   q22 = q21 * M(c+2) + q22 * M(c+3);
   q11 = t11;
   q21 = t21;
end
out = [q11, q12; q21, q22];
end

This 64 times faster than the direct approach "out * M(:,:,ii)":

 Elapsed time is 0.061287 seconds.  Unrolled

Obviously Matlab calls very smart highly optimized libraries for the matrix multiplication, which treat the tiny input with the same hammer method as a 1000x1000 matrix.

But this unrolled version is such ugly, that I would hesitate to use it in productive code. For x = rand(2, 2, 100000) I get the timings for 1000 iterations:

 Elapsed time is 1.377695 seconds.  C-mex
 Elapsed time is 2.872356 seconds.  M with unrolled mtimes

Only a factor 2! Another example, that loops are not such bad in Matlab compared to C.

[Andrei Bobrov]

@Jan: It's great! +1

[Jos (10584)]

haha, I really liked my anonymous function approach, and did expect it to perform poorly, but that poor ... haha

Answer 2

[Andrei Bobrov]

s = size(M)
out = M(:,:,1);
for ii = 2:s(3)
    out = out*M(:,:,ii);
end

[Edgars Stegenburgs]

Thanks for this answer but I already knew this approach, I want to find more vectorized method.

[Matt J]

I don't believe there is a vectorized method, however, a parfor loop could be used if you have the Parallel Computing Toolbox

parfor ii = 2:s(3)
  out = out*M(:,:,ii);
end

[Edgars Stegenburgs]

parfor cannot be used for matrix multiplication since it has to be done in the original sequence, the following relation is not true with matrices: A*B=B*A.

[Jan]

+1: This is the nicest solution. That the multiplication of 2x2 matrices is much faster with hard coded algorithm is not a problem of this solution.

Although the C-Mex approach is faster, it would be very hard to generalize it for inputs beside 2x2xN arrays.

[Matt J]

Although the C-Mex approach is faster, it would be very hard to generalize it for inputs beside 2x2xN arrays.

Just wanted to note that, while my solution based on MTIMESX is not as fast as Jan's for the 2x2xN case, it is applicable to arbitrary MxMxN arrays,

Answer 3

[Matt J]

The following is not a one-line solution (for that just stick it in a function file) and requires MTIMESX from the File Exchange. However, I do see a few factors speed-up over a conventional for-loop,

out=M;

while size(out,3)>1

      n=size(out,3);

      if mod(n,2)
          n=n-1;
       A=out(:,:,1:2:n);
       B=out(:,:,2:2:n);         
        out=cat(3,mtimesx(A,B),out(:,:,n+1));  
      else
       A=out(:,:,1:2:n);
       B=out(:,:,2:2:n);
       out=mtimesx(A,B);

      end

end

[Edgars Stegenburgs]

I think instead of:

out=cat(3,mtimesx(A,B),out(:,:,n+1));

you should have:

out=mtimesx(A,B)*out(:,:,n+1);

to end up with a 2D array.

I also thought of a similar way. Will have to check mtimesx finally.

[Matt J]

you should have: out=mtimesx(A,B)*out(:,:,n+1);

No, that will give you a dimensions mismatch error. I tested the code for agreement with a plain for-loop.

However, maybe you meant this

       A=out(:,:,1:2:n);
       B=out(:,:,2:2:n);
         B(:,:,end)=B(:,:,end)*out(:,:,n+1);

       out=mtimesx(A,B);

If so, I agree that that should be more efficient.

[Matt J]

When I use the mtimesx SPEED flag, things actually get considerably faster:

 N=1e7;
 M=rand(2,2,N);
   M=M./sum(M,2);

 %%%Baseline implementation
  tic; 
    out0=eye(2);
    for k=1:N
        out0=out0*M(:,:,k);
    end
  toc;

Elapsed time is 12.285774 seconds.

%%%With mtimesx

mtimesx SPEED

tic
out=M;
n=size(out,3);
while n>1

      if mod(n,2)
        C=out(:,:,n);  
         n=n-1;
        A=out(:,:,1:2:n);
        B=out(:,:,2:2:n); 
          B(:,:,end)=B(:,:,end)*C;
        out=mtimesx(A,B);

      else  
       A=out(:,:,1:2:n);
       B=out(:,:,2:2:n); 
       out=mtimesx(A,B);
      end

      n=size(out,3);

end
toc

Elapsed time is 0.785595 seconds.

[James Tursa]

Side Note: MTIMESX by default calls BLAS library routines for matrix multiply so that it matches MATLAB for-loop m-code result, whereas MTIMESX with the 'SPEED' option will use hand-coded inline matrix multiply code for up to 5x5 size slices which may not match MATLAB for-loop m-code result exactly.

Sometime back I had a beta version of MTIMESX that implemented the matrix equivalent versions of 'prod' and 'cumprod'. Maybe it is time I dust that off and finish the implementation/testing so I can publish it.

[Matt J]

That is strange, since I still see significant speed-up even with

mtimesx MATLAB

Answer 4

[Steven Lord]

If you're using release R2020b or later, take a look at the pagemtimes function introduced in that release.

Answer 5

[Jos (10584)]

Here is one using recursion without a for-loop; not faster though, and somewhat mysterious, but just nice :) ...

M = randi(5,[2 2 4]) ;  % data
iif = @(varargin) varargin{2*find([varargin{1:2:end}], 1, 'first')}() ;
mprodf = @(F,M,n) iif (n < 2, M(:,:,1), true, @() F(F,M,n-1) * M(:,:,n)) ; 
out = mprodf(mprodf,M,size(M,3))  % voila, it works!

[Edgars Stegenburgs]

Yep, this is really mysterious :D

[Jos (10584)]

It is the inline version of this recursive m-file:

    function X = mprod(M,n)
    % X = mprod(M) returns M(:,:,1) * M(:,:,2) * ... * M(:,:,end)
    % where M is a 3D array
    if nargin==1
       X = mprod(M,size(M,3)) ;
    elseif n < 2
      X = M(:,:,1) ;
    else
      X = mprod(M,n-1) * M(:,:,n) ; 
    end

[Edgars Stegenburgs]

Look very nice, albeit, I guess this will be slower due to recursivity.