如何在不使用 for 循环的情况下对不同大小的矩阵的各个部分求和?
How to sum parts of a matrix of different sizes, without using for loops?
我有一个相对较大的矩阵 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
您可以按如下方式使用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
[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;
我想指出那些对另一个论坛上提供的答案感兴趣的人
S=sparse(1:N,idxM,1);
sumM=S.'*(M*S);
学分(和有用的讨论):
我有一个相对较大的矩阵 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
您可以按如下方式使用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
[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;
我想指出那些对另一个论坛上提供的答案感兴趣的人
S=sparse(1:N,idxM,1);
sumM=S.'*(M*S);
学分(和有用的讨论):