Pfthb

What I have is a cell J of size 1xs filled with sparse matrices of size mxn (m>=n), and two full matrices r and l of size mxcxp and sxcxp, respectively. The dimensions are roughly

s = 4; % or 9
m = 10000; % can increase up to 300k
n = 36; % can increase up to m
c = 3;
p = 25;

What I do until now is this (see code below), but this seems to be quite inefficient. I’m looking for a way on how to speed up things in a scalable way, as I need to do this operation many times (also m can increase up to 300k), thus having this step as fast as possible would be great, as this seems to be the bottleneck till now. Here the code on what I need to do:

J = cell(s, 1);
% just fill J with sparse matrices of size mxn. Each sparse matrix is different for each cell, but all have the same nnz.
J(:) = sprand(m,n,0.1);
r = rand(m, c, p);
l = rand(s, c, p);

% preallocate vectors
row_vec = zeros(nnz(J1),c*p);
col_vec = zeros(nnz(J1),c*p);
val_vec = zeros(nnz(J1),c*p);

% do computation
for pi = 1:p
 for ci = 1:c
 J_ = 0;
 for si=1:s % multiply each sparse matrix in cell with scalar l(si,ci,pi) and sum them up 
 J_ = J_ + Jsi * l(si,ci,pi);
 end
 % multiply resulting sparse matrix with diagonal matrix (resulting from vector r(:,ci,pi)) and get final indices for later
 [row_vec_temp, ...
 col_vec(:,(pi-1)*c+ci), ...
 val_vec(:,(pi-1)*c+ci)] = find(spdiags(r(:,ci,pi),0,m,m) * J_);
 row_vec(:,(pi-1)*c+ci) = row_vec_temp + row_vec(end,max(1,(pi-1)*c+ci-1));
 end
end
% build final stacked sparse matrix of size m*c*pxn using calculated indices.
J_final = sparse(row_vec, col_vec, val_vec);

Additionally, I found this way without nested for-loops, but this seems to be even less efficient:

% create cell 1xcxp cell from r, where each cell is mx1 vector
R = mat2cell(r, m, ones(c,1), ones(p,1));
% make each cell a sparse diagonal matrix
R = repmat(cellfun(@(x) spdiags(x,0,m,m), R, 'UniformOutput', false), s, 1, 1);
% do point-wise computation
C = cellfun(@(x,y,z) z*(x.*y), repmat(J',1,c,p), mat2cell(l,ones(s,1),ones(c,1),ones(p,1)), R, 'UniformOutput',false);
% sum over each row of resulting cell C
J_ = cell(1,c,p);
J_(:) = 0;
for si=1:s
 J_(1,:,:) = cellfun(@(x,y) (x+y), J_(1,:,:), C(si,:,:), 'UniformOutput',false);
end
% stack final cell and convert it to sparse matrix
J_final = cell2mat(cat(1,J_(:)));

The first version takes roughly 0.27s and the second version takes about 0.3s.