В поисках эффективного способа вычисления - Matlab

Question

У меня есть скалярная функция f([x,y],[i,j])= exp(-norm([x,y]-[i,j])^2/sigma^2), которая получает два 2-мерных векторы в качестве входных данных (здесь норма реализует евклидову норму). Значения x, i находятся в диапазоне 1:w, а значения y,j - в диапазоне 1:h. Я хочу создать массив ячеек X так, чтобы X{x,y} содержал матрицу awxh, такую ​​что X{x,y}(i,j) = f([x,y],[i,j]). Очевидно, что это можно сделать, используя 4 вложенных цикла:

for x=1:w;
for y=1:h;
    X{x,y}=zeros(w,h);
    for i=1:w
        for j=1:h
            X{x,y}(i,j)=f([x,y],[i,j])
        end 
    end
end
end

Это, однако, крайне неэффективно. Я был бы очень признателен за эффективный способ создания X.

Answer 1

Один из способов сделать это - удалить 2 самые внутренние петли и заменить их векторной версией. Судя по вашей функции f это не должно быть слишком плохо

Сначала нам нужно построить две матрицы, содержащие от 1 до w в каждой строке и от 1 до h в каждом столбце, например:

wMat=repmat(1:w,h,1);
hMat=repmat(1:h,w,1)';

Это будет представлять две внутренние петли, и транспонирование позволит нам получить все комбинации. Теперь мы можем векторизовать вычисление (f([x,y],[i,j])= exp(-norm([x,y]-[i,j])^2/sigma^2)):

for x=1:w;
    for y=1:h;
        temp1=sqrt((x-wMat).^2+(y-hMat).^2);
        X{x,y}=exp(temp1/(sigma^2));
    end
end

Где мы вычислили евклидову норму для всех пар узлов во внутренних циклах одновременно.

Answer 2

Некоторое обсуждение и код

Хитрость здесь в том, чтобы выполнить norm-calculations с числовыми массивами и сохранить результаты в версии массива ячеек как можно позже. Для выполнения norm-calculations Вы можете воспользоваться помощью ndgrid, bsxfun и немного permute + reshape чтобы придать ему "форму", необходимую для окончательной версии массива ячеек. Итак, вот векторизованный подход для выполнения этих задач -

%// Create x-y/i-j values to be used for calculation of function values
[xi,yi] = ndgrid(1:w,1:h);

%// Get the norm values
normvals = sqrt(bsxfun(@minus,xi(:),xi(:).').^2 + ...
                                bsxfun(@minus,yi(:),yi(:).').^2);
%// Get the actual function values
vals = exp(-normvals.^2/sigma^2); 

%// Get the values into blocks of a 4D array and then re-arrange to match
%// with the shape of numeric array version of X
blks = reshape(permute(reshape(vals, w*h, h, []), [2 1 3]), h, w, h, w);
arranged_blks = reshape(permute(blks,[2 3 1 4]),w,h,w,h);

%// Finally get the cell array version
X = squeeze(mat2cell(arranged_blks,w,h,ones(1,w),ones(1,h)));

Бенчмаркинг и время выполнения

После улучшения исходного зацикленного кода с предварительным выделением для X и функциональный f, runtime-benchmarks были выполнены с ним против предложенного векторизованного подхода с размерами данных как w, h = 60 и результаты выполнения, полученные таким образом, были:

----------- With Improved loopy code
Elapsed time is 41.227797 seconds.
----------- With Vectorized code
Elapsed time is 2.116782 seconds.

Это наводит на мысль о том, что 20x Ускорение с предложенным решением!

---

Для чрезвычайно больших размеров данных

Если вы имеете дело с огромными размерами данных, по сути, вы не даете достаточно памяти для bsxfun работать и bsxfun Известно, что он использует много памяти для создания эффективного векторизованного решения. Таким образом, для таких случаев с огромным размером данных вы можете использовать следующий циклический подход для замены normvals расчеты, которые были перечислены в предыдущем bsxfun решение на основе -

%// Get the norm values
nx = numel(xi);
normvals = zeros(nx,nx);
for ii = 1:nx
    normvals(:,ii) = sqrt( (xi(:) - xi(ii)).^2 + (yi(:) - yi(ii)).^2 );
end

Answer 3

Мне кажется, что когда вы пробегаете цикл для x=w, y=h, вы рассчитываете все значения, которые вам нужны одновременно. Поэтому вам не нужно пересчитывать их. Как только у вас есть это:

for i=1:w
    for j=1:h
        temp(i,j)=f([x,y],[i,j])
    end 
end

Тогда, например X{1,1} просто temp(1,1), X{2,2} просто temp(1:2,1:2), и так далее. Если вы можете векторизовать расчет f (norm вот только евклидова норма этого вектора?) тогда станет еще проще.