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循环的情况下实现这个目标，以便更快地计算结果？

- ftxx

---

你想要快速的结果还是避免使用 for 循环？对于这个问题，在现代版本的MATLAB中，我强烈怀疑使用 for 循环是最快的解决方案。如果你有性能问题，我怀疑你正在根据过时的建议寻找解决方案而不是正确的地方。 - Will

@Will 我希望结果能够快速得出。我的矩阵非常大。使用for循环在我的电脑上运行时间约为30秒。我想知道是否有一些聪明的向量化操作，或者mapReduce等方法可以提高速度。 - KevinKim

我怀疑你做不到。向量化仅适用于精确描述的向量和矩阵，但是你的结果允许长度不同的向量。因此，我的假设是，你总会有一些显式循环或类似“cellfun”的伪装循环。 - HansHirse

@ftxx 多大？一行通常有多少个 1？我不会期望一个 find 循环需要花费接近 30 秒的时间，对于任何小到足以适应物理内存的东西。 - Will

@ftxx 请查看我的更新答案，我已经进行了编辑，以实现轻微的性能改进。 - Wolfie

Answer 1

在这个答案的底部有一些基准测试代码，因为您明确表示您对性能感兴趣而不是任意避免for循环。

实际上，我认为for循环可能是最有效的选择。自从引入了“新”的（2015b）JIT引擎以来（source），for循环并不是本质上缓慢的--事实上它们在内部进行了优化。

您可以从基准测试中看到，由ThomasIsCoding here提供的mat2cell选项非常缓慢...

如果我们删除那行代码以使比例更清晰，那么我的splitapply方法就会变得相当慢，obchardon的accumarray选项稍微好一些，但最快（且可比较）的选项要么是使用arrayfun（如Thomas所建议的），要么是使用for循环。请注意，对于大多数用例，arrayfun基本上是一个伪装成for循环的函数，因此这不是一个令人惊讶的平局！

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

编辑：

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

具体而言：

• 使 M 逻辑。如下图所示，对于相对较小的 M，这可能更快，但与类型转换相比，对于大型 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

- Wolfie

---

已经看到并点赞了。 :-) 仍在等待Luis; 他肯定有一些黑色的MATLAB魔法来解决这个问题。 - HansHirse

@Hans 哈哈，是的，尽管他通常使用的技巧（隐式扩展、巧妙的索引等）通常会将事物保持为矩阵，但这里的瓶颈在于对单元格进行汇总。 - Wolfie

请注意，这些时间强烈依赖于 M 的稀疏程度。例如，如果只有5％的元素被填充为 M = randi([0,20],N) == 20;，那么 for 循环远远比您的 arrayfun 方法慢。 - Will

@HansHirse :-) 我的方法是使用 accumarray 而不是 ind2sub，但它比 for 循环慢。 - Luis Mendo

Answer 2

您可以尝试使用以下方式的arrayfun，它会遍历M的行。

A = arrayfun(@(r) find(M(r,:)),1:size(M,1),'UniformOutput',false)

A =
{
  [1,1] =  2
  [1,2] =

     1   3

  [1,3] =

     1   2   3

}

或者（通过mat2cell的较慢方法）

[i,j] = find(M.');
A = mat2cell(i,arrayfun(@(r) sum(j==r),min(j):max(j)))

A =
{
  [1,1] =  2
  [2,1] =

     1
     3

  [3,1] =

     1
     2
     3

}

- ThomasIsCoding

---

虽然arrayfun基本上是一个伪装成循环的循环，因此这可能会在避免循环和像OP所希望的那样快速两方面都失败。 - Wolfie

Answer 3

使用 accumarray：

M = [0 1 0
     1 0 1
     1 1 1];

[val,ind] = find(M.');

A = accumarray(ind,val,[],@(x) {x});

- obchardon

---

在Octave和MATLAB Online中的执行时间大约是一个简单的for循环的两倍，例如：MM{I} = find(M(I, :))。 - HansHirse

@Hans 你可能想看一下我的回答。 - Wolfie

是的，由于每个单元格的大小不同，因此这个问题无法完全向量化（或者有一个我没有看到的技巧）。这只是一个隐藏了for循环的解决方案。 - obchardon

不需要使用ind2sub函数：[ii, jj] = find(M); accumarray(ii, jj, [], @(x){x}) - Luis Mendo

Answer 4

你可以使用strfind函数：strfind。

A = strfind(cellstr(char(M)), char(1));

- rahnema1

---

我（懒惰地）甚至没有查看文档，但是使用实际的“字符串”类型会更快吗？有很多针对字符串的优化，这就是它们存在的原因... - Wolfie

@Wolfie 我认为数值数组与字符数组更相似，因此将数值数组转换为字符数组应该比转换为字符串更直接。 - rahnema1

Answer 5

编辑：我添加了一个基准测试，结果显示for循环比accumarray更有效率。

---

你可以使用find和accumarray:

[c, r] = find(A');
C = accumarray(r, c, [], @(v) {v'});

矩阵被转置（A'）是因为find按列分组。

示例：

A = [1 0 0 1 0
     0 1 0 0 0
     0 0 1 1 0
     1 0 1 0 1];

%  Find nonzero rows and colums
[c, r] = find(A');

%  Group row indices for each columns
C = accumarray(r, c, [], @(v) {v'});

% Display cell array contents
celldisp(C)

输出：

C{1} = 
     1     4

C{2} = 
     2

C{3} =
     3     4

C{4} = 
     1     3     5

基准测试：

m = 10000;
n = 10000;

A = randi([0 1], m,n);

disp('accumarray:')
tic
[c, r] = find(A');
C = accumarray(r, c, [], @(v) {v'});
toc
disp(' ')

disp('For loop:')
tic
C = cell([size(A,1) 1]);
for i = 1:size(A,1)
    C{i} = find(A(i,:));
end
toc

结果：

accumarray:
Elapsed time is 2.407773 seconds.

For loop:
Elapsed time is 1.671387 seconds.

使用for循环比使用accumarray更高效...

- Eliahu Aaron

---

这基本上就是obchardon已经提出的方法，对吧？ - Wolfie

是的，我有点慢，我在发布我的答案后才看到他的回答。 - Eliahu Aaron