Separable Block-wise Operations

Efficiently performs separable operations (e.g., sum, mean,prod, min, max,...) on array sub-blocks

This contribution was inspired by a series of posts by Bruno Luong and Jan Simon. Often, we wish to divide an array into equal sized sub-blocks and perform an operation on these blocks, reducing them to scalars. For general block operations, MATLAB has made available functions like BLOCKPROC. It is often also possible to do this kind of processing by using MAT2CELL to subdivide the array into cells, each containing a block, and then applying CELLFUN.
For certain specific and common kinds of operations, however, it is possible to do the computation in a particularly efficient manner. These functions include sum(), prod(), mean(), max(), min(),... which can be done separably along each dimension of the block, first along columns, then along rows, etc... For example,

sum(B)=sum(sum(sum(B,1),2),3)

By decomposing these functions into separable calls, it is possible to do them blockwise in an array with a minimum of data copying and a high degree of vectorization and sequential memory access. The mfunction SEPBLOCKFUN in this submission optimizes the computation in this way for any separable function the user cares to supply.

USAGE:

Y=sepblockfun(X,blockdims,fun)

IN:


X: A full array. If the ndSparse class defintion is on the path, then X
can also be a regular sparse matrix or ndSparse array. Performance might
not be as strong as for full arrays, however.

blockdims: a vector of integers specifying the dimensions of the
sub-blocks. The array X must partition evenly into blocks
of this size. If blockdims(i) is set to Inf then it will be
replaced with blockdims(i)=size(X,i).

fun: function handle to an operation assumed to be separable
(Examples: max,min,sum,prod,mean, etc...). The function must
accept the input syntax fun(B,DIM) where B is an input array
and DIM is a dimension along which to operate. Alternatively,
fun can be one of the following strings 'max','min','sum','mean',
'prod'.


OUT:

Y: the output array. Y(i)=fun(Xi(:),1) where Xi is the i-th sub-block of
the input array X.


EXAMPLE 1:

Divide a 400x400x400 array into 10x10x10 blocks. Return the blockwise,
mean, max., and min. of each block, each organized as a 40x40x40 array.

A=rand(400,400,400);
Ameans=sepblockfun(A,[10,10,10],@mean);
Amins=sepblockfun(A,[10,10,10],'min' );
Amaxs=sepblockfun(A,[10,10,10], @(B,d) max(B,[],d) );

EXAMPLE 2:

Not all operations satisfy the separability property, but
sometimes inseparable operations can be decomposed into separable ones. As
an example, we take the blockwise standard deviations of the same array
from Example 1.

Astds=sqrt( sepblockfun(A.^2,[10,10,10],'mean') - Ameans.^2 );

It is also possible to use SEPBLOCKFUN for sparse matrices and arrays if you happen to have my ndSparse class <http://www.mathworks.com/matlabcentral/fileexchange/29832-n-dimensional-sparse-arrays> installed. However, sparse arrays may benefit less from the computation strategy than full arrays, in part because RESHAPE operations are less efficient for sparse matrices. Also, it is only possible to apply separable functions that have been overloaded for the ndSparse class.