如何在不使用 for 循环的情况下对不同大小的矩阵的各个部分求和？答案

【问题标题】：How to sum parts of a matrix of different sizes, without using for loops?如何在不使用 for 循环的情况下对不同大小的矩阵的各个部分求和？
【发布时间】：2018-12-06 04:54:41
【问题描述】：

我有一个相对较大的矩阵 NxN (N~20,000) 和一个 Nx1 向量，用于标识必须组合在一起的索引。

我想将矩阵的各个部分相加，原则上可以有不同数量的元素和不相邻的元素。我很快写了一个可以正常工作的双 for 循环，但它当然效率低下。分析器将这些循环识别为我的代码中的瓶颈之一。

我试图找到一种智能矢量化方法来解决这个问题。我探索了arrayfun、cellfun 和bsxfun 函数，并寻找类似问题的解决方案……但我还没有找到最终解决方案。

这是带有两个 for 循环的测试代码：

M=rand(10); % test matrix
idxM=[1 2 2 3 4 4 4 1 4 2]; % each element indicates to which group each row/column of M belongs
nT=size(M,1);
sumM=zeros(max(idxM),max(idxM));
for t1=1:nT
    for t2=1:nT
        sumM(t1,t2) = sum(sum(M(idxM==t1,idxM==t2)));
    end
end

【问题讨论】：

也许值得注意的是 arrayfun 和 cellfun 基本上是变相的循环，很可能优化后的 for 循环也一样快（如果不简洁的话）。
感谢您指出这一点。您在答案中使用的 bsxfun 怎么样？
我假设您的意思是 Luis 的回答...bsxfun 是另一种野兽，通常可以提高速度（取决于它的使用方式）。将cellfun 和arrayfun 视为“我想遍历一个数组并对每个元素应用一个函数”，而bsxfun 是“然后将这些向量扩展为等效大小执行矩阵运算” 所以不是逐个元素。它们有不同的有用应用，bsxfun 适用于数组和矩阵之间的数值运算，而循环函数适用于简洁的代码，而无需初始化输出等。
@user9998992 bsxfun 是 very fast，implicit expansion 也是如此
是的，同意。行：nT=size(M,1); sumM=zeros(max(idxM),max(idxM)); 应与行交换：nT=max(idxM); sumM=zeros(nT,nT);（在我的原始代码中，它们是正确的，但我在测试代码中输入了拼写错误）

标签： matlab matrix vectorization bsxfun

【解决方案1】：

你可以使用accumarray如下：

nT = size(M,1); % or nT = max(idxM)
ind = bsxfun(@plus, idxM(:), (idxM(:).'-1)*nT); % create linear indices for grouping
sumM = accumarray(ind(:), M(:), [nT^2 1]); % compute sum of each group
sumM = reshape(sumM, [nT nT]); % reshape obtain the final result

【讨论】：

从 R2016b 开始，您可以依赖 MATLAB 的隐式扩展，不需要 bsxfun... 第 2 行变为 ind = idxM(:) + (idxM(:).'-1)*nT
在重新起草之前，我的评论中确实有一些关于你有多喜欢它的评论！公平地说，它确实增加了关于该行的作用的明确性。
@Wolfie 不要忘记向后兼容 :)
@Wolfie 是的，这也是一个重要原因。事实上，当我在代码中使用隐式扩展时，我觉得有必要添加注释，因为很容易忽略正在发生的事情。所以bsxfun更清晰
@LuisMendo：这是我第一次看到“bsxfun 更清晰”——这是他们能想到的最晦涩、最迟钝的函数名称。从字面上看，我花了三年时间才记住这个名字。我很高兴我们不再需要它了。 :)

【解决方案2】：

使用cumsum 和diff 的解决方案。

[s,is] = sort(idxM);
sumM  = M(is,is);
idx = [diff(s)~=0 ,true];
CS = cumsum(sumM);
CS = cumsum(CS(idx,:),2);
n=sum(idx);
result = diff([zeros(n,1) diff([zeros(1,n); CS(:,idx)])],1,2);
sumM (:)=0;
sumM (s(idx),s(idx))=result;

【讨论】：

建议的解决方案非常聪明。然而，正如对其他提议的解决方案所观察到的那样，对于我的笔记本电脑上尺寸小于约 10,000x10,000 的矩阵 M，它仅与原始丑陋的 for 循环一样快，但使用了约 50% 的 RAM。对于较大的矩阵 M，我的笔记本电脑上的计算时间会激增，这可能是由于 RAM 的限制。

【解决方案3】：

我想向那些对此感兴趣的人指出另一个论坛上提供的答案

S=sparse(1:N,idxM,1); sumM=S.'*(M*S);

致谢（和有用的讨论）：

https://www.mathworks.com/matlabcentral/answers/407634-how-to-sum-parts-of-a-matrix-of-different-sizes-without-using-for-loops

【讨论】：

非常聪明..！