This can be done on the GPU I think, but you will definitely need to get rid of your global variables. This doc page should help here - you need to parameterise your backProjectionKernel function to avoid it needing to use global.
The main trick here is to split the loop into two pieces - firstly, an arrayfun portion that performs the independent calls to backProjectionKernel. This relies on the implicit dimension expansion to "loop" over two dimensions. Then, the second piece uses accumarray to build up the result.
szI = 5;
szJ = 10;
iVec = gpuArray(1:szI);
jVec = gpuArray(1:szJ)';
% First, perform a dummy computation that calculates the indices of J
% to accumulate into, along with the indices of I from which to take
% the values. Note that the gpuArray version of arrayfun performs
% implicit dimension expansion, so this is effectively a double-loop
% over all combinations of iVec and jVec
[indJ, indI] = arrayfun(@iDummyCalc, iVec, jVec);
% Build dummy I
I = rand(szI, 'gpuArray');
% Compute J using accumarray - first as a vector...
J = accumarray(indJ(:), I(indI(:)), [szJ * szJ, 1]);
% ... then reshape to the matrix
J = reshape(J, [szJ, szJ]);
function [indJ, indI] = iDummyCalc(iVal, jVal)
% Compute some dummy values that are valid linear indices into
% I and J.
indJ = randi([jVal 100]);
indI = randi([iVal 25]);
end
