使用向量索引矩阵而不是线性索引

Question

你好， 我正在尝试找到一种方法，使用[x，y]点的向量在MATLAB中从大矩阵中进行索引。 通常，我会将下标点转换为矩阵的线性索引。（例如：使用向量作为矩阵索引）但是，该矩阵是4维的，并且我想获取所有具有相同第1和第2维的第3和第4维元素。让我用一个例子来演示：

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;

现在，我想要使用下标（1,2）和（3,4）等进行索引，并且不仅返回Matrix(:,:,1,1)中存在的999和888，还要返回Matrix(:,:,1,2)，Matrix(:,:,2,1)和Matrix(:,:,2,2)等位置上存在的内容，以此类推（实际情况中，Matrix的维度可能更像size(Matrix) = (300 250 30 200)）。我不想使用线性索引，因为我希望结果以类似向量的方式呈现。例如，我希望得到一个类似以下的结果：

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

我还想补充一点，由于我处理的矩阵很大，速度是一个问题 - 这就是为什么我想使用下标索引来索引数据的原因。

我还应该提到（与此问题不同：Accessing values using subscripts without using sub2ind），由于我希望所有信息都存储在第i和j个下标的额外维度3和4中，我认为稍微快一点的sub2ind版本仍无法解决问题。

- David_G

Answer 1

我可以想到三种方法来处理这个问题。

简单循环

只需循环所有的2D索引，并使用冒号访问其余维度：

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

向量化计算线性索引

不使用sub2ind，向量化计算线性索引：

% 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

您可以使用Matlab自带的现成工具：

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

速度

那么哪一个最快呢？让我们来一探究竟：

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

结果：

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

结论：使用循环。

尽管上述测试在很大程度上受到JIT未能编译最后两个循环以及对4D数组的非特定性影响（最后两种方法也适用于ND数组）。为4D制作专门版本无疑会更快。

然而，使用简单循环进行索引是最简单、最容易理解和非常快的，这要归功于JIT。

- Rody Oldenhuis

---

非常详细的回答，我会进行一些测试。谢谢！ - David_G

非常好的答案 - 尽管我害怕for循环，但第一种方法似乎运行良好。此外，我不知道JIT加速器，所以感谢您提到它！ - David_G

Answer 2

所以，这里有一个可能的答案...但它很混乱。我怀疑它比更直接的方法更加计算昂贵...而且这绝对不是我的首选答案。如果我们能够在没有任何for循环的情况下得到答案，那将是很棒的！

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

结果矩阵result的大小应为size(result) = (4 N tt)。

我想这应该能运行，即使Matrix不是正方形。然而，正如我上面所说，这并不理想。

- David_G