在 MATLAB 中从非常大的数据集生成直方图的有效方法？

Question

我有两个最大为35,000*35,000 的二维数组：indices 和dotPs。由此，我想创建两个一维数组，使得pop 包含每个数字在indices 中出现的次数，nn 包含dotPs 中与这些数字相对应的元素的总和。我想出了以下（非常愚蠢）的方式：

dotPs = [81.4285    9.2648   46.3184    5.7974    4.5016    2.6779   16.0092   41.1426;
      9.2648   24.3525   11.4308   14.6598   17.9558   23.4246   19.4837   14.1173;
     46.3184   11.4308   92.9264    9.2036    2.9957    0.1164   26.5770   26.0243;
      5.7974   14.6598    9.2036   34.9984   16.2352   19.4568   31.8712    5.0732;
      4.5016   17.9558    2.9957   16.2352   19.6595   16.0678    3.5750   16.7702;
      2.6779   23.4246    0.1164   19.4568   16.0678   25.1084    6.6237   15.6188;
     16.0092   19.4837   26.5770   31.8712    3.5750    6.6237   61.6045   16.6102;
     41.1426   14.1173   26.0243    5.0732   16.7702   15.6188   16.6102   47.3289];

indices = [3     2     1     1     2     1     2     1;
           2     2     1     2     2     1     2     2;
           1     1     3     3     2     2     2     2;
           1     2     3     4     3     3     4     2;
           2     2     2     3     3     1     3     2;
           1     1     2     3     1     8     2     2;
           2     2     2     4     3     2     4     2;
           1     2     2     2     2     2     2     2];


nn = zeros(1,8);
pop = zeros(1,8);
uniqueInd = unique(indices);
for k=1:numel(uniqueInd)
    j = uniqueInd(k);
    [I,J]=find(indices==j);
    if j == 0 || numel(I) == 0
        continue
    end

    pop(j) = pop(j) + numel(I);
    nn(j) = nn(j) + sum(sum(dotPs(I,J)));
end

由于find 函数，这非常慢。我怎样才能更聪明地做到这一点，让它在几秒钟而不是几分钟内运行？

编辑：添加了用于测试代码的小型虚拟矩阵。

Answer 1

这两个任务都可以通过accumarray 函数完成：

pop = accumarray(indices(:), 1, [max(indices(:)) 1]).';
nn = accumarray(indices(:), dotPs(:), [max(indices(:)) 1]).';

这假定indices 只包含正整数。

---

编辑：

从 cmets，只应使用不带对角线的 indices 矩阵的下部，并且保证包含正整数。在这种情况下：

mask = tril(true(size(indices)), -1);
indices_masked = indices(mask);
dotPs_masked = dotPs(mask); 
pop = accumarray(indices_masked, 1, [max(indices_masked) 1]).';
nn = accumarray(indices_masked, dotPs_masked, [max(indices_masked) 1]).';

Answer 2

计算pop，可以使用hist，计算nn，我找不到聪明的解决方案（但我找到了不使用find的解决方案）：

pop = hist(indices(:), max(indices(:)));

nn = zeros(1,8);
uniqueInd = unique(indices);
for k=1:numel(uniqueInd)
    j = uniqueInd(k);
    nn(j) = sum(dotPs(indices == j));
end

必须有更好的解决方案来计算nn。

---

我发现了一个应用排序的更智能的解决方案。

我不确定它是否更快，因为对 35,000*35,000 个元素进行排序可能需要很长时间。

1. 排序indices 只是为了获取索引以将dotPs 按indices 排序。
2. 根据上一次排序返回的索引对dotPs进行排序。
3. cumsumPop = 计算pop 的累积和（indices 的直方图的累积和）。

4. cumsumPs = 计算已排序 dotPs 的累积总和。

5. 现在 cumsumPop 的值可以用作 cumsumPs 中的索引。
   因为 cumsumPs 是累积和，所以我们需要使用diff 以获得解决方案。

这里是“智能”解决方案：

pop = hist(indices(:), max(indices(:)));

[sortedIndices, I] = sort(indices(:));
sortedDotPs = dotPs(I);

cumsumPop = cumsum(pop);
cumsumPs = cumsum(sortedDotPs);

nn = diff([0; cumsumPs(cumsumPop)]);
nn = nn';

Answer 3

首先，请注意indices 的维度无关紧要（例如，如果indices 和dotPs 都是1D 数组或3D 数组，结果将是相同的）。

pop可以通过histcount函数计算出来，但是由于还需要计算dotPs数组对应元素的总和，问题就变得更难了。

这是一个可能的解决方案，带有for 循环。这个方案的好处是我没有在循环中调用find函数，所以应该更快：

%Example input
indices=randi(5,3,3);
dotPs=rand(3,3);

%Solution
[C,ia,ic]=unique(indices);
nn=zeros(size(C));
pop=zeros(size(C));
for i=1:numel(indices)
    nn(ic(i))=nn(ic(i))+1;
    pop(ic(i))=pop(ic(i))+dotPs(i);
end

此解决方案使用向量 ic 对每个输入值进行分类。之后，我遍历每个元素并更新nn(ic) 和pop(ic)。