Как векторизовать цикл 'for' в Matlab

Question

Я пишу приложение Matlab, которое вычисляет какую-то матрицу. Я пытаюсь заменить циклы for в моей программе на основе векторных вычислений, однако, я застрял.

До сих пор я понял, что это просто образец часть кода:

kernel = zeros(5,5,5); 
offset = 3; 
sig=1; 

for x=-2:2 
    for y=-2:2 
     for z=-2:2 
      IN = x.^2/(sig^2) + y.^2/(sig^2) + z.^2/(sig^2); 
      kernel(offset+x, offset+y, offset+z) = exp(-IN/2); 
     end 
    end 
end 

можно заменить такой конструкции:

[x,y,z] = ndgrid(-2:2,-2:2,-2:2); 
IN = x.^2/(sig^2) + y.^2/(sig^2) + z.^2/(sig^2); 
kernel = exp(-IN/2); 

и дают одинаковые результаты. Но что, если мне нужно сделать некоторые небольшие изменения:

kernel = zeros(5,5,5); 
offset = 3; 
sig=1; 

%sample 3x3 matrix 
R=magic(3)/10; 

for x=-2:2 
    for y=-2:2 
     for z=-2:2 

      % calculation of new x, y, z 
      point = [x,y,z]*R; 

      IN = (point(1)^2)/(sig^2) + (point(2)^2)/(sig^2) + (point(3)^2)/(sig^2); 
      kernel(offset+x, offset+y, offset+z) = exp(-IN/2); 
     end 
    end 
end 

Как я могу ускорить эту конструкцию? Можно ли легко векторизовать? Я совершенно новичок в Matlab, поэтому я был бы признателен за любую помощь. Большое спасибо!

Answer 1

Один из вариантов - использовать arrayfun.

sig=1; 

%sample 3x3 matrix 
R=magic(3)/10; 

[x,y,z] = ndgrid(-2:2,-2:2,-2:2); 
kernel = arrayfun(@(x, y, z) exp(-(norm([x,y,z]*R/sig)^2)/2), x,y,z); 

Пояснение:

arrayfun принимает функцию, которая действует на скалярной вход (ы) для получения скалярного выхода (ов), а также массивы входов, чтобы перейти к функции. Затем он перебирает входные массивы, запускает вашу функцию на каждом из них и помещает вывод каждого в соответствующую запись в выходной матрице (-es). Итак, arrayfun в основном делает вложенный цикл, который замедляет все для вас.

В этом примере я также использовал anonymous function (также называемую лямбда-функцией) для выполнения работы во внутреннем цикле. Поскольку лямбда-функции должны быть одним выражением в Matlab, мне пришлось немного переписать внутренний цикл (используя простые алгебраические манипуляции). Если вам нужно использовать arrayfun с функцией, которая не легко выражается в виде лямбда, вы всегда можете записать эту функцию в отдельный файл .m и передать ее в arrayfun.

EDIT: обратите внимание, что вам больше не нужно назначать kernel, и вам также не нужно offset.