Reduce execution time in a code having 'while' loop

I have written a code for finding galaxy pairs. Array 'v1' contains 8 columns in which each column contains different properties of galaxies. 'v1' has 93000 rows and 8 columns. The code is attached below:-

c=3*10^5;
H0=67.4;
f=@(x) (0.315.*(1+x).^3+0.685).^(-1/2);
i=1;j=0;
while i<=size(v1,1);
    k=1;
    while k<=size(v1,1);
        if k~=i;
            
            dvel=c*abs(v1(i,2)/(1+v1(i,2))-v1(k,2)/(1+v1(k,2)))
            theta=acos(cosd(v1(i,4))*cosd(v1(k,4))*cosd(v1(i,3)-v1(k,3))+sind(v1(i,4))*sind(v1(k,4)));
            r=(c/H0)*integral(f,0,0.5*(v1(i,2)+v1(k,2)));
            rp=theta*r;
            if rp<=0.15;
                j=j+1;
                pair(j,:)=[v1(i,1) v1(k,1)];   
            end
            
        end
        k=k+1;
    end
    i=i+1
end

(1)c,H0 and f are some declared constant and functions respectively.

(2) size(v1,1)=93000

(3) For each value of 'i' here, 'k' loop is running (93000-1) times and therefore considering all values of 'i' , 'i' loop is running entirely for 93000 X (93000-1) times which actually takes a huge huge amount of time. Can anyone please suggest how to reduce the execution time of this loop, it would be of great help?

[Walter Roberson]

Question: what do you do with dvel ? You compute it, but I do not see any use of it.

Also, is there any column that the data is sorted on? With some algorithms, sorting on the second column might help speed things up, potentially allowing us to cut out some of the integrations. Not sure yet it is worth the overhead.

Also, what order of magnitude are the values in the second column ?

Answer 1

[Walter Roberson]

            if rp<=0.15;
                j=j+1;
                pair(j,:)=[v1(i,1) v1(k,1)];   
            end
            

Do not do that. You are growing the pair array in-place every time, and the total amount of time that takes increases with the sqquare of the number of entries in the array.

Instead, preallocate

 pairik = false(size(v1,1), size(v1,1));
 

and inside the loop,

            if rp<=0.15;
                pairok(i,k) = true;   
            end

After that, after the loops,

   [row, col] = find(pairok);
   pair = [row, col];
   

Also, I recommend converting the while to for loops.

I suspect you can do even better by using vectorization, but this change should make a really noticeable difference.

[Apashanka Das]

Sir thanks for this answer. Tried that but still it takes very very long time (e.g intitially for i=1 it takes nearly 1 min.). Sir can you please suggest any other method or algorithm for finding pairs as such execution time will be small.

[Walter Roberson]

The below is only partly complete ... I have added a clarifying question as a comment above

c = 3*10^5;
H0 = 67.4;
f=@(x) (0.315.*(1+x).^3+0.685).^(-1/2);

v12 = v1(:,2);
v13 = v1(:,3);
cv2o2i = c .* v12./(1+v12);
cv2o2k = cv2o2i.';
v13 = 
cosdv3i_j = cosd(v13 - v13.');
cosdv4i = cosd(v1(:,4));
cosdv4k = cosdv4i .';
sindv4i = sind(v1(:,4));
sindv4k = sindv4i .';
dvel = abs(cv2o2i - cv2o2k);  %needs 64 gigabytes
theta = acos( cosdv4i .* cosdv4k .* cpsdv3i_j + sindv4i .* sindv4k );

[Apashanka Das]

Sir thanks for this answer. Tried this but for this line 'cosdv3i_j = cosd(v13 - v13.')' error message is showing. Also for each element in 'v1' array (e.g ith term), the terms like theta,rp are to be checked with other elements except with the ith term itself.

[Walter Roberson]

please answer the questions in

https://www.mathworks.com/matlabcentral/answers/1583094-reduce-execution-time-in-a-code-having-while-loop#comment_1825964

I stopped trying to write the solution when I realized that I needed the answer to those questions.

Answer 2

[Apashanka Das]

Sir actually in the loop ' if rp<=0.15 ' should be replaced with ' if rp<=0.15 & dvel<=300 ' thats why dvel is calculated.

For every ith element dvel, theta, r, rp is calculated with the kth element (where k~=i), where i runs from 1 to 93000 and for each i, k also runs from 1 to 93000.

The second column values runs from 0 to 0.115.

[Apashanka Das]

Sir, it would be of great help if you suggest something in implementing this code in a much lesser time.

[Walter Roberson]

What accuracy levels do you need? Some of your constants only have single digit precision, but some of them have 3 digits precision, and your rp boundary only has two digits precision.

The reason I ask is that over that range of values in column 2, the integral is pretty much a straight line, and could be replaced by a very low order polynomial, if you are willing to accept that the occasional point might not be perfectly classified if it is right on the boundary. Remembering that your speed of light is not exactly accurate anyhow...

[Apashanka Das]

Sir for rp 4 digit precision will be fine.