Fast way of computing all possible products of n matrices

Hi, I have some n matrices (stored in a cell array) and want to compute all possible n! (max-plus-algebra-)products you can get by reordering them and everytime save the element in the bottom left corner. At the moment, I'm using the straight forward approach by iterating through all permutations and just calculating the product from start to end, which looks something like this:

for j = 1:n                             %initialise the n matrices               
    matrices{j} = some matrix;
end

interestingValue = zeros(1, factorial(n));  %array where to save the values

for i = 1:factorial(n)                    %iterate through all schedules 
    schedule = oneperm(n, i);             %get the current schedule (order)
    interestingValue(i) = computeProduct(matrices, schedule, n);         %compute product for current schedule
end

with the computeProduct function looking like this:

function value = computeProduct(matrices, schedule, n)

    product = matrices{schedule(1)};                  %initialise with first matrix
    for j = 2:n
        product = max_plus_product(product, matrices{schedule(j)});              %compute product of current product with next matrix
    end
    
    value = product(end, 1);    

My question would be how to speed up the calculation because performance is important for the project. I thought this approach I'm using right now is not really efficient because e.g. for all (n-2)! orders that would start with the matrices 1 and 2, you always compute the product of matrix 1 with matrix 2, instead of computing it just once because its always the same. So I thought about storing some intermediate results and using them later but I dont really know how to do this the most efficient way. Thats why I wanted to ask if there maybe already is a MATLAB function which is useful for this or if there is an other clever way.

Edit: mathemactical clarification for max-plus-multiplication: -> product of two n*n matrices A and B

Thank you very much for your help ^^

[Fifteen12]

Another way to state this question:

Given n square matrices, what is the fastest way to compute the products of all possible combinations? For instance, given matrices A, B, and C, what is the fastest way to compute {A*B*C, A*C*B, B*A*C, B*C*A, C*A*B, and C*B*A}?

[Bjorn Gustavsson]

If you only need the lower-left element then that is something that ought to significantly help speeding things up. You should really sit down and sketch this out on paper and think about what calculations you really need, then how to organize them efficiently.

[Matt J]

Given n square matrices, what is the fastest way to compute the products of all possible combinations?

It really depends on n, and on the sizes of the matrices.

[Matt J]

In a max-plus algebra, aren't products the same as sums? If so, why does the order of the matrices matter?

[Jonas Peschel]

Yes, you're right, standard multiplication becomes addition and standard addition becomes max operation. Still if you change the order in a matrix multiplication you compute different things just as in standard algebra.

[Matt J]

Still if you change the order in a matrix multiplication you compute different things just as in standard algebra.

I don't see how that's true, since if matrix multiplication is now standard matrix addition and since the latter is commutative - the order shouldn't matter. Additionally, if you are only interested in the lower left element of the result, it should be enough to sum the lower left elements of all the matrices independently of all other matrix entries.

I admit this is dealing with mathematics I'm not well-versed in, but I would guess almost all other forum contributors are going to be in the same boat. So, to get useful responses, you would probably have to provide some mathematical background.

[Jonas Peschel]

Hi, I updated the question to clarify the max-plus matrix multiplication a bit. It's not like you assumed that matrix multiplication becomes just adding the matrices element-wise. Anyway I have a local function to take care of that multiplication and I just wanted to ask for an efficient solution MATLAB-wise to compute all the possible products (using the local function and not standard MATLAB matrix multiplication) without maybe computing some specific products over and over as I stated in the question. I already tried to simply add memoization for the multiplication function to store the outputs and avoid redundant computations but thats not faster for me because one single multiplication for relatively small matrices is already really fast.

I hope my question makes a bit sense now thank you for any kind of help :)

[Bruno Luong]

Is max-plus associative? Meaning

(A*B)*C is equal to A*(B*C)

[Jonas Peschel]

Hi @Bruno Luong, yes max-plus matrix multiplication is associative.

[Bruno Luong]

My impression is if you bookeeping well you could reduce to two (and not n-1) max_plus_products by combinations using older results.

The problem is to store them if n is large. WHat are typical value of n?

