What are efficient ways to extract lower dimensional slices from a high-dimensional array ?

I am looking for an efficient way to extract sub-arrays from a larger array. There are (at least) two standard solutions both having certain disadvantages.

First solution: Use a loop Disadvantage: Might be slow if nidx (see code below) is large

Second solution: use sub2ind Disadvantage: Interim indices created are large if nidx * ndim3 is large

Any other (especially better) ideas?

Kind regards, gg

   %create 3dim testarray
   ndim1 = 3;
   ndim2 = 4;
   ndim3 = 5;
   testarr = rand([ndim1,ndim2,ndim3]);

   %pick some indices into the first two dimensions
   nidx = 10;
   idx1 = randi([1,ndim1],[nidx,1]);
   idx2 = randi([1,ndim2],[nidx,1]);

   %now extract slices from testarr according to the indices

   %First solution: by loop
   newarr = zeros(nidx, ndim3);
   for kidx = 1:nidx
      newarr(kidx,:) = testarr(idx1(kidx), idx2(kidx), :);
   end

   %Second solution using sub2ind
   refidx1 = repmat(idx1, [ndim3,1]);
   refidx2 = repmat(idx2, [ndim3,1]);
   tmp = repmat((1:ndim3),[nidx, 1]);
   refidx3 = tmp(:);
   lidx = sub2ind(size(testarr), refidx1, refidx2, refidx3);
   newarr_alt = reshape(testarr(lidx), nidx, ndim3);

   %are they equal this time?
   isequal(newarr, newarr_alt)

[James Tursa]

What are you doing with these slices downstream? There may be ways to process them without physically copying them to a new variable.

[gg]

James, assume I would like to square each element of a slice, or set them to zero. What would this change? The problem is accessing (addressing) those elements efficiently. Having solved that I could something "in situ" on them or copy into a new variable.

In addition I would like the solution be "generic" in a sense so that I do not need to fiddle with index arrays everytime this problem occurs.

[James Tursa]

Copying the data can easily dominate the run time over the calculations involved, particularly for something as simple as squaring elements or setting elements to a value. E.g., a mex routine can often significantly outperform m-code for these cases.

[gg]

Now I understand! Yes, of course you are completely right. But then on a MEX level this problem doesn't arise. Andrei's way of doing index arithmetic is the way to do it in C, isn't it?

So I was really asking (not entirely precise obviously) on a MATLAB intrinsic way of doing it, without using MEX.

Note that if you are using Andrei's method you create an index (somewhat implicitely though). So even if you do in situ calculations you may need to allocate memory for this index. Actually it is not clear to me under what conditions, which solution is the "best". Except for MEX maybe which seems to be straightforward and thus quick. I have never used it but then there might be overheads involved as well?

Answer 1

[Andrei Bobrov]

s = size(testarr);
m = prod(s(1:2));
out = testarr(bsxfun(@plus,idx2*(s(1)-1)+idx1,0:m:(s(3)-1)*m))

[gg]

Andrei,

Thanks!

I consider this an optimissed variant of sub2ind. It does suffer from creating an auxiliary index as well. Although the size is only ndim3 as opposed to ndim3 * nidx, which is certainly a huge improvement. It will be interesting to compare speed.