Vectorizing MATLAB function

2
votes

I have double summation over m = 1:M and n = 1:N for polar point with coordinates rho, phi, z:

Double summ over m and n

I have written vectorized notation of it:

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)';

Now I need to rewrite this function for accepting vectors (columns) of coordinates rho, phi, z. I tried arrayfun, cellfun, simple for loop - they work too slow for me. I know about "MATLAB array manipulation tips and tricks", but as MATLAB beginner I can't understand repmat and other functions.

Can anybody suggest vectorized solution?

2
matlabsumvectorizationbessel-functions
Can you elaborate on the dimensions of each variable. I.e. rho is 1xA, etc. and what dimensions you expect for the output. First of all, this helps us to help you and secondly, this helps you to help yourself as proper dimensions are the first thing to consider when using MATLAB. – Egon

2 Answers

1
votes

I think your code is already well vectorized (for n and m). If you want this function to also accept an array of rho/phi/z values, I suggest you simply process the values in a for-loop, as I doubt any further vectorization will bring significant improvements (plus the code will be harder to read).

Having said that, in the code below, I tried to vectorize the part where you compute the various components being multiplied {row N} * { matrix N*M } * {col M} = {scalar}, by making a single call to the BESSELJ and COS functions (I place each of the row/matrix/column in the third dimension). Their multiplication is still done in a loop (ARRAYFUN to be exact):

%# 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) )

I did another benchmark test using the TIMEIT function (gives more fair timings). The result agrees with the previous:

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

Note that as you increase the size of the input vectors, the two methods will start to have similar timings (the FOR-loop even wins on some occasions)

0
votes

You need to create two matrices, say m_ and n_ so that by selecting element i,j of each matrix you get the desired index for both m and n.

Most MATLAB functions accept matrices and vectors and compute their results element by element. So to produce a double sum, you compute all elements of the sum in parallel by f(m_, n_) and sum them.

In your case (note that the .* operator performs element-wise multiplication of matrices)

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(:) );