Стабильный кластер в MATLAB

Question

Встроенная функция MATLAB accumarray принимает функцию fun в качестве четвертого аргумента.

A = accumarray(subs,val,sz,fun);

Это относится fun для каждого подмножества элементов в val которые имеют идентичные подписки в subs, Документация однако заявляет:

Если подписки в subs не отсортированы по своим линейным показателям, fun не должно зависеть от порядка значений во входных данных.

Как мы можем реализовать стабильную версию accumarray, который не имеет этого ограничения, но гарантирует, что подмножества принимают тот же порядок, что и val?

Пример:

subs = [1:10,1:10];
val = 1:20;
accumarray(subs(:), val(:), [], @(x)x(end)).'

Ожидаемый результат этого будет 11:20 если accumarray были стабильны. На самом деле вывод:

ans =
    11    12    13    14     5     6     7    18    19    20

Наша реализация должна дать:

accumarrayStable(subs(:), val(:), [], @(x)x(end)).'`
ans =
    11    12    13    14    15    16    17    18    19    20

Answer 1

Мы можем использовать sortrows как шаг предварительной обработки, чтобы сначала отсортировать индексы и соответствующие значения, как указано в документации:

SORTROWS использует стабильную версию быстрой сортировки.

Как подписчики в subs должны быть отсортированы по их линейным индексам, нам нужно отсортировать их в обратном лексикографическом порядке. Это может быть достигнуто путем переключения порядка столбцов до и после использования sortrows,

Это дает нам следующий код для стабильной версии accumarray:

function A = accumarrayStable(subs, val, varargin)
[subs(:,end:-1:1), I] = sortrows(subs(:,end:-1:1));
A = accumarray(subs, val(I), varargin{:});

Альтернатива:

Как полагает Луис Мендо, вместо sortrows можно также генерировать линейные индексы из индексов и использовать sort вместо.

function A = accumarrayStable(subs, val, varargin)
if numel(varargin)>0 && ~isempty(varargin{1})
    sz = varargin{1};
else
    sz = max(subs,[],1);
end
[~, I] = sort(subs*cumprod([1,sz(1:end-1)]).');
A = accumarray(subs(I,:), val(I), sz, varargin{:});

Обратите внимание, что мы должны использовать 1+(subs-1)*cumprod([1,sz(1:end-1)]).' для перевода в линейные индексы. Мы опускаем +1 а также -1 в результате sort останется прежним; что экономит нам несколько циклов.

Какое из приведенных выше решений быстрее, зависит от вашей машины и версии MATLAB. Вы можете проверить, например, через:

A = randi(10, 1e4, 5); 
timeit(@()accumarrayStable(A(:,1:end-1), A(:,end), [], @(x)x(1))