is there any way to vectorize the code to speed up the calculation?

I have a function to calculation some lineshapes which was called thousands of times. it is pretty slow. I'm wondering how to make it run faster.maybe vectorizing the code to eliminate for loop? the relevant code is as following. thanks,

       % input variables, x and y,  with x =3001 points, fixed length 
x= -1.5:0.001:1.5; 
y = 1.18;
% function doing the calculation
N = 33; 
K = zeros(size(x));
a = zeros(1,N);
 summation = 0; 

    for n = 1:N
      a(n) = 2/9*exp(-((n-1)*pi/9)^2);
        first = (1i*(n-1)*pi*9+9^2*y)*(1-exp(-(1i*(n-1)*pi+9*y))*cos(9.*x)) + exp(-(1i*n*pi+12*y))*12^2.*x.*sin(9.*x);

        second = (1i*(n-1)*pi*9-9^2*y)*(1-exp(1i*(n-1)*pi-9*y).*cos(9.*x)) - exp(-(1i*n*pi+12*y))*12^2.*x.*sin(9.*x);

        summation = summation + a(n)*(first - second);
    end

    third = -(y-exp(-(9*y))*(y.*cos(9.*x)-x.*sin(9.*x)))./(x.^2+y^2);
    K = summation - a(1)*third;

[Marc]

What is tau?? I set it to 0.3 arbitrarily and got 0.01 seconds with the profiler (~0.008 to be exact) 2013b, MacOS 64bit. Seems fast to me.

[Chong Tao]

I changed my code. it is true, it dosen't take much time for one peak calculation. it is pretty slow, tens of minutes if it is called tens of thousands times.

Answer 1

[Walter Roberson]

You could use meshgrid() or ndgrid() to create meshes of n and x, and use those meshes in your code.

[xM, nM] = ndgrid(x, 1:N);

a = 2/9 * exp(-((nM-1)*pi/9).^2);
first = (1i*(nM-1)*pi*9+9^2*y).*(1-exp(-(1i*(nM-1)*pi+9*y)).*cos(9.*xM)) + exp(-(1i*nM*pi+12*y))*12^2.*xM.*sin(9.*xM);
second = (1i*(nM-1)*pi*9-9^2*y).*(1-exp(1i*(nM-1)*pi-9*y).*cos(9.*xM)) - exp(-(1i*nM*pi+12*y))*12^2.*xM.*sin(9.*xM);

Now (first - second) will be a matrix rather than a vector. Your code would then become something like

summation = sum( repmat(a, size(xM,1), 1) .* (first - second), 1);

but you might need to transpose "a" or repmat it along the second dimension instead of the first.

[Chong Tao]

Thank you Walter. sorry for the later reponse. I tried what you suggested. I used

        summation = sum( bsxfun(@times,a,bsxfun(@minus,first,second)), 2);

for the summation instead of repmat. I then tested the speed. It is about twice slower than the original for loop, which suprised me. Is that something you would expect?