Векторизовать/оптимизировать этот код в MATLAB?

Я создаю свою первую крупномасштабную программу MATLAB, и мне удалось написать оригинальный векторный код для всего, пока я не попытался создать изображение, представляющее плотность вектора в стереографической проекции. После пары неудачных попыток я пошел на сайт обмена файлами Mathworks и нашел программу с открытым исходным кодом, которая соответствует моим потребностям Malcolm Mclean. С тестовой матрицей его функция производит что-то вроде этого:

И хотя это почти то, что я хотел, его код основывается на тройном вложенном for-loop. На моей рабочей станции тестовая матрица данных размером 25000x2 заняла 65 секунд в этом разделе кода. Это неприемлемо, так как я буду масштабироваться до матриц данных размером 500000x2 в моем проекте.

До сих пор я мог векторизовать самый внутренний цикл (который был самым длинным/худшим), но я хотел бы продолжить и полностью избавиться от циклов, если это возможно. Вот оригинальный код Malcolm, который мне нужно оцифровать:

dmap = zeros(height, width); % height, width: scalar with default value = 32
for ii = 0: height - 1          % 32 iterations of this loop
    yi = limits(3) + ii * deltay + deltay/2; % limits(3) & deltay: scalars
    for jj = 0 : width - 1      % 32 iterations of this loop
        xi = limits(1) + jj * deltax + deltax/2; % limits(1) & deltax: scalars
        dd = 0;
        for kk = 1: length(x)   % up to 500,000 iterations in this loop
            dist2 = (x(kk) - xi)^2 + (y(kk) - yi)^2;
            dd = dd + 1 / ( dist2 + fudge); % fudge is a scalar
        end
        dmap(ii+1,jj+1) = dd;
    end
end

И вот это касается изменений, которые я уже внес в самый внутренний цикл (который был самым большим утечкой эффективности). Это сокращает время с 65 секунд до 12 секунд на моей машине для той же тестовой матрицы, которая лучше, но все же намного медленнее, чем хотелось бы.

     dmap = zeros(height, width);
    for ii = 0: height - 1
        yi = limits(3) + ii * deltay + deltay/2;
        for jj = 0 : width - 1
            xi = limits(1) + jj * deltax + deltax/2;
            dist2 = (x - xi) .^ 2 + (y - yi) .^ 2;
            dmap(ii + 1, jj + 1) = sum(1 ./ (dist2 + fudge));
        end
    end

So my main question, are there any further changes I can make to optimize this code? Or even an alternative method to approach the problem? I've considered using C++ or F# instead of MATLAB for this section of the program, and I may do so if I cannot get to a reasonable efficiency level with the MATLAB code.

Также обратите внимание, что на данный момент у меня нет ЛЮБЫХ дополнительных панелей инструментов, если бы я это сделал, тогда я знаю, что это было бы тривиально (например, используя hist3 из инструментария статистики).

Ответ 1

Запасное решение

yi = limits(3) + deltay * ( 1:height ) - .5 * deltay;
xi = limits(1) + deltax * ( 1:width  ) - .5 * deltax;
dx = bsxfun( @minus, x(:), xi ) .^ 2;
dy = bsxfun( @minus, y(:), yi ) .^ 2;
dist2 = bsxfun( @plus, permute( dy, [2 3 1] ), permute( dx, [3 2 1] ) );
dmap = sum( 1./(dist2 + fudge ) , 3 );

ИЗМЕНИТЬ

обработка чрезвычайно больших x и y путем разбиения операции на блоки:

blockSize = 50000; % process up to XX elements at once
dmap = 0;
yi = limits(3) + deltay * ( 1:height ) - .5 * deltay;
xi = limits(1) + deltax * ( 1:width  ) - .5 * deltax;
bi = 1;
while bi <= numel(x)
    % take a block of x and y
    bx = x( bi:min(end, bi + blockSize - 1) );
    by = y( bi:min(end, bi + blockSize - 1) );
    dx = bsxfun( @minus, bx(:), xi ) .^ 2;
    dy = bsxfun( @minus, by(:), yi ) .^ 2;
    dist2 = bsxfun( @plus, permute( dy, [2 3 1] ), permute( dx, [3 2 1] ) );
    dmap = dmap + sum( 1./(dist2 + fudge ) , 3 );
    bi = bi + blockSize;
end

Ответ 2

Это хороший пример того, почему запуск цикла из 1 вопроса. Единственная причина, по которой ii и jj инициируются в 0, - это убить термины ii * deltay и jj * deltax, которые, однако, вводят последовательность в индексировании dmap, предотвращают распараллеливание.. p >

Теперь, переписав петли, вы можете использовать parfor() после открытия matlabpool:

dmap = zeros(height, width);
yi   = limits(3) + deltay*(1:height) - .5*deltay;
matlabpool 8
parfor ii = 1: height
    for jj = 1: width
        xi    = limits(1) + (jj-1) * deltax + deltax/2;
        dist2 = (x - xi) .^ 2 + (y - yi(ii)) .^ 2;
        dmap(ii, jj) = sum(1 ./ (dist2 + fudge));
    end
end
matlabpool close

Имейте в виду, что открытие и закрытие пула имеет значительные накладные расходы (10 секунд на моем Intel Core Duo T9300, Vista 32 Matlab 2013a).

PS. Я не уверен, что внутренний цикл вместо внешнего может быть значимо распараллелен. Вы можете попытаться переключить parfor на внутренний и сравнить скорости (я бы рекомендовал перейти к большой матрице сразу же, поскольку вы уже работаете за 12 секунд, а накладные расходы почти такие же).

Ответ 3

В качестве альтернативы эту проблему можно решить с помощью методов оценки плотности ядра. Это часть инструментария Statistics Toolbox, или эта реализация KDE Здравко Ботева (никаких инструментов не требуется).

В приведенном ниже примере кода я получаю 0,3 секунды для N = 500000 или 0,7 секунды для N = 1000000.

N = 500000;
data = [randn(N,2); rand(N,1)+3.5, randn(N,1);];  % 2 overlaid distrib
tic; [bandwidth,density,X,Y] = kde2d(data); toc;
imagesc(density);