TL;DR - 这是我的建议：
nI = numel(idx);
sz = size(m);
x = sparse( 1:nI, idx, m(sub2ind( size(m), 1:numel(idx), idx )), sz(1), sz(2), nI);
帖子的其余部分讨论了为什么它会更好。
看到所需的输出矩阵主要由零组成，这实际上要求使用sparse matrices！这不仅应该提高性能（尤其是对于较大的矩阵），而且应该对内存更加友好。
我要给Wolfie's benchmark添加两个函数：
function x = f_sp_loop( m, idx )
nI = numel(idx);
sz = size(m);
x = spalloc( sz(1), sz(2), nI ); % Initialize a sparse array.
for indI = 1:nI
x( indI, idx(indI) ) = m( indI, idx(indI) ); % This generates a warning (inefficient)
end
end
function x = f_sp_sub2ind( m, idx )
nI = numel(idx);
sz = size(m);
x = sparse( 1:nI, idx, m(sub2ind( size(m), 1:numel(idx), idx )), sz(1), sz(2), nI);
end
本质上的区别在于，我们没有将输出预分配为一个零数组，而是作为一个稀疏数组。基准1的结果如下：
...这提出了一个问题 - sparse 方法为什么要快一个数量级？
要回答这个问题，我们应该查看基准函数内部的实际运行时分布，我们可以从profiler 中获得。要获得更多信息，我们可以make the profiler output memory consumption info，使用profile('-memory','on')。在运行一个较短版本的基准2 之后，它只运行k 的最高值，我们得到：
所以我们可以总结几件事：
- 运行时的绝大部分时间都花在分配和释放内存上，这就是算法看起来几乎具有相同性能的原因。因此，如果我们减少内存分配，我们会直接节省时间（
sparse!）。
- 即使
sub2ind 和loop 方法看起来相同，我们仍然可以在两者之间建立一个“赢家”（见下图紫色框）-sub2ind！ sub2ind 是 32 毫秒，而循环是 41 毫秒。
-
稀疏循环方法的速度并不令人意外，正如mlint 警告我们的那样：
说明
代码分析器检测到可能很慢的稀疏数组的索引模式。更改稀疏数组的非零模式的赋值可能会导致此错误，因为此类赋值会导致相当大的开销。
建议的操作
如果可能，使用sparse 构建稀疏数组，如下所示，并且不要使用索引赋值（例如C(4) = B）来构建它们：
- 创建单独的索引和值数组。
- 调用 sparse 来组装索引和值数组。
如果您必须使用索引分配来构建稀疏数组，您可以通过首先使用 spalloc 预分配稀疏数组来优化性能。
如果代码只更改已经非零的数组元素，那么开销是合理的。按照Adjust Code Analyzer Message Indicators and Messages 中的说明隐藏此消息。
有关详细信息，请参阅“Constructing Sparse Matrices”。
- 结合了两全其美的方法，即
sparse 的内存节省和sub2ind 的矢量化，似乎是运行时间仅为 3 毫秒的最佳方法！
1 图表制作代码：
function q51605093()
K = 15; % Max loop variable
T = zeros( K, 4 ); % Timing results
for k = 1:K
M = 2^(k-1); N = 2^k; % size of matrix
m = rand( M, N ); % random matrix
idx = randperm( N, M ); % column indices
% Define anonymous functions with no inputs, for timeit, and run
f = cell(4,1);
f{1} = @() f_sub2ind( m, idx );
f{2} = @() f_loop( m, idx );
f{3} = @() f_sp_loop( m, idx );
f{4} = @() f_sp_sub2ind( m, idx );
T(k,:) = cellfun(@timeit, f);
if k == 5 % test equality during one of the runs
R = cellfun(@feval, f, 'UniformOutput', false);
assert(isequal(R{:}));
end
end
% Plot results
figure();
semilogy( (1:K).', T, 'linewidth', 2 ); grid on; xticks(0:K);
legend( {'sub2ind', 'loop', 'sp\_loop', 'sp\_sub2ind'}, 'Location', 'NorthWest' );
xlabel( 'k, where matrix had 2^{(k-1)} rows and 2^k columns' );
ylabel( 'function time (s)' )
end
function x = f_sub2ind( m, idx )
% Using the in-built sub2ind to generate linear indices, then indexing
lin_idx = sub2ind( size(m), 1:numel(idx), idx );
x = zeros( size(m) );
x( lin_idx ) = m( lin_idx );
end
function x = f_loop( m, idx )
% Directly indexing in a simple loop
x = zeros( size(m) );
for ii = 1:numel(idx)
x( ii, idx(ii) ) = m( ii, idx(ii) );
end
end
function x = f_sp_loop( m, idx )
nI = numel(idx);
sz = size(m);
x = spalloc( sz(1), sz(2), nI ); % Initialize a sparse array.
for indI = 1:nI
x( indI, idx(indI) ) = m( indI, idx(indI) ); % This generates a warning (inefficient)
end
end
function x = f_sp_sub2ind( m, idx )
nI = numel(idx);
sz = size(m);
x = sparse( 1:nI, idx, m(sub2ind( size(m), 1:numel(idx), idx )), sz(1), sz(2), nI);
end
2 分析代码：
function q51605093_MB()
K = 15; % Max loop variable
M = 2^(K-1); N = 2^K; % size of matrix
m = rand( M, N ); % random matrix
idx = randperm( N, M ); % column indices
% Define anonymous functions with no inputs, for timeit, and run
f = cell(4,1);
f{1} = f_sub2ind( m, idx );
f{2} = f_loop( m, idx );
f{3} = f_sp_loop( m, idx );
f{4} = f_sp_sub2ind( m, idx );
% assert(isequal(f{:}));
end
... the rest is the same as above