If it is reasonable I might take a closer look.

[Jonas Peschel]

Hi, at the moment a typical value of n would be maybe 9,10 but it would be desirable to also be able to go a bit higher maybe.

Answer 1

[Matt J]

This solution (note the attachments) does it in about 1 sec for a 20x20x10 set of single float matrices.

matrices=rand(20,20,10,'single');

tic;

    M=matrices;
    Mt=pagetranspose(M);
    
    [N,~,n]=size(matrices);
    
    s=ceil(n/2); t=n-s;
    
    e=uint8(1:n);
    Pleft=nchoosek(e,s); 
    Pright=flipud(nchoosek(e,t));
    
    
    p=height(Pleft);
    
    result=cell(1,p);
    
    for i=1:p
    
      matricesLeft=M(:,:,Pleft(i,:));  
      matricesRight=Mt(:,:,Pright(i,:)); 
    
      partialLeft=permProducts(matricesLeft,N);
      partialRight=permProducts(matricesRight,1).';
    
      tmp=pageMaxPlus(partialLeft,partialRight);
    
      result{i}=tmp(:).';
    end


   result=cell2mat(result).';
   
toc;

Elapsed time is 1.043404 seconds.

[Jonas Peschel]

Hi Matt, thank you really much for the answer :)

I'm sorry I didnt stress this much in the question but to me it looks like the pagemtimes function only uses standard multiplication where I wanted to use a local funtion (i.e. max-plus-algebra matrix multiplication). Is there a similar way to yours but where you can use your own function for the multiplication?

Also you said the answer to the question depends on the number of matrices and their size, for me the maximum number of matrices would be maybe 10 or so and the size about 20*20. Would this be still small enough for your approach to be efficient or why not?

Thanks a lot!

[Jonas Peschel]

Thank you very much for the answer and your effort. I'm sorry I dont really fully understand your code and also when I try to run it it gives me an error for this line: p=reshape(p,[],P/chunksize,[]); . But overall to me it looks like your method is still basically doing what a loop would do and compute every possible permutation from start to end, isnt it? Just in a more fancy MATLAB way, sorry if this is completly not the case and I just dont understand it. So I was wondering if this is really much faster than the loop and also like you stated the amount of matrices to store is really huge. What I had in mind was more like how to avoid doing all the redundant computations by storing some intermediate results e.g. start by computing every schedule starting with matrices 1, 2 but not every time recomputing matrix1 * matrix2 and so on. I dont know if this makes any sense and I'm sorry if I just completly misunderstood your code.

[Matt J]

I'm sorry I dont really fully understand your code and also when I try to run it it gives me an error for this line: p=reshape(p,[],P/chunksize,[]);

Sorry, I've updated my answer with a fully tested version.

But overall to me it looks like your method is still basically doing what a loop would do and compute every possible permutation from start to end, isnt it?

I believe mine is much more vectorized, since your method has a loop of factorial(n) iterations, whereas mine has loops of only O(n) iterations (for n<=16).

In any case, how fast is your implementation? How fast where you hoping it would be? For a 20x20x10 set of single float matrices, mine runs in about 1 second.

What I had in mind was more like how to avoid doing all the redundant computations by storing some intermediate results

I've implemented that in my other answer, but I find that it is considerably slower. The thing to remember is that we live in an age of multicore computing. So, it is not how many operations you do that determines performance exclusively, but also how well you can parallelize the computations. In Matlab, it can be beneficial to do redundant computation in the interest of being able to take advantage of its highly multi-threaded matrix arithmetic.

[Jonas Peschel]

Wow, thank you so much for your very detailed answers and your effort, thank you. This is completly awesome, I didnt realise your implementation is that much faster, for n=9 already my looping takes like 20s and yours is like 0.5s this is very impressive to me. I get that your version is somehow more vectorized and more "MATLAB-style" or something but still I didnt know the loop was that bad and I dont really understand why yours is soo much better. I guess I'm going to spend some time now on trying to understand what your code does exactly ^^ because for the project I may need some slightly modified version. Thank you so much again for helping, this really surprised me that you can improve the performance of this so much, also I'm feeling a little bad about myself now knowing that my whole code could probably be much faster if written by you xD

