How to vectorize a for loop?

How can i make that code run faster?

 data_file=importdata('Hz.txt');
xA=data_file(:,1); 
yA=data_file(:,2);
zA=data_file(:,3); 
[ndat,yyy]=size(xA);  #ndat=540

aA=1./(xA+1); %Number of randoms Ran=10000;

 X1=randi([-100,150],[1,Ran])/100;
s1=randi([-100,150],[1,Ran])/100;
s2=randi([-100,150],[1,Ran])/100;
sw=randi([50,130],[1,Ran]);
T=zeros([ndat,Ran]);
H=zeros([ndat,Ran]);
XsA=zeros([1,ndat]);
Chit=zeros([1,Ran]);
t=0:0.0000000001:0.04;

 for j=1:Ran

    x1=X1(j);
    y1=s1(j);
    y2=s2(j);
    w=(sw(j));
    x2=-x1-1/4*y1^2-3/2*y1*y2-1/4*y2^2;

    As=x1*y1*y2*exp(-w*t)+x2;

    for i=1:ndat
         dA=abs(As-aA(i));

        [v T(i,j)]=min(dA);
        I=T(i,j);

        H(i,j)=x1*y2+exp(-w*t(I))+exp(x2); 

        XsA(i)=(H(i,j)-yA(i)).^2 ./(zA(i).^2);
        Chit(j)=sum(XsA);
    end

end

end

[Jan]

An optimization needs a look in the profiler (although this disables the important JIT!) and the possibility to run the code. Even advanced programmers cannot predict, how the JIT will handle modifications of the code. Therefore it would be helpful, if you provide test data. E.g. RAND with a fixed seed is a good and cheap method in the forum - if such data are valid for the tests.

Answer 1

[Jan]

At first follow the standard rule to avoid all repeated calculations by moving them befor the loop:

 c1 = exp(-w * t);
 c3 = 1 ./ (zA .^ 2);
 for j=1:Ran
    x1 = X1(j);
    y1 = s1(j);
    y2 = s2(j);
    w  = sw(j);
    x2 = -x1 - 0.25 * y1^2 - 1.5 * y1 * y2 - 0.25 * y2^2;
    As = x1 * y1 * y2 * c1 + x2;
    c4 = x1 * y2 + exp(x2);
    for i = 1:ndat
      dA     = abs(As - aA(i));
      [v, I] = min(dA);
      T(i,j) = I;
      H(i,j) = c4 + c1(I); 
      XsA(i) = (H(i,j) - yA(i)) .^ 2 .* c3(i);
    end
    Chit(j) = sum(XsA);  % I guess this should be here
  end

From this point a vectorizeation would mean to replace "i" inside the loop by "1:ndat".

[Tony]

the Chit is ok. and also c4 and c3. But w,t are vectors so c1 is not working. Also i have the same code for another problem where inside the 1:ndat loop there is a quad command. i.e. H=quad(As,0,t(I))*aA(i).

Answer 2

[Doug Hull]

400 Million long vector t seems excessive to me.

What if you make it significantly shorter? I would lower your number of iterations on everything. See if that gives good results. Then start cranking up the number of iterations, and resolution. If it is not fine enough, then run it through the profiler to see where the bottlenecks are. My guess is that doing the above will get the results you need without code changes.

[Tony]

When i define a smaller vector t , i don't have accuracy in the dA. But if that cost, i think that will be ok if i can set 4M long t.

Answer 3

[Robert Cumming]

you over write:

   Chit(j)=sum(XsA);

in every i loop - are you sure you want to do that?

Profile the code and you will see which lines are the most expensive.

[Tony]

I have correct that.

[Jan]

Where have you corrected this? Is sum(XsA) calculated inside or outside the for i loop?