How to code a faster loop

I have a program that running loop, i=1:n, where n in the magnitude of 10'000. I found out that this part of the loop could be the main cause of slowing down the program significantly. Before adding this part overall I only need on average less than 3min. Now, I need 3.5hrs.

Part of the formulation basically requires us to deduct current value (i), with the value before (i-1). While running, I see repeatiton of e.g s(i) = columns 5 through 10 all the way to 1000 (n value). Same goes to ds(i). Not sure if this is avoidable and seems to be the root of the long computation time.

clc, clear all
ntrial=1000;
for ntl = 1:ntrial
    int=0.8
    ds1=0
    ds2=0.2
    n=1000. 
    x(1) = 0;
    y(1) = 0;
    
    xt(1)=2;
    yt(1)=3;
    
    un=0;
    en=0;
    ms=0;
    
    % s(1)=0
    SDEG=90;
    for i=1:n;
        dist(i) = sqrt(((x(i)-xt(i))^2)+((y(i)-yt(i))^2));
        % dist(i) = 5+x(i)+y(i);
        %******* get angle between 2 points *****
        angleDEG=atand((yt(i)-y(i))/(xt(i)-x(i)))
        angleRAD=angleDEG/180*pi
        
        %******* get s value *************************************************
        if angleDEG<=SDEG
            s(i)=(900/dist(i))*exp(-2*(angleRAD)^2)+int
            un=un+1
        else
            s(i)=int
            en=en+1
        end
        if (i==1) % When i=1, i-1 = 0 dont exist.we use this to defined the initial condition
            % s(i-1)=s(i)
            ds(1)=0
            s1=ds1+s(i)
            s2=ds2+s(i)
        else
            ds(i)= s(i)-s(i-1)  % Now i>1, we can procedd with the usual formulation
            s1=ds1+s(i-1)       % As ds= s(now)-s(previous)
            s2=ds2+s(i-1)
        end
        x(i+1)=x(i);
        y(i+1)=y(i);
        
        xt(i+1)=xt(i);
        yt(i+1)=yt(i);
        if ds(i)>=0.2;
            ms=1;
        end
    end
    data(ntl,:)=[ms un en];
end
toc

Thank you for the help.

Answer 1

[dpb]

The biggest slowdown in your code above is you haven't preallocated any of the arrays x, y, xt, yt, s, data so they're all being dynamically allocated and reallocated every time the next loop index value is created.

Also, it doesn't look to me as though the inner loop is any different the 2:1000 times it runs than it is the first...you reset x,y and xt, yt back to the same starting values each pass through the outer loop. What's the point in that?

>> data

data =

     0    10     0
     0    10     0
     0    10     0
     0    10     0
     0    10     0
     0    10     0
     0    10     0
     0    10     0
     0    10     0
     0    10     0

>> 

I cut the two loops down to 10 for the above to illustrate -- the output array data is constant. That will be the same for your above code for all 1000 trials.

[Julian Tan]

Many thanks for the reply.

Sorry for the confusion. It is part of the much larger program with many condition set in for x, y, xt, yt, s, all of which depending on a certain outcome and random generator during each of the inner loop 'i'.

The outer loop ntrial is basically running the same simulation for 1000 times to collect the 'probablitic' outcome.

The codes runs well before introduction of s(i) and ds(i). Just manage to cut down the time by half by defining 'ds' as a constant number during that given 'i' instead of ds(i) which it will keep as an array.

Now left the s(i) which I cant define it as a constant number during that particular 'i' since I will need the s(i-1) during (i) to define 'ds'.

I am stil new to matlab and learning the code, thinking of coding something like this to avoid the need of s(i), e.g

During i=1

kk= s %value of s when i=1

Then when i=2,

ds= s -kk

%where s now is value of s when i=2 and kk was the value of s when i=1 ,hence achieve my needs for

ds= s(i)-s(i-1)

Still trying out the code/logic.

Cheers

[dpb]

OK, I figured probably something was generating a PRNV somewhere, but wasn't shown in posted code.

"The codes runs well before introduction of s(i) and ds(i)."

Again, the biggest problem still in run time is that you haven't preallocated all the arrays you're creating dynamically before the loop.

As you note, anything you can compute dynamically as a single value that isn't mandatory to be retained for the entire loop is a plus; altho you could see the difference wouldn't be nearly so much if that array was preallocated (as long as you have enough memory to not have to page memory, anyway).

