Use a vector to index a matrix without linear index

Question

G'day, I'm trying to find a way to use a vector of [x,y] points to index from a large matrix in MATLAB. Usually, I would convert the subscript points to the linear index of the matrix.(for eg. Use a vector as an index to a matrix) However, the matrix is 4-dimensional, and I want to take all of the elements of the 3rd and 4th dimensions that have the same 1st and 2nd dimension. Let me hopefully demonstrate with an example:

Matrix = nan(4,4,2,2); % where the dimensions are (x,y,depth,time)
Matrix(1,2,:,:) = 999; % note that this value could change in depth (3rd dim) and time (4th time) 
Matrix(3,4,:,:) = 888; % note that this value could change in depth (3rd dim) and time (4th time) 
Matrix(4,4,:,:) = 124;

Now, I want to be able to index with the subscripts (1,2) and (3,4), etc and return not only the 999 and 888 which exist in Matrix(:,:,1,1) but the contents which exist at Matrix(:,:,1,2),Matrix(:,:,2,1) and Matrix(:,:,2,2), and so on (IRL, the dimensions of Matrix might be more like size(Matrix) = (300 250 30 200)

I don't want to use linear indices because I would like the results to be in a similar vector fashion. For example, I would like a result which is something like:

ans(time=1)
999 888 124
999 888 124
ans(time=2)
etc etc etc
etc etc etc

I'd also like to add that due to the size of the matrix I'm dealing with, speed is an issue here - thus why I'd like to use subscript indices to index to the data.

I should also mention that (unlike this question: Accessing values using subscripts without using sub2ind) since I want all the information stored in the extra dimensions, 3 and 4, of the i and jth indices, I don't think that a slightly faster version of sub2ind still would not cut it..

Answer 1

I can think of three ways to go about this

Simple loop

Just loop over all the 2D indices you have, and use colons to access the remaining dimensions:

for jj = 1:size(twoDinds,1)
    M(twoDinds(jj,1),twoDinds(jj,2),:,:) = rand;
end

Vectorized calculation of Linear indices

Skip sub2ind and vectorize the computation of linear indices:

% generalized for arbitrary dimensions of M

sz = size(M);
nd = ndims(M);

arg = arrayfun(@(x)1:x, sz(3:nd), 'UniformOutput', false);

[argout{1:nd-2}] = ndgrid(arg{:});

argout = cellfun(...
    @(x) repmat(x(:), size(twoDinds,1),1), ...
    argout, 'Uniformoutput', false);

twoDinds = kron(twoDinds, ones(prod(sz(3:nd)),1));

% the linear indices
inds = twoDinds(:,1) + ([twoDinds(:,2) [argout{:}]]-1) * cumprod(sz(1:3)).';

Sub2ind

Just use the ready-made tool that ships with Matlab:

inds = sub2ind(size(M), twoDinds(:,1), twoDinds(:,2), argout{:});

Speed

So which one's the fastest? Let's find out:

clc

M = nan(4,4,2,2);

sz = size(M);
nd = ndims(M);

twoDinds = [...
    1 2
    4 3
    3 4
    4 4
    2 1];

tic
for ii = 1:1e3
    for jj = 1:size(twoDinds,1)
        M(twoDinds(jj,1),twoDinds(jj,2),:,:) = rand;
    end
end
toc


tic
twoDinds_prev = twoDinds;
for ii = 1:1e3

    twoDinds = twoDinds_prev;

    arg = arrayfun(@(x)1:x, sz(3:nd), 'UniformOutput', false);

    [argout{1:nd-2}] = ndgrid(arg{:});

    argout = cellfun(...
        @(x) repmat(x(:), size(twoDinds,1),1), ...
        argout, 'Uniformoutput', false);

    twoDinds = kron(twoDinds, ones(prod(sz(3:nd)),1));
    inds = twoDinds(:,1) + ([twoDinds(:,2) [argout{:}]]-1) * cumprod(sz(1:3)).';

    M(inds) = rand;

end
toc


tic
for ii = 1:1e3

  twoDinds = twoDinds_prev;

    arg = arrayfun(@(x)1:x, sz(3:nd), 'UniformOutput', false);

    [argout{1:nd-2}] = ndgrid(arg{:});

    argout = cellfun(...
        @(x) repmat(x(:), size(twoDinds,1),1), ...
        argout, 'Uniformoutput', false);

    twoDinds = kron(twoDinds, ones(prod(sz(3:nd)),1));

    inds = sub2ind(size(M), twoDinds(:,1), twoDinds(:,2), argout{:});

    M(inds) = rand;
end
toc

Results:

Elapsed time is 0.004778 seconds.  % loop
Elapsed time is 0.807236 seconds.  % vectorized linear inds
Elapsed time is 0.839970 seconds.  % linear inds with sub2ind

Conclusion: use the loop.

Granted, the tests above are largely influenced by JIT's failure to compile the two last loops, and the non-specificity to 4D arrays (the last two method also work on ND arrays). Making a specialized version for 4D will undoubtedly be much faster.

Nevertheless, the indexing with simple loop is, well, simplest to do, easiest on the eyes and very fast too, thanks to JIT.

Answer 2

So, here is a possible answer... but it is messy. I suspect it would more computationally expensive then a more direct method... And this would definitely not be my preferred answer. It would be great if we could get the answer without any for loops!

Matrix = rand(100,200,30,400);
grabthese_x = (1 30 50 90);
grabthese_y = (61 9 180 189);
result=nan(size(length(grabthese_x),size(Matrix,3),size(Matrix,4));
for tt = 1:size(Matrix,4)
subset = squeeze(Matrix(grabthese_x,grabthese_y,:,tt));
 for NN=1:size(Matrix,3)
 result(:,NN,tt) = diag(subset(:,:,NN));
 end
end

The resulting matrix, result should have size size(result) = (4 N tt).

I think this should work, even if Matrix isn't square. However, it is not ideal, as I said above.