Matlab 向量化 - 单元格的非零矩阵行索引

Question

我正在使用 Matlab。

我有一个二元方阵。对于每一行，有一个或多个 1 的条目。我想遍历此矩阵的每一行并返回这些 1 的索引并将它们存储在单元格的条目中。

我想知道是否有一种方法可以在不循环此矩阵的所有行的情况下执行此操作，因为在 Matlab 中 for 循环非常慢。

例如，我的矩阵

M = 0 1 0
    1 0 1
    1 1 1 

然后最终，我想要类似的东西

A = [2]
    [1,3]
    [1,2,3]

A细胞也是如此。

有没有办法在不使用 for 循环的情况下实现这个目标，目的是更快地计算结果？

Answer 1

这个答案的底部是一些基准测试代码，因为您澄清说您对性能感兴趣而不是随意避免for循环。

事实上，我认为for循环可能是这里性能最高的选项。由于引入了“新”（2015b）JIT 引擎（源代码），for循环本身并不慢——事实上它们在内部进行了优化。

您可以从基准测试中看到mat2cellThomasIsCoding此处提供的选项非常慢...

如果我们去掉那条线以使比例更清晰，那么我的splitapply方法相当慢，obchardon 的accumarray 选项好一点，但最快（和可比）的选项要么使用arrayfun（也由 Thomas 建议）或for循环。请注意，对于大多数用例来说，这arrayfun基本上是一个for伪装的循环，所以这并不奇怪！

我建议您使用for循环来提高代码可读性和最佳性能。

编辑：

如果我们假设循环是最快的方法，我们可以围绕find命令进行一些优化。

具体来说

• 让M逻辑。如下图所示，这对于相对 small 来说可能更快M，但在对 large 进行类型转换的权衡时会更慢M。

• 使用逻辑M来索引数组1:size(M,2)而不是使用find. 这避免了循环中最慢的部分（find命令）并超过了类型转换的开销，使其成为最快的选择。

以下是我对最佳性能的建议：

function A = f_forlooplogicalindexing( M )
    M = logical(M);
    k = 1:size(M,2);
    N = size(M,1);
    A = cell(N,1);
    for r = 1:N
        A{r} = k(M(r,:));
    end
end

我已将此添加到下面的基准测试中，这是循环式方法的比较：

基准代码：

rng(904); % Gives OP example for randi([0,1],3)
p = 2:12; 
T = NaN( numel(p), 7 );
for ii = p
    N = 2^ii;
    M = randi([0,1],N);

    fprintf( 'N = 2^%.0f = %.0f\n', log2(N), N );

    f1 = @()f_arrayfun( M );
    f2 = @()f_mat2cell( M );
    f3 = @()f_accumarray( M );
    f4 = @()f_splitapply( M );
    f5 = @()f_forloop( M );
    f6 = @()f_forlooplogical( M );
    f7 = @()f_forlooplogicalindexing( M );

    T(ii, 1) = timeit( f1 ); 
    T(ii, 2) = timeit( f2 ); 
    T(ii, 3) = timeit( f3 ); 
    T(ii, 4) = timeit( f4 );  
    T(ii, 5) = timeit( f5 );
    T(ii, 6) = timeit( f6 );
    T(ii, 7) = timeit( f7 );
end

plot( (2.^p).', T(2:end,:) );
legend( {'arrayfun','mat2cell','accumarray','splitapply','for loop',...
         'for loop logical', 'for loop logical + indexing'} );
grid on;
xlabel( 'N, where M = random N*N matrix of 1 or 0' );
ylabel( 'Execution time (s)' );

disp( 'Done' );

function A = f_arrayfun( M )
    A = arrayfun(@(r) find(M(r,:)),1:size(M,1),'UniformOutput',false);
end
function A = f_mat2cell( M )
    [i,j] = find(M.');
    A = mat2cell(i,arrayfun(@(r) sum(j==r),min(j):max(j)));
end
function A = f_accumarray( M )
    [val,ind] = ind2sub(size(M),find(M.'));
    A = accumarray(ind,val,[],@(x) {x});
end
function A = f_splitapply( M )
    [r,c] = find(M);
    A = splitapply( @(x) {x}, c, r );
end
function A = f_forloop( M )
    N = size(M,1);
    A = cell(N,1);
    for r = 1:N
        A{r} = find(M(r,:));
    end
end
function A = f_forlooplogical( M )
    M = logical(M);
    N = size(M,1);
    A = cell(N,1);
    for r = 1:N
        A{r} = find(M(r,:));
    end
end
function A = f_forlooplogicalindexing( M )
    M = logical(M);
    k = 1:size(M,2);
    N = size(M,1);
    A = cell(N,1);
    for r = 1:N
        A{r} = k(M(r,:));
    end
end