Can I vectorize this any further?

Hello,

I have such a code:

maxN = 100;
xElements = 200;
XY = ones(xElements, xElements); % for simplicity, if fact there are more complicated doubles
%
for i = 1 : xElements
    for j = 1 : xElements
        nVec = 0:maxN;
        besselForCurrentN = besselj(nVec+1, 2*pi*XY(i, j));
%         (...)
    end
end

How to vectorize it also for i and j?

Best regards, Alex

Answer 1

[Roger Stafford]

My understanding is that - I admit I haven’t tried this - if you use the ‘bsxfun’ with its ‘fun’ defined as ‘besselj’ and the two arguments “reshape(nVec+1,1,1,[])” and “2*pi*XY”, you will get a 3D array for 'besselForCurrentN' as the result in accordance with nVec+1 varying along the third dimension. An alternative would be to use ‘repmat’ and ‘reshape’ on both nVec+1 and 2*pi*XY in such a way that both become 3D arrays of the same size and then call on ‘besselj’ with these as inputs - I believe ‘besselj’ will accept two 3D arrays if they are the same size.

Whether you can use such a 3D result in a vectorized manner in your further computations is something only you can determine.

[Alex Kurek]

Thanks. The first (bsxfun) approach is not much faster, than loops:

Elapsed time is 1.659047 seconds.

Elapsed time is 1.560014 seconds.

I have problems veryfying if I did the second correctly.

[Roger Stafford]

What I had in mind is:

   [m,n] = size(XY);
   B = besselj(repmat(reshape(nVec+1,1,1,[]),m,n,1),repmat(2*pi*XY,1,1,length(nVec)));

However, I doubt that this will be as fast as using ‘bsxfun’.

[Alex Kurek]

I tried, the speed is identical.