加速迭代 - MATLAB

Question

考虑两个向量 A = [20000000 x 1] 和 B = [20000000 x 1 ]。

我需要找到所有与B中每个唯一元素对应的A的总和。

尽管看起来非常简单，但在MATLAB中却需要很长时间。

目前，我正在使用：

u = unique(B);
length_u = length(u);
C = zeros(length_u,1);

for i = 1:length_u
   C(i,1) = sum(A(B==u(i)));
end

有没有任何方法可以使其运行更快？我尝试了将循环拆分并使用并行计算工具箱中的2个parfor循环来运行（因为我只有2个内核）。仍需要几个小时。

注：是的，我应该买一台更好的电脑。

- enigmae

---

瓶颈在哪里？是“unique”还是循环？ - Shai

你能发一下你的A和B矩阵吗？例如，A=randi(100,[20000000,1]);等。 - Autonomous

@Shai 这是一个循环。运行了如此多的迭代，Matlab可能会感到困惑。 - enigmae

如果我对你的代码进行profile，它对我来说是独一无二的。 - bdecaf

@bdecaf，你用什么来表示A和B？如果A的唯一元素数量相对较小，那么循环就不再是瓶颈了，但如果A很大且有许多唯一元素，则unique变得可以忽略。 - Shai

Answer 1

您首先必须查看这个答案。
如果必需，您可以使用histc和accumarray的组合。

A = randi( 500, 1, 100000 );
B = randi( 500, 1, 100000 );

ub = unique( B );

[ignore idx] = histc( B, [ub-.5 ub(end)+.5] );
C = accumarray( idx', A' )';

可以在ideone中查看玩具比较和naive for-loop实现。

它是如何工作的？

我们使用histc的第二个输出来将B（以及后来的A）的元素映射到由ub（B的唯一元素）定义的箱中。
然后，使用accumarray根据idx定义的映射对A的所有条目进行求和。
注意：假设B的唯一元素至少相差0.5。

- Shai

---

@Shai，你的解决方案中唯一不太好的地方就是你在“注释”中写的事实。但真的做得很好。我没有想到histc。+1 - The Minion

@TheMinion 可以通过查找“ub”中最小的绝对差来避免 +0.5，但我想保持代码简单。 - Shai

您可以（并且应该）使用 ~ 而不是忽略。 - Ander Biguri

@AnderBiguri 你说得对。然而，我习惯于旧版本的Matlab，它不支持~ :( ... 老习惯难改。 - Shai

Answer 2

根据Shai建议，进一步简化代码：

A = randi( 500, 1, 100000 );
B = randi( 500, 1, 100000 );

[~,~,idb] = unique( B );

C = accumarray( idb', A' )';

这里的 "idb" 给出了一个与 "idx" 向量相同的向量，这与 Shai 提供的代码相同。

- Naveen

Answer 3

如果B只包含整数，你可以轻松地在一行代码中完成它，使用sparse添加相同索引的元素的事实:

C = nonzeros(sparse(B,1,A));

- Luis Mendo

---

整数的好解决方案。我不知道sparse可以这样做。 - The Minion

@TheMinion 谢谢！是的，这不是 Matlab 的一个很出名的特性。它更多地将 sparse 转换为带有 issparse 选项的 accumarray。 - Luis Mendo

在这方面，accumarray 是 sparse 的扩展。使用 accumarray 与 histc 一起可以推广到非整数。 - Shai

Answer 4

我修改了求和的方式。不再需要检查每个元素是否符合条件 (B==u(i))，而是对数组排序并在元素改变时停止遍历。从那个元素开始下一个求和过程。这样我只需要循环一次A中的每个元素，而不是length_u次。以下是我使用的代码:

A= rand(100000,1);
B= round(rand(100000,1)*25000);
u = unique(B);
length_u = length(u);
C = zeros(length_u,1);
E = zeros(length_u,1);
tic;
for k = 1:length_u
   C(k,1) = sum(A(B==u(k)));
end
t_OP=toc;

tic
D= sortrows([A,B],2);
n=1;
for l=1:numel(u)
    m=n;
    while m<numel(B) && D(m+1,2)==u(l) 
        m=m+1;
    end
    E(l,1) = sum(D(n:m,1));
    n=m+1;
end
t_trial=toc;
display(t_OP)
display(t_trial)

我也使用了你的代码。你的代码运行时间为：t_OP=10.9398，我的修改后运行时间为：t_trial=0.0962。希望这有所帮助。我通过构建sum(E-C)来确保代码能够正常工作，结果为0。
编辑：速度测试
我还将其与@Shai的解决方案进行了比较。结果是：

t_OP =

   10.8147
t_trial =

    0.0984
t_Shai =

    0.0154

编辑：@moarningsun的评论
如果您在构建总和之前对数组进行排序，则可以使用unique的第二个输出，而无需使用while循环。

tic
A = randi( 25000, 1, 100000 );
B = randi( 25000, 1, 100000 );
D= sortrows([A',B'],2);
[u, idx] = unique(D(:,2));
idx = [idx; numel(D(:,2))+1];
for l=1:numel(u)
    E(l,1) = sum(D(idx(l):idx(l+1)-1,1));
end
t_trial=toc;

- The Minion

---

你的方法与 histc 和 accumarray 相比如何？ - Shai

@Shai 我编辑了我的帖子。我的修改导致了一个100倍的因素，而你的修改导致了额外的10倍因素。所以比OP快了大约1000倍。 - The Minion

你也可以先将 B 排序，然后从 [i, idx] = unique(B_sorted) 获取“边界”索引。这样可以避免使用 while 循环。 - user2379410

@moarningsun 你说得对。我尝试了，但实际上并没有改变运行时间。我会将其作为第二次编辑发布以供日后参考。 - The Minion