使用向量索引没有线性索引的矩阵

Question

G'day,我试图找到一种方法来使用[x,y]点的矢量来从MATLAB中的大矩阵中索引.通常,我会将下标点转换为矩阵的线性索引.(例如,使用向量作为矩阵的索引)但是,矩阵是4维的,我想要采用的所有元素第3和第4维度具有相同的第1和第2维度.让我希望通过一个例子来证明:

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)等,并返回不仅是999和888存在于Matrix(:,:,1,1),但其中存在的内容Matrix(:,:,1,2),Matrix(:,:,2,1)并且Matrix(:,:,2,2),等等(IRL ,尺寸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

我还想补充一点,由于我正在处理的矩阵的大小,速度是一个问题 - 因此我想使用下标索引来索引数据.

我还应该提一下(不像这个问题:使用下标而不使用sub2ind访问值)因为我想要存储在i和jth索引的额外维度3和4中的所有信息,我不认为稍快一点版sub2ind仍然不会削减它..

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,使用简单循环的索引最简单,最简单,也很快.