On shape-agnostic slicing of ndarrays

13
votes

In this post, I'm using the term slice to refer to a subarray B_i of an n-dimensional array A such that size(B_i, d) is 1, for some dimension d. A consists of size(A, d) such slices, concatenated along dimension d.

For example, if ndims(A) is 6 and d is 3, then the expressions of the form 

A(:, :, i, :, :, :)

for i in 1:size(A, d) represent all the slices (along dimension d) that make up A.

The problem with an expression like A(:, :, i, :, :, :) is that it cannot be generalized symbolically to slices along a dimension different from 3 in arrays having a number of dimensions different from 6. E.g., to get A's slices along dimension 2, one would need a different expression, A(:, i, :, :, :, :). This means that such expressions are useless in code that is agnostic about the shape of some array from which slices are to be extracted.

The function below is my matlab-noob attempt to implement shape-agnostic slicing. (The name slice is already taken, hence I called the function hslice, short for hyperslice.) The function's strategy is to reshape the input array into a suitable 3-d array, take the desired slice along the reshaped array's second dimension, and reshape the result to have the shape of a slice from the original input array.

function out = hslice(ndarray, d, i)
    sz = size(ndarray);
    pfx = sz(1:d-1);    % dimensions before d
    sfx = sz(d+1:end);  % dimensions after d

    tmp = reshape(ndarray, prod(pfx), sz(d), prod(sfx));
    out = reshape(tmp(:, i, :), [pfx 1 sfx]);
end

Is there a built-in, or at least a more efficient, way to achieve the same result (in a shape-agnostic way)?

2
matlabmultidimensional-array
why do you think your method is not efficient? have you tried to compare it with explicit slicing? – Shai
I actually like yours because no data is rewritten prior to slicing. – chappjc
@Shai: I was mostly hoping for some low-level, built-in (and presumably very efficient) function for this; also since I don't know much about what's efficient in MATLAB, my assumption is that any implementation by me is likely to be inefficient (though maybe I got lucky this time). IMO, profiling (and not just in MATLAB) is tricky. In this case any profiling scheme I can think of (as MATLAB noob) would require putting the explicit slicing expression inside a function, which is likely to distort the comparison significantly. IOW, I wouldn't trust any profiling I could do on this. – kjo
Yeah, Matlab performance is hard to get a handle on because the cost of operations can vary widely, so O(n) analysis is often less useful than figuring out how to "vectorize" stuff with fast Matlab built-ins. But all these indexing, reshaping, and copy operations you're using are quite cheap, relatively speaking. And doing it in a function is fine: Because you're making a new array, there's no in-place optimization possible, and copy-on-write means passing arrays between functions is cheap. Your method should perform fine. – Andrew Janke
And oh yeah - reshape() in particular is uber-fast, because it doesn't actually rearrange your raw primitive data in memory: it just twiddles the "dimensions" info in the header of Matlab's internal data structure (the "mxArray"), and retains a pointer to the original underlying primitive data. So your implementation here is only making a single extraction/copy of the primitive data, just like the static A(:, :, i, :, ...) indexing expression. – Andrew Janke

2 Answers

17
votes

Yeah. You can use the equivalence between dereferenced cell arrays and "comma-separated lists", and the fact that you can use a char ':' as an index to dynamically construct that A(:, :, :, i, :, ...) invocation for arbitrary dimensions.

function out = slice(A, ix, dim)

subses = repmat({':'}, [1 ndims(A)]);
subses{dim} = ix;
out = A(subses{:});

This will generalize fully, and will do exactly the same "slice" indexing operation as the original static A(:, :, i, :, ...) expression, aside from the overhead of twiddling those strings to set it up.

Or if you wanted to, you could just use sprintf to construct that A(:, :, i, :, ...) as a string and then call eval() on it. But I like to avoid eval if at all possible.

Note that your original implementation is using fast operations and should perform just fine, about as fast as this one. I'm just posting this because I think it's very readable, does answer your question as originally posed, and it could be applied to other useful stuff.

Assignment to Slices

You can also use this same subscripts-in-a-cell technique as an lvalue to assign in to slices of an array. You can't reuse the slice function directly, though, since it returns an extracted subset of the array, and not an lvalue reference. So you can make a very similar function that does the assignment itself.

function A = slice_assign(A, ix, dim, B)
%SLICE_ASSIGN Assign new values to a "slice" of A

subses = repmat({':'}, [1 ndims(A)]);
subses{dim} = ix;
A(subses{:}) = B;

In practice, you might also want a function that just returns the computed indexes in a cell array, so you could carry those around and use them repeatedly for assignment and referencing.

function out = slice_subs(A, ix, dim)

subses = repmat({':'}, [1 ndims(A)]);
subses{dim} = ix;
out = subses;
3
votes

You can try permute and setdiff to move that dimension to a consistent position:

function out = hslice(ndarray, d, i)
subdims = setdiff(1:ndims(ndarray),d);
sz = size(ndarray);
outsz = sz(subdims);
order = [d subdims];
ndarray = permute(ndarray,order);
out = reshape(ndarray(i,:),outsz);
end

For example:

d = 3; i = 2;
nd = randi(23,3,3,3,2);
out = hslice(nd,d,i);   % out = squeeze(nd(:,:,i,:)) for d=3

BUT, data is rewritten prior to slicing here, and not with the code in the question. So, I'd actually go with the OP's!