Thanks

[Matt J]

already my looping takes like 20s

I'm suprised plain vanilla looping is that fast!

[Bruno Luong]

Very iimpressive timing

it takes Elapsed time is 0.473249 seconds 

on my PC.

However I have two questions:

• why not push the divide and conquer to more than first level?
• It seems there are still redundancy between partialLeft and partialRight, can it be resuced?

[Matt J]

why not push the divide and conquer to more than first level?

Not exactly sure how you would do that, or if it would be worthwhile if you're not going to go much beyond n=10.

It seems there are still redundancy between partialLeft and partialRight, can it be resuced?

The only way I can think of to combine divide-and-conquer with redundancy-elimination is to re-implement permProducts with the recursive method in my other answer. But I consistently see that recursion costs more than it saves.

[Bruno Luong]

For some reason this code returns 120 data instead of 6 for n=3.

[Jonas Peschel]

I guess you figured it out yourself but maybe for anyone else reading this or wanting to use the code, I think there is a little bug in the code, so that it doesnt work for number of matrices <= 3. I think its because the squeeze function at the end of permProducts doesnt work as intended in these cases if the third dimension is also 1, so that results in 'partialRight' or 'partialLeft' being the transposed of what you really want. So maybe just add these two lines after computation of partialLeft and partialRight (which doesnt do anything in case n>3) and then everything works as intended:

partialLeft = reshape(partialLeft, [], N);
partialRight = reshape(partialRight, N, []);

Still the code doesnt work in this form for n=1 because nchoosek(1,0) returns 1 instead of zero but if you like you could easily fix this as well.

[Bruno Luong]

Yes I do, the usage of squeeze must be banned. Everytime I use it a hdden bug is just wait to bite me later.

[Jonas Peschel]

Hi, one quick question, I dont understand why you get different results when choosing a different partition in the beginning. E.g. when I use floor instead of ceil, so one matrix gets shifted to the matricesRight I obtain different results at the end. I dont get why this should be the case isnt the partition into matricesRight and matricesLeft arbitrary?

[Bruno Luong]

May be the order of the results are different but as set they are identical?

[Jonas Peschel]

Yeah, I thought so as well but they are really different not only in order.

[Bruno Luong]

Small experience with 3/2 split vs 2/3 split. The results look identical to me.

n=5;
% Fake data
M=rand(20,20,n,'single');
matrices=M;
tic;
M=matrices;
Mt=pagetranspose(M);

[N,~,n]=size(matrices);

s=ceil(n/2); t=n-s;

e=uint8(1:n);
Pleft=nchoosek(e,s); 
Pright=flipud(nchoosek(e,t));


p=height(Pleft);

result=cell(1,p);

for i=1:p
    
    matricesLeft=M(:,:,Pleft(i,:));  
    matricesRight=Mt(:,:,Pright(i,:)); 
    
    partialLeft=permProducts(matricesLeft,N);
    partialRight=permProducts(matricesRight,1).';
    
    tmp=pageMaxPlus(partialLeft,partialRight);
    
    result{i}=tmp(:).';
end
result=cell2mat(result).';

toc;

Elapsed time is 0.094892 seconds.

tic
M=matrices;
Mt=pagetranspose(M);

[N,~,n]=size(matrices);

s=floor(n/2); t=n-s;

e=uint8(1:n);
Pleft=nchoosek(e,s); 
Pright=flipud(nchoosek(e,t));


p=height(Pleft);

result2=cell(1,p);

for i=1:p
    
    matricesLeft=M(:,:,Pleft(i,:));  
    matricesRight=Mt(:,:,Pright(i,:)); 
    
    partialLeft=permProducts(matricesLeft,N);
    partialRight=permProducts(matricesRight,1).';
    
    tmp=pageMaxPlus(partialLeft,partialRight);
    
    result2{i}=tmp(:).';
end
result2=cell2mat(result2).';
toc

Elapsed time is 0.017183 seconds.

close all
plot(sort(result(:)),'bo')
hold on
plot(sort(result2),'r+')