Insert the lines

    % s(1)=0
    ...
    % preallocate -- and any others you must have full vectors for
    s=zeros(n,1); ds=s; x=s; xt=s; y=s; yt=s;
    for i=1:n;
        d=sqrt(((x(i)-xt(i))^2)+((y(i)-yt(i))^2));
        % dist(i) = 5+x(i)+y(i);
        angleRAD=deg2rad(atand((yt(i)-y(i))/(xt(i)-x(i))));
        
        %******* get s value *************************************************
        if angleDEG<=SDEG
            s(i)=(900/d)*exp(-2*(angleRAD)^2)+int;
            un=un+1;
        else
            s(i)=int;
            en=en+1;
        end
        ...

Didn't look like any need to save the dist array; I converted it to d above...again, look at the final algorithm and eliminate anything that doesn't have to have the full array as the result from subscripting -- you can avoid for the difference between two terms by just keeping the previous and current and reassigning them each pass.

W/O the full specification of what the simulation is doing it's not simple to see if the whole loop could be vectorized if you just first generated the vector or random values to use. For example, terms like (yt(i)-y(i) can be written as yt-y for the whole vector if one can build y, yt in place. But, the present code just duplicates those values from preceding so that it's not possible to know what must really happen--it's just a repmat(xt(1),n,1) now to duplicate the sample code. That, more than likely, isn't the real answer, though... :)

Here's example of the difference --

This is a very basic timings loop that looked something like

while N>100
  ds=zeros(N,1);
  tic;ds(1)=0;
  for i=2:N,ds(i)=s(i)-s(i-1);end;dt=toc;
  fprintf('%10d %e\n',N,dt);
  N=N/10;
end

for the preallocated case where as the other had

  clear ds;

in lieu of the allocation step. For the simple variable, it's approximately linear for each but the timing ratio shows a bigger hit with size, varying from 20X to 40X as fast as N goes up. The variation can become much more extreme with more complex data types as well.

  >> timings(:,2)./timings(:,4)

ans =

   18.4897
   19.9982
   18.2296
   22.8117
   35.4246
   39.7234

>> 

Answer 2

[Julian Tan]

WOW, many thanks for the effort and graph presentation! Will definitely learn a great deal from your sharing. Will go through it and learn all the tricks! Learn a great deal on the important of 'preallocation'.

Yes you are right, I tried as much as possible to eliminate those don't need to be kept. And yes the dist(i) too. Manage to further reduced the time, at first to around 3hours then 2hrs and now only need 15-18mins! Amazing!

This is the 'rought' test version of the code. I removed the outer loop. Previously take around 10s to compute the same code. now only at around 1.7s.

clc, clear all
tic
    
int=0.8
ds1=0
ds2=0.2
n=1000. 

x(1) = 0;
y(1) = 0;

xt(1)=2;
yt(1) =3;

SDEG=90;

                    %*** Initial condition for Xsi (to be value of previous 's')
                    dist = sqrt(((x(1)-xt(1))^2)+((y(1)-yt(1))^2));
                    angleDEG=atand((yt(1)-y(1))/(xt(1)-x(1)))
                    angleRAD=angleDEG/180*pi

                    if angleDEG<=SDEG
                        s(1)=(900/dist)*exp(-2*(angleRAD)^2)+int
                    else
                        s(i)=int
                    end
                    
                    Xsi=s
                    %*** Initial condition for Xsi ****** end
          
        for i=1:n
 
                    dist = sqrt(((x(i)-xt(i))^2)+((y(i)-yt(i))^2));
                    %******* get angle between UAV and target *****
                    angleDEG=atand((yt(i)-y(i))/(xt(i)-x(i)))
                    angleRAD=angleDEG/180*pi
                    
                    %******* get s value *************************************************
                    if angleDEG<=SDEG
                        s=(900/dist)*exp(-2*(angleRAD)^2)+int
                    else
                        s=int
                    end
                  
                    if (i==1)    % When i=1, i-1 = 0 dont exist. so we need to defined the initial condition
                        ds=0
                        s1=ds1
                        s2=ds2
                    else
                        ds= s-Xsi
%                         ds(i)= s(i)-s(i-1)  % Now i>1, we can procedd with the usual formulation
                        s1=ds1+Xsi       % As ds= s(now)-s(previous)
                        s2=ds2+Xsi
                    end
                    
                   AA_DATAcheck(i,:)=[ds s Xsi]; %DATA checker
                    
                    Xsi=s % Renew Xsi for next loop, anything above this is previous value of Xsi
                    
                    

        x(i+1)=x(i);
        y(i+1)=y(i);
        
        xt(i+1)=xt(i);
        yt(i+1)=yt(i);
        
        
        
        end
        toc