Efficient summing of parts of an array

I have a 4D dataset that I'm expressing as a sparse 2D array of 1s and 0s. The 4D data has dimensions NPixHeight by NPixWidth by NBins by NFrames, and as a sparse array has the dimensions NPixHeight*NFrames by NPixWidth*NBins.

As part of my pre processing of I want to compress the 2D array by summing along the F dimension in groups, the size of the group is defined by nExps in the below code. To realise this I create and transform matrix T that after multiplication with the original array outputs a NPixHeight*nExps by NPixWidth*NBins array.

This last step is pretty slow for our planned application and I was hoping someone might be able to point me in the right direction to speed it up.

Here's an excerpt from the code

%% Data Reshape
% number of exposures for macro time
nExps = NFrames/macroFrames;
% generates sparse array of counts Ypix*NFrames by Xpix*Bins where
% NFrames is limited by macrotime
yNew = y+(f-1).*NPixHeight;
xNew = x+(b-1).*NPixWidth;
SpCounts = sparse(yNew,xNew,1,NPixHeight*NFrames,NPixWidth*NBins); % SLOW
% transform matrix that sums a number frames set by macro time
T = repmat(speye(NPixHeight,NPixHeight),1,macroFrames);
TCell = repmat({T}, 1, nExps);
clearvars T
T = blkdiag(TCell{:});
% multiplication and reshaping to output 
% 4D array of Y x X x Bins x NumberMacroExposures
camData.counts = sparse(size(SpCounts,2),size(T,1));
camData.counts = full(T*SpCounts).'; % SLOW
camData.counts = reshape(camData.counts,NPixWidth*NBins,NPixHeight,nExps);
camData.counts = reshape(pagetranspose(camData.counts),NPixHeight,NPixWidth,NBins,nExps);  

(I know the double reshape is awkard but its not enough of a bottleneck to bother fixing)

[James Tursa]

What are the typical values of your size variables? What is the typical sparsity of your sparse matrix? I would have to experiment with a small example set of data to confirm what is going on, but often it is the deep data copies that can drag the performance, and sometimes the only way to avoid that is to resort to a mex routine.

[Ife]

It depends on the user input but for a "worst case" with:

NPixHeight = 192

NPixWidth = 126

NBins = 344

NFrames = 247360

macroFrames = 30920

nExps = 8

gives a sparse array of 47493120 x 43344 with a sparsity of 6.2715e-05.

Forgive the stupid question but what is a mex routine...?

[the cyclist]

A MEX file is a function that calls C,C++, or Fortran from within MATLAB.

Answer 1

[Matt J]

I would recomend that you maintain the 4D array in ndSparse form instead,

https://www.mathworks.com/matlabcentral/fileexchange/29832-n-dimensional-sparse-arrays

at which point you can do,

temp = reshape(yourArray,[],NFrames)*kron(speye(nExps),ones(macroFrames,1));

result=reshape(temp,NPixHeight, NPixWidth , NBins, []);

[Ife]

Thanks! Due to what I'm using this code for I can skip the 4D sparse representation then reshape steps and readout my data directly into the correct size/orientation. My code is running 2x - 5x faster depending on the size of the datasets now!