Функция векторизации MATLAB - PullRequest
Новая
       42

Функция векторизации MATLAB

2 голосов
/ 09 сентября 2011

У меня двойное суммирование по м = 1: М и n = 1: N для полярной точки с координатами rho , phi , z :

Double summ over m and n

Я написал ее векторизованную запись: 

N = 10;
M = 10;
n = 1:N;
m = 1:M;

rho = 1;
phi = 1;
z = 1;

summ =  cos (n*z)  * besselj(m'-1, n*rho) * cos(m*phi)';

Теперь мне нужно переписать эту функцию для приема векторов (столбцов) координат rho , phi , z .Я попробовал arrayfun, cellfun, simple for loop - они работают слишком медленно для меня.Я знаю о «советах и ​​приемах манипулирования массивами MATLAB», но, как новичок в MATLAB, я не могу понять repmat и другие функции.

Кто-нибудь может предложить векторизованное решение?

Ответы [ 2 ]

1 голос
/ 10 сентября 2011

Я думаю, ваш код уже хорошо векторизован (для n и m). Если вы хотите, чтобы эта функция также принимала массив значений rho / phi / z, я предлагаю вам просто обработать значения в цикле for, так как я сомневаюсь, что любая дальнейшая векторизация принесет значительные улучшения (плюс код будет сложнее читать).

Сказав, что в приведенном ниже коде я попытался векторизовать часть, в которой вы вычисляете различные умножаемые компоненты {row N} * { matrix N*M } * {col M} = {scalar}, путем одного вызова функций BESSELJ и COS (я помещаю каждую строку / матрицу / столбец в третьем измерении). Их умножение все еще выполняется в цикле (точнее, ARRAYFUN): 

%# parameters
N = 10; M = 10;
n = 1:N; m = 1:M;

num = 50;
rho = 1:num; phi = 1:num; z = 1:num;

%# straightforward FOR-loop
tic
result1 = zeros(1,num);
for i=1:num
    result1(i) = cos(n*z(i)) * besselj(m'-1, n*rho(i)) * cos(m*phi(i))';
end
toc

%# vectorized computation of the components
tic
a = cos( bsxfun(@times, n, permute(z(:),[3 2 1])) );
b = besselj(m'-1, reshape(bsxfun(@times,n,rho(:))',[],1)');             %'
b = permute(reshape(b',[length(m) length(n) length(rho)]), [2 1 3]);    %'
c = cos( bsxfun(@times, m, permute(phi(:),[3 2 1])) );
result2 = arrayfun(@(i) a(:,:,i)*b(:,:,i)*c(:,:,i)', 1:num);            %'
toc

%# make sure the two results are the same
assert( isequal(result1,result2) )

Я провел еще один тест с использованием функции TIMEIT (дает более точные значения времени). Результат согласуется с предыдущим: 

0.0062407    # elapsed time (seconds) for the my solution
0.015677     # elapsed time (seconds) for the FOR-loop solution

Обратите внимание, что при увеличении размера входных векторов оба метода будут иметь одинаковые временные характеристики (цикл FOR даже выигрывает в некоторых случаях)

0 голосов
/ 09 сентября 2011

Вам нужно создать две матрицы, скажем, m_ и n_, чтобы, выбрав элемент i,j каждой матрицы, вы получили желаемый индекс как для m, так и n.

Большинство функций MATLAB принимают матрицы и векторы и вычисляют их результаты поэлементно.Таким образом, чтобы получить двойную сумму, вы вычисляете все элементы суммы параллельно на f(m_, n_) и суммируете их.

В вашем случае (обратите внимание, что оператор .* выполняет поэлементное умножение матриц) 

N = 10;
M = 10;
n = 1:N;
m = 1:M;

rho = 1;
phi = 1;
z = 1;

% N rows x M columns for each matrix
% n_ - all columns are identical
% m_ - all rows are identical
n_ = repmat(n', 1,  M);
m_ = repmat(m , N,  1);

element_nm =  cos (n_*z) .* besselj(m_-1, n_*rho) .* cos(m_*phi);
sum_all = sum( element_nm(:) );
...