%%
function AxB=pageMaxPlus(A,B)
%Page-wise implementation of max-plus multiplication
[mA,nA,pA]=size(A);
[mB,nB,pB]=size(B);
A = reshape(A, mA, nA,  1, pA);
B = reshape(B,  1, mB, nB, pB);
AxB = max(A+B,[],2);
AxB = reshape(AxB,mA,nB,max(pA,pB));
end

function result=permProducts(matrices, keyrow)

    [N,~,n]=size(matrices);
   
    if nargin<2, keyrow=N; end

    numchunks=max(1,factorial(n)/factorial(8)); 
    
    p=perms(1:n)';
      p=reshape(p,n,[],numchunks);
      P=width(p);
    
    result=cell(1,numchunks);  
    for c=1:numchunks
    
        M=reshape( matrices(:,:,p(:,:,c)) , N,n*N,P);
        
        M=mat2cell(M,N,N*ones(1,n),P);
    
        M{1}=M{1}(keyrow,:,:);

        for s=1:n-1
        
            M{s+1}=pageMaxPlus(M{s}, M{s+1});
        
        end
        
        result{c}=M{end}; %final result
    
    end
    
    result=squeeze(cell2mat(result)).';


end

[Jonas Peschel]

Yes, really sorry to bother you. It works now for me also as it should be, I deleted a transpose which was necessary ^^

Answer 2

[Matt J]

Here's an alternative method that eliminates redundant computation, like you were hoping, using recursion. However, it doesn't seem to give a net performance advantage. On my laptop, it takes 380 sec. for a 20x20x10 set of single float matrices, which is quite a bit slower than my other answer.

clear, clc

matrices=rand(20,20,10,'single');


tic;
n=size(matrices,3);


result=cell(1,n);

for i=1:n
 result{i}=recursionStep(matrices,matrices(:,1,i),i);
end

 result=cell2mat(result);

 result=result(end,:);

toc

Elapsed time is 380.316553 seconds.

function out=recursionStep(matrices,accum,sofar)

  n=size(matrices,3);

 if n==numel(sofar), out=accum; return ; end  %end recursion


  S=setdiff(1:n,sofar);  
  s=numel(S);


  accum=pageMaxPlus(matrices(:,:,S),  accum);

  Sofars=[repmat(sofar,s,1),S(:)];

  out=cell(1,s);
  for i=1:s

   out{i}=recursionStep(matrices,accum(:,:,i),Sofars(i,:));

  end

  out=cell2mat(out); out=out(:,:);


end

function out=pageMaxPlus(L,R)
%Page-wise implementation of max-plus multiplication

    [mL,nL,~]=size(L);
    [mR,nR,~]=size(R);   

    assert(mR==nL,'Size mismatch');

    if mL<nR, implem='left'; else, implem='right'; end

   switch implem

       case 'left'

         L=pagetranspose(L);

         out=R(1:mL,:,:);

         for i=1:mL
             out(i,:,:)=max(L(:,i,:)+R,[],1);
         end

       case 'right'

         R=pagetranspose(R);

         out=L(:,1:nR,:);

         for i=1:nR
             out(:,i,:)=max(L+R(i,:,:),[],2);
         end        

   end


end

[Bruno Luong]

I don't know if it's matter but this implementation of pageMaxPlus avoid the transpose

function out=pageMaxPlus(L,R)
%Page-wise implementation of max-plus multiplication
[mL,nL,~]=size(L);
[mR,nR,~]=size(R);   
assert(mR==nL,'Size mismatch');
L = reshape(L, mL, nL, 1, []);
R = reshape(R,  1, mR,nR, []);
out = max(L+R,[],2);
out = reshape(out,mL,nR, []);
end

On my computer the code takes 121 sec

Answer 3

[Bruno Luong]

The more I try my idea, the more Matt's solution seem amazing to me.

Here is an algorithm that that build a set of product of all combinations (order matter) of k+1 matrices; starting from k=1, 2, up to n-1. If I count correctly there are exp(1)*factorial(n) matrix products to carried out, which seem to be an advance from the brute force algoritm, which requires (n-1)*factorial(n), for n>=4. However the timing is far, far away from Matt's solution.

