Эффективный способ умножения ячеек и укладки разреженных матриц в MATLAB - PullRequest
Новая
       14

Эффективный способ умножения ячеек и укладки разреженных матриц в MATLAB

0 голосов
/ 12 ноября 2018

У меня есть ячейка J размера 1xs, заполненная разреженными матрицами размера mxn (m>=n) и двумя полными матрицами r и l размера mxcxp и sxcxp соответственно. Размеры примерно 

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

То, что я делаю до сих пор, так это (см. Код ниже), но это, похоже, довольно неэффективно. Я ищу способ, как ускорить процесс масштабируемым образом, так как мне нужно многократно выполнять эту операцию (также m может увеличиться до 300 КБ), поэтому было бы здорово иметь этот шаг как можно быстрее. , поскольку это кажется узким местом до сих пор. Вот код того, что мне нужно сделать: 

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(J{1}),c*p);
col_vec = zeros(nnz(J{1}),c*p);
val_vec = zeros(nnz(J{1}),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_ + J{si} * 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);

Кроме того, я нашел этот путь без вложенных for -циклов, но это, кажется, еще менее эффективно: 

% 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_(:)));

Первая версия занимает примерно 0,27 с, а вторая - около 0,3 с.

...