MATLAB中每一行出现的索引

Question

我有两个矩阵，A 和 B。 （B 和1:n 一样是连续的）

我需要在A 中找到B 的每一行的所有出现，并将这些行索引相应地存储在单元格数组C 中。请参阅下面的示例。

A = [3,4,5;1,3,5;1,4,3;4,2,1]
B = [1;2;3;4;5]

因此，

C = {[2,3,4];[4];[1,2,3];[1,3,4];[1,2]}

注意C 不需要在我的应用程序的元胞数组中。我只建议它，因为C 的行向量长度不等。如果您可以提出解决方法，这也很好。

我尝试对B 的每一行使用运行ismember 的循环，但是当矩阵A 和B 很大时，这太慢了，大约有一百万个条目。向量化代码表示赞赏。

（为了给你上下文，这样做的目的是在网格中识别那些附加到单个顶点的面。注意我不能使用函数 edgeattachments 因为我的数据不是“TR”形式三角剖分表示。我只有一个面列表和顶点列表。）

Answer 1

我们可以在不对B 做任何假设的情况下做到这一点。试试bsxfun 和mat2cell 的这种用法：

M = squeeze(any(bsxfun(@eq,A,permute(B,[3 2 1])),2)); % 4x3x1 @eq 1x1x5 => 4x3x5
R = sum(M); % 4x5 -> 1x5
[ii,jj] = find(M);
C = mat2cell(ii,R)

上面C 中的单元格将是列向量，而不是您的示例中的行。要使单元格包含行向量，请改用C = mat2cell(ii',1,R)'。

我唯一担心的是 mat2cell 对于数百万个 R 的值可能会很慢，但是如果您希望在单元格中输出，我不确定您能做得更好。 编辑：如果您可以像在 Werner 的第一个解决方案中那样使用循环处理稀疏矩阵，请将上面的最后一行替换为以下内容：

>> Cs = sparse(ii,jj,1)
Cs =
   (2,1)        1
   (3,1)        1
   (4,1)        1
   (4,2)        1
   (1,3)        1
   (2,3)        1
   (3,3)        1
   (1,4)        1
   (3,4)        1
   (4,4)        1
   (1,5)        1
   (2,5)        1

不幸的是，如果 size(A,1) 和 numel(B) 都很大，bsxfun 可能会耗尽内存！如果内存成为问题，您可能必须遍历 A 或 B 的元素。这是通过循环遍历B 中的顶点的一种方法：

for i=1:numel(B), C{i} = find(any(A==B(i),2)); end

是的，就这么简单。在 MATLAB 中，元胞数组的增长速度非常快，因为它类似于存储对数据的连续引用的序列容器，而不是保持数据本身是连续的。也许ismember 是您测试的瓶颈。

Answer 2

嗯，这个问题的最佳答案需要了解如何填充 A。如果 A 是稀疏的，也就是说，如果它的列值很少并且 B 很大，那么我认为节省内存的最佳方法可能是使用稀疏矩阵而不是单元格。

% No fancy stuff, just fast and furious 
bMax = numel(B);
nRows = size(A,1);

cLogical = sparse(nRows,bMax);

for curRow = 1:nRows
  curIdx = A(curRow,:);
  cLogical(curRow,curIdx) = 1;
end

答案：

cLogical =

   (2,1)        1
   (3,1)        1
   (4,1)        1
   (4,2)        1
   (1,3)        1
   (2,3)        1
   (3,3)        1
   (1,4)        1
   (3,4)        1
   (4,4)        1
   (1,5)        1
   (2,5)        1

如何阅读答案。对于每一列，行显示列索引出现在A中的索引。即1出现在行[2 3 4]，2出现在行[4]，3行[1 2 3]，4行[1 3 4]，5 在行中[1 2]。

那么您以后可以根据需要使用cLogical 代替单元格作为索引矩阵。

另一种方法是为 C 分配索引应在 C 中出现多少次的预期值。

% Fancier solution using some assumed knowledge of A
bMax = numel(B);
nRows = size(A,1);
nColumns = size(A,2);

% Pre-allocating with the expected value, an attempt to reduce re-allocations.
% tic; for rep=1:10000; C = mat2cell(zeros(bMax,nColumns),ones(1,bMax),nColumns); end; toc 
% Elapsed time is 1.364558 seconds.
% tic; for rep=1:10000; C = repmat({zeros(1,nColumns)},bMax,1); end; toc
% Elapsed time is 0.606266 seconds.
% So we keep the not fancy repmat solution
C = repmat({zeros(1,nColumns)},bMax,1);
for curRow = 1:nRows
  curIdxMsk = A(curRow,:);
  for curCol = 1:nColumns
    curIdx = curIdxMsk(curCol);
    fillIdx = ~C{curIdx};
    if any(fillIdx) 
      fillIdx = find(fillIdx,1);
    else
      fillIdx = numel(fillIdx)+1;
    end
    C{curIdx}(fillIdx) = curRow;
  end
end

% Squeeze empty indexes:
for curRow = 1:bMax
  C{curRow}(~C{curRow}) = [];
end

答案：

>> C{:}

ans =

     2     3     4


ans =

     4


ans =

     1     2     3


ans =

     1     3     4


ans =

     1     2

哪种解决方案效果最好？您在代码中进行性能测试，因为它取决于 A、bMax 的大小、计算机的内存大小等等。然而，我仍然对其他人可以为这个 x) 做的解决方案感到好奇。我喜欢 chappjc 的解决方案，尽管它有他指出的缺点。

对于给定的示例（10k 次）：

Solution 1: Elapsed time is 0.516647 seconds. 
Solution 2: Elapsed time is 4.201409 seconds (seems that solution 2 is a bad idea hahaha, but since it was created to the specific issue of A having many rows it has to be tested in those conditions).
chappjc' solution: Elapsed time is 2.405341 seconds.