Hat off for him.

n=9;
% Fake data
M=rand(20,20,n,'single');

n = size(M,3);
% C wil contain the result, it will have size m x m x factorial(n)
% We start with C of "product" of single matrices
C = M;
% This logical array will growed in row and keeps track of which matrix that has
% been accumulated in C
tobedo = (1:n).'~=(1:n);
for k=1:n-1
    p = size(tobedo,1);
    CTmp = cell(p,1);
    for i=1:p % possible to do parfor
        tbd = tobedo(i,:);
        CTmp{i} = pageMaxPlus(C(:,:,i),M(:,:,tbd));
    end
    C = cat(3,CTmp{:});

    m = n-k;
    [i,j] = find(tobedo);
    [i,is] = sort(i);
    j = j(is);
    h = repmat((1:m).',p,1);
    tobedo = repelem(tobedo, m, 1);
    idxdone = h + m*(i-1+p*(j-1));
    tobedo(idxdone) = false;

end


%%
function AxB=pageMaxPlus(A,B)
%Page-wise implementation of max-plus multiplication
[mA,nA,~]=size(A);
[mB,nB,~]=size(B);   
A = reshape(A, mA, nA,  1, []);
B = reshape(B,  1, mB, nB, []);
AxB = max(A+B,[],2);
AxB = reshape(AxB,mA,nB, []);
end

[Matt J]

If I count correctly there are exp(1)*factorial(n) matrix products to carried out, which seem to be an advance from the brute force algoritm, which requires (n-1)*factorial(n), for n>=4. However the timing is far, far away from Matt's solution.

I'm not completely sure of the book-keeping, but here is one thing I notice:

% C wil contain the result, it will have size m x m x factorial(n)

So, my code doesn't compute an mxm matrix for each permutation. The OP said he is only interested in the lower left entry of the final matrices, so I throw away the first m-1 rows of the first matrix in the sequence, thereby reducing all the matrix products to vector-matrix products. This requires at least a factor of O(m) less computation and memory.

Also, because of the divide and conquer, there are (n-1)*factorial(n)/factorial(n/2) vector-matrix products.

Also, because the products are vector-matrix products, my implementation of pageMaxPlus with pagetranspose might be cheaper than the approach you take with implicit expansion. Transposition of a multi-page vector is really only a reshape() operation.

[Bruno Luong]

Oh thanks @Matt J; I miss the trick of not computing the whole matrix, here is the version with that simplification

n=9;
% Fake data
M=rand(20,20,n,'single');

n = size(M,3);
% This logical array will growed in row and keeps track of which matrix that has
% been accumulated in C
tobedo = true(1,n);
% Define what elements of the final cumulative matrices we want
selected_row = size(M,1); % EDIT fix bug
selected_col = 1;
% C wil contain the result, it will have size m x 1 x factorial(n)
% We start with C of "product" of single matrices
C = M(selected_row,:,:);
p = 1;
for k=1:n-1
    % Update the tobedo array
    m = n-k+1;
    [j,i] = find(tobedo.');
    h = repmat((1:m)',p,1);
    idxdone = h + m*(i-1+p*(j-1));
    tobedo = repelem(tobedo, m, 1);
    tobedo(idxdone) = false;

    p = p*m;
    CTmp = cell(p,1);
    for i=1:p % possible to do parfor
        tbd = tobedo(i,:);
        CTmp{i} = pageMaxPlus(C(:,:,i),M(:,:,tbd));
    end
    C = cat(3,CTmp{:});
end
% extract and reshape to get final result, row vector of length factorial(n)
result = reshape(C(:,selected_col,:), 1, []);


%%
function AxB=pageMaxPlus(A,B)
%Page-wise implementation of max-plus multiplication
[mA,nA,pA]=size(A);
[mB,nB,pB]=size(B);
A = reshape(A, mA, nA,  1, pA);
B = reshape(B,  1, mB, nB, pB);
AxB = max(A+B,[],2);
AxB = reshape(AxB,mA,nB,max(pA,pB));
end