Efficient method of summing the values in multiple images on GPU?

I have a series of images sitting on GPU in a 71x71x5000 matrix.

Is there an efficient method of finding the sum of all elements in each 71x71 image on the GPU?

I wanted to do something similar to this but it isn't so simple.

sums = pagefun(@sum, images)

Answer 1

[Brendan Hamm]

I don't have a machine with the GPU capabilities with me right now, but I think this should work:

sums = pagefun(@(x) sum(x(:)),images);

[Christopher Desmond]

Error using gpuArray/pagefun

The GPU does not support the specified function handle using PAGEFUN.

[Brendan Hamm]

Ah yes, this supports only element-wise operations. I'll have to check this later when I get home, but we might be able to get this done with arrayfun.

Another possibility is that we could do this with 2 calls to pagefun by performing an mtimes with a vector of all ones (essentially using matrix multiplication to replicate sum(sum(Q(:,:,k))).

[Brendan Hamm]

So arrayfun will also not work and you are left with 3 options that I could think of (not sure if this is efficient enough):

A = rand(3,3,4,'gpuArray');

• You can reshape the gpuArray so each page becomes a column and perform the sum:

sum(reshape(A,[],size(A,3))

• As I mentioned in the previous post we could express this as a multiplication with pagefun.

pageSum = pagefun(@mtimes,A,ones(size(A,2),1))
% One pagefun gives you the sums of the rows
% repeating again will give the sums for pages
pageSum = pagefun(@mtimes,ones(1,size(pageSum,1)),pageSum)

If you would rather you can write this in one line:

pageSum = ...
pagefun(@mtimes,ones(1,size(A,1)),pagefun(@mtimes,A,ones(size(A,2),1)))

• Lastly a simple for loop can work:

for k = size(A,3):-1:1
    tmp = A(:,:,k);
    sumk(k) = sum(tmp(:));
end

I wrote a function to time these methods:

function t = timeGpu(m,n)
% Test the time different methods take to calculate the sum of the elements
% for each page on a GPU array.

% Create a gpuArray m-by-m-by-n
A = rand(m,m,n,'gpuArray');
% Allocate space for timing
t = nan(3,1);
% Time each method
t(1) = gputimeit(@() reshapeIt);
t(2) = gputimeit(@() pageIt);
t(3) = gputimeit(@() loopIt);

t = table(t,'RowNames',{'Reshape','pageIt','loopIt'});

function reshapeIt
    % Method one: use reshaping and sum
    ret = sum(reshape(A,[],size(A,3)));
end
function pageIt
    % Nested calls to pagefun with a matrix multiplication
    % ones'*A*ones, where ones is a column vector.
    vec1 = ones(size(A,2),1);
    ret = pagefun(@mtimes,vec1.',pagefun(@mtimes,A,vec1));
end
function loopIt
    % Just a regular loop through the pages
    for k = size(A,3):-1:1
    tmp = A(:,:,k);
    sumk(k) = sum(tmp(:));
end
end

end

Here is some sample results:

>> t = timeGpu(100,1000)
t = 
                 t    
             _________

    Reshape    0.0010045
    pageIt     0.0079246
    loopIt       0.27821
>> t = timeGpu(100,100)
t = 
                   t     
               __________

    Reshape    0.00035554
    pageIt      0.0014899
    loopIt       0.028113
>> t = timeGpu(1000,100)
t = 
                  t    
               ________

    Reshape    0.009075
    pageIt     0.052044
    loopIt      0.03972

It looks like reshape wins.

[Joss Knight]

The reshape method is the recommended method. reshape doesn't actually create a new array it just modifies header information. pagefun(@mtimes, ...) is too blunt an instrument.

[Christopher Desmond]

Thanks, really excellent work.

I went with the reshape method as reshape doesn't actually change anything, only how you interact with the matrix.

All in all, I had the series of images (71x71x5000), instead of reading them in that way, as I read them in I reshaped them into row vectors (5000x5041), and then I was easily able to run sum(images,2) to sum every row.

It's actually quite impressive how much time you can save by simply changing the representation of the matrix.

[Christopher Desmond]

If you have time, would you mind taking a look at another question I have. Its the only other problem I've encountered working with the GPU.

http://www.mathworks.com/matlabcentral/answers/293410-help-with-efficient-collection-of-sub-matrices-on-gpu

Thanks for all of your help.

[Brendan Hamm]

Unfortunately, I will not have access to my GPU machine for a few weeks here. I'd be more than happy to make some suggestions on your other post after this though. Likely someone else will get back to you before then.