# Multiply matrices in cell array by another matrix

11 views (last 30 days)
Ivan Bioli on 27 Oct 2023
Commented: Matt J on 29 Oct 2023
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

James Tursa on 27 Oct 2023
Edited: James Tursa on 27 Oct 2023
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 on 27 Oct 2023
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.

Matt J on 27 Oct 2023
Edited: Matt J on 27 Oct 2023
AU=blkReshape( vertcat(A{:})*U, [n,r] ,1,1,[]);

Matt J on 27 Oct 2023
Edited: Matt J on 27 Oct 2023
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 on 29 Oct 2023
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 on 29 Oct 2023
Yes, indeed.