Construct product of matrices without for loop

Hi, everyone:

Suppose A is a 2 by 2 by N matrix, is it possible to construct a matrix B, such that:

B=[Prod(A)_1^{N-n},   Prod(A)_2^{N-n},   ... , Prod(A)_{N-n}^{N-n};  

   Prod(A)_1^{N-n+1}, Prod(A)_2^{N-n+1}, ... , Prod(A)_{N-n}^{N-n+1};
           .                       .          .              .        
           .                       .          .              .         
   Prod(A)_1^N,       Prod(A)_2^N,       ... , Prod(A)_{N-n}^N]

where Prod(A)_m^n = A(:,:,m)*A(:,:,m+1)*A(:,:,m+2)* ... * A(:,:,n); (1 << n << N)

without any loop?? (maybe use arrayfun? cellfun?)

because for loop is too slow for me.

Thanks

[Jan]

When you claim, that the FOR loop is too slow, do not hope, that CELLFUN or ARRAY is faster. But use loops internally also, but have the additional overhead for calling Matlab from the Mex level.

When your code is too slow for your needs, you could post the relevant part of it and ask, if somebody knows improvements. For testing the ideas before answering, provided test data are very useful.

"n" appears twice with different meaning in your example. How large is it in the definition of B and what are typical sizes of N?

[Walter Roberson]

This appears to duplicate part of http://www.mathworks.co.uk/matlabcentral/answers/84498-construct-product-of-matrices-without-for-loop

Answer 1

[Jan]

I assume a FOR loop is the best way to go. Perhaps the slow speed is caused by a forgotten pre-allocation or if you calculate the power operation by the expensive ^ operator instead of a cumulative multiplication.

C = zeros(2, 2, n+1, N-n);
for i2 = 1:N-n
  AA = A(:, :, i2);
  PA = AA ^ (N-n-1);
  for i1 = 1:n+1
    PA = PA * AA;
    C(:, :, i1, i2) = PA;
  end
end

But if N is large and n is small, the power operation will be the bottleneck again. Then try these functions written by James: FEX: MTimesX or FEX: MPower2.

Another approach could be to avoid calling the powerful BLAS function for the matrix multiplication for a tiny 2x2 matrix:

function MatrixMultTest
x = rand(2, 2);
y = rand(2, 2)
tic; for k = 1:1e6; v = mult1(x,y); end, toc
tic; for k = 1:1e6; v = mult1(x,y); end, toc

function v = mult1(x, y)
v = x * y

function v = mult2(x, y)
v(2,2) = x(2,1)*y(1,2) + x(2,2)*y(2,2);  % Implicit pre-allocation
v(1,1) = x(1,1)*y(1,1) + x(1,2)*y(2,1);
v(1,2) = x(1,1)*y(1,2) + x(1,2)*y(2,2);
v(2,1) = x(2,1)*y(1,1) + x(2,2)*y(2,1);

Well, it seems like the manual product is 10% slower under R2009a. But it was worth to try.

[Andrei Bobrov]

+1