arrayfun with function with inputs of different dimensions

Question

I'm trying to create a matrix that contains the averages of the kxk submatrices of a larger nxn matrix, where n is divisible by k. I can accomplish this fairly efficiently with something like this

mat = mat2cell(mat, k*ones(1,n/k), k*ones(1,n/k))
mat = cellfun(@mean,mat,'UniformOutput',false);
mat = cellfun(@mean,mat,'UniformOutput',false); %repeated to collapse cells to 1x1
mat = cell2mat(mat)

However, since I have a very large amount of data all in very large matrices, repeating this process can still take a long time even on a cluster, and combining the matrices is not an option due to memory limitations. I'm wondering if it's possible to rewrite this code using arrayfun instead so I can utilize its GPU capability (as GPU cannot process cell arrays), but I'm running into problems since the helper function

function avg = blockavg(mat,i,j,k)
i1 = (i-1)*k+1;
i2 = i*k;
j1 = (j-1)*k+1;
j2 = j*k
avg = mean(mean(mat(i1:i2,j1:j2)));
end

calls for two broadcast variables mat, an nxn matrix, and k, a scalar, and these are not array inputs when plugged into arrayfun. When running

ii = 1:(n/k);
jj = 1:(n/k);
matavg = arrayfun(@blockavg,mat,ii,jj,k)

An error message returns stating that input arguments must be arrays of the same size and shape, as only ii and jj are the array inputs, while mat and k do not vary across elements. I'm not quite sure how to work around this issue, or even if this type of operation is possible at all for arrayfun. Any suggestion is appreciated.

Thanks!

Answer 1

Another nice possibility without using the image processing toolbox:

squeeze(mean(mean(reshape(X,k,n/k,k,n/k),1),3))

This code reshapes the matrix in a way, that each block can be accessed by Y(a,:,b,:), which means they have identical first and third dimension index. Then mean is applied to first and third dimension, finally singleton dimensions are removed to get a 2d-result.

Answer 2

What you are trying to do is exactly what blockproc from the image processing toolbox provides:

blockproc(X,[k k],@(x)(mean(mean(x.data))))

To speed up, you may try setting the UseParallel-parameter.