Эффективный способ создания гистограммы из очень большого набора данных в 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 содержит только натуральные числа.

---

РЕДАКТИРОВАТЬ:

Из комментариев следует использовать только нижнюю часть матрицы 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 были одномерными или трехмерными массивами, результат будет одинаковым).

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).