Как я могу ускорить этот вызов квантиля в Matlab? - PullRequest
Новая
       41

Как я могу ускорить этот вызов квантиля в Matlab?

8 голосов
/ 22 декабря 2011

У меня есть процедура MATLAB с одним довольно очевидным узким местом.Я профилировал функцию, в результате чего 2/3 вычислительного времени используется в функции levels:

enter image description here

Функция levels принимает матрицуof float и разбивает каждый столбец на nLevels сегментов, возвращая матрицу того же размера, что и входные данные, с каждой записью, заменяемой номером сегмента, в который она попадает.Функция 1011 * для получения пределов сегментов и цикл для назначения записей в сегменты.Вот моя реализация: 

function [Y q] = levels(X,nLevels)
% "Assign each of the elements of X to an integer-valued level"

p = linspace(0, 1.0, nLevels+1);

q = quantile(X,p);
if isvector(q)
    q=transpose(q);
end

Y = zeros(size(X));

for i = 1:nLevels
    % "The variables g and l indicate the entries that are respectively greater than
    % or less than the relevant bucket limits. The line Y(g & l) = i is assigning the
    % value i to any element that falls in this bucket."
    if i ~= nLevels % "The default; doesnt include upper bound"
        g = bsxfun(@ge,X,q(i,:));
        l = bsxfun(@lt,X,q(i+1,:));
    else            % "For the final level we include the upper bound"
        g = bsxfun(@ge,X,q(i,:));
        l = bsxfun(@le,X,q(i+1,:));
    end
    Y(g & l) = i;
end

Что я могу сделать, чтобы ускорить это?Можно ли векторизовать код?

Ответы [ 5 ]

4 голосов
/ 22 декабря 2011

Если я правильно понимаю, вы хотите знать, сколько предметов выпало в каждом ведре.Используйте:

n = исторических (Y, nbins)

Хотя я не уверен, что это поможет в ускорении.Это просто чище.

Редактировать: После комментария:

Вы можете использовать второй выходной параметр histc

[n, bin] = histc (...) также возвращает индексную матрицу bin.Если x - вектор, n (k) => sum (bin == k).bin - ноль для значений вне диапазона.Если x - это матрица M-by-N, то

2 голосов
/ 22 декабря 2011

Как насчет этого 

function [Y q] = levels(X,nLevels)

p = linspace(0, 1.0, nLevels+1);
q = quantile(X,p); 
Y = zeros(size(X));
for i = 1:numel(q)-1    
    Y = Y+ X>=q(i);
end

Это приводит к следующему: 

>>X = [3 1 4 6 7 2];
>>[Y, q] = levels(X,2)

Y =

     1  1  2  2  2  1

q =

     1  3.5  7

Вы также можете изменить логическую строку, чтобы значения были меньше, чем начало следующего бина. Однако я не думаю, что это необходимо.

2 голосов
/ 22 декабря 2011

Я думаю, вы должны использовать histc

[~,Y] = histc(X,q)

Как вы можете видеть в документе Matlab:

Описание

n = histc (x, ребра) подсчитывает количество значений в векторе x, попадающих между элементами в векторе ребер (которые должны содержать монотонно неубывающие значения).n - вектор длины (ребер), содержащий эти числа.Элементы x не могут быть сложными.

1 голос
/ 22 декабря 2011

Вы можете sort столбцы и разделить + округлить обратные индексы: 

function Y = levels(X,nLevels)
% "Assign each of the elements of X to an integer-valued level"
[S,IX]=sort(X);
[grid1,grid2]=ndgrid(1:size(IX,1),1:size(IX,2));
invIX=zeros(size(X));
invIX(sub2ind(size(X),IX(:),grid2(:)))=grid1;
Y=ceil(invIX/size(X,1)*nLevels);

Или вы можете использовать tiedrank: 

function Y = levels(X,nLevels)
% "Assign each of the elements of X to an integer-valued level"
R=tiedrank(X);
Y=ceil(R/size(X,1)*nLevels);

Удивительно, но оба эти решения немного медленнее, чем решение quantile + histc.

1 голос
/ 22 декабря 2011

Я сделал пару уточнений (в том числе вдохновленный Aero Engy в другом ответе), который привел к некоторым улучшениям.Чтобы проверить их, я создал случайную матрицу из миллиона строк и 100 столбцов для запуска улучшенных функций: 

>> x = randn(1000000,100);

Сначала я запустил мой неизмененный код со следующими результатами:

enter image description here

Обратите внимание, что из 40 секунд около 14 из них тратятся на вычисление квантилей - я не могу ожидать, чтобы улучшить эту часть процедуры (я предполагаю, что Mathworks уже оптимизировал ее,хотя я предполагаю, что для предположения получается ...)

Затем я изменил подпрограмму следующим образом, что должно быть быстрее и иметь преимущество в том, что меньше строк! 

function [Y q] = levels(X,nLevels)

p = linspace(0, 1.0, nLevels+1);
q = quantile(X,p);
if isvector(q), q = transpose(q); end

Y = ones(size(X));

for i = 2:nLevels
    Y = Y + bsxfun(@ge,X,q(i,:));
end

Результаты профилирования с этим кодом:

enter image description here

Так что это на 15 секунд быстрее, что представляет собой 150% -ное ускорение части моего кода, а не MathWorks.

Наконец, следуя совету Андрея (снова в другом ответе), я изменил код, чтобы использовать второй вывод функции histc, которая назначает записи для корзин.Он не обрабатывает столбцы независимо, поэтому мне пришлось циклически перебирать столбцы, но, похоже, он работает очень хорошо.Вот код: 

function [Y q] = levels(X,nLevels)

p = linspace(0,1,nLevels+1);

q = quantile(X,p);
if isvector(q), q = transpose(q); end
q(end,:) = 2 * q(end,:);

Y = zeros(size(X));

for k = 1:size(X,2)
    [junk Y(:,k)] = histc(X(:,k),q(:,k));
end

И результаты профилирования:

enter image description here

Теперь мы тратим всего 4,3 секунды в кодах вне функции quantile, котораяпримерно на 500% быстрее, чем я писал изначально.Я потратил немного времени на написание этого ответа, потому что я думаю, что он превратился в хороший пример того, как вы можете использовать профилировщик MATLAB и StackExchange в комбинации, чтобы получить гораздо лучшую производительность из вашего кода.

Ядоволен этим результатом, хотя, конечно, я продолжу рад слышать другие ответы.На этом этапе основное увеличение производительности произойдет за счет увеличения производительности части кода, которая в данный момент вызывает quantile.Я не вижу, как это сделать немедленно, но, может быть, кто-то еще здесь может.Еще раз спасибо!

...