Multiply matrices in cell array by another matrix

Hi,

I have a matrix U of size [n, r] and cell array A of length m such that A{i} is a matrix of size [n, n] for every i. I would like to obtain a matrix AU of size [n, r, m] such that AU(:, :, i) = A{i} * U. I cannot save A as a 3d matrix and use pagemtimes because each A{i} is a sparse matrix. For the moment I just used a naive for loop

AU = zeros(n, r, m)
for i = 1:m
    AU(:, :, i) = A{i} * U;
end

Is there a more efficient and/or more compact way of doing this?

Thanks,

Ivan

Answer 1

[James Tursa]

You may be stuck with the loop. There are ways to rearrange and stack things so that you can do everything in a single matrix multiply, but this will involve deep data copies of the sparse matrices and maybe you don't want that. What are the sizes involved? I.e., what are typical values of n, r, and m? Is U full or sparse?

*** EDIT ***

E.g.,

n = 3;
m = 3;
r = 2;
% generate sample inputs
A = arrayfun(@(k)sparse(rand(n)),1:m,'uni',false);
U = rand(n,r);
% Looping method
AU = zeros(n, r, m);
for i = 1:m
    AU(:, :, i) = A{i} * U;
end
% Single matrix multiply method
AC = vertcat(A{:});
ACU = AC * U;
C = mat2cell(ACU,n*ones(m,1),r);
ACU = cat(3,C{:});
disp(AU)

(:,:,1) =

    0.7145    0.7748
    0.2497    0.2794
    0.4050    0.4254


(:,:,2) =

    0.8667    0.8813
    0.5480    0.6060
    0.2374    0.2379


(:,:,3) =

    0.1883    0.1916
    0.9469    1.0044
    0.5184    0.5348

disp(ACU)

(:,:,1) =

    0.7145    0.7748
    0.2497    0.2794
    0.4050    0.4254


(:,:,2) =

    0.8667    0.8813
    0.5480    0.6060
    0.2374    0.2379


(:,:,3) =

    0.1883    0.1916
    0.9469    1.0044
    0.5184    0.5348

disp(max(abs(AU(:)-ACU(:))))

     0

So both methods can get the same result, but both methods require some deep data copying. You would have to run this with your actual variable sizes and sparsity to see if there is any benefit for the 2nd method. But my gut is there will be more data churning in the 2nd method giving you no performance benefit. The mat2cell( ) and cat( ) stuff could probably be combined into one step if I was more clever (easy to do in a mex routine in one step), saving you some time.

[Ivan Bioli]

Thanks for your answer. To be more precise, n and m range from 1'000 to 50'000, roughly speaking, and you can assume that n=m. r is much lower than both n and m, say from 10 to 200. U is full.

Answer 2

[Matt J]

Using this FEX download,

https://www.mathworks.com/matlabcentral/fileexchange/115340-block-transposition-permutation-and-reshaping-of-arrays

AU=blkReshape( vertcat(A{:})*U,  [n,r] ,1,1,[]);

Answer 3

[Matt J]

If instead of creating A{i}, you can instead create a cell array B such that B{i}=A{i}.' , then it would be much more efficient,

AU = pagetranspose( reshape(  U.'*[B{:}]  ,[r,n,m])  )

because column-wise concatenation of sparse matrices is always faster.

[Ivan Bioli]

Thanks for your answer. The matrices involved are symmetric, hence B = A in my case. By doing [B{:}], aren't we copying the matrices B under the hood?

[Matt J]

Yes, indeed.