Optimising Nearest Neighbor Program

Hi Guys,

I am trying to optmimise this code so that it runs in under 10 seconds for N=20k.

It currently takes around 40 seconds to run.

I think that I need to vecotirse some or all of the loops so that the calculations are done at the same time, but I cannot figure out how to do it.

Any help would be much appreciated as we are learning from home with next to no support from the University.

Here is the code.

s = zeros(N,N);
for n = 1:3
    for a = 1:N
        count = 0;
        for b=1:N
            off=0;
            if (pos(n,a)<=0.25 && pos(n,b) >= 0.75)
                off=1;
            elseif pos(n,a)>=0.75 && pos(n,b)<=0.25
                off=-1;
            end
            s(a,b)=s(a,b)+(pos(n,a)-pos(n,b)+off)^2.0;
            if n == 3
                s(a,b)=sqrt(s(a,b));
            end
        end
    end
end
match = zeros(1,N);
for a=1:N
    mindist=1e10;
    for c=1:N
        if  (a~=c)
            mindist=min(s(a,c),mindist);
            if (mindist==s(a,c))
                match(a)=c;
            end
        end
    end
end

[KSSV]

Read about knnsearch.

[Samuel Leeney]

Thanks, I know there are functions that specifically do the job out there, but this is a learning excersise so if possible I would like to practice optimising my own code.

[KSSV]

Refer this, it might help you: https://in.mathworks.com/matlabcentral/fileexchange/44334-nearest-neighboring-particle-search-using-all-particles-search

[Samuel Leeney]

Thanks a lot I have had a look through the link and it does help me understand how the nearest neighbor works - the purpose of my question though is to learn about optimisation so I am really just trying to find a way to improve the efficiency of the code above (or any future code I write) rather than write a nearest neighbor program from scratch.

Answer 1

[Bruno Luong]

I put my comments as answer here, so you can accept if it helps

N=20000;
nd=3;
pos=rand(nd,N);

pos_r = reshape(pos,[nd 1 N]);
s = zeros(N);
for b=1:N
    posb = pos_r(:,:,b);
    off = (pos<=0.25 & posb>=0.75) - ...
          (pos>=0.75 & posb<=0.25); 
    sb = sum((pos-posb+off).^2,1);
    s(:,b) = sb(:);
end
s = sqrt(s);

s(1:N+1:end) = Inf;
[~,match] = min(s,[],2);
match = match.'; % row vector

Still if one dosn't have to deal with the odd OFFSET, the delaunay approach is much faster.

Answer 2

[KSSV]

N = 1000;
seed = 1234;
%seed=input('Random number seed = ');
rng(seed)
pos=rand(3,N);

a = 1:N ;
b = 1:N ;

[a,b] = meshgrid(a,b) ; 

S = zeros(N,N) ; 
for n = 1:3
  off = zeros(N,N) ;
  P =  pos(n,:) ; 
  off(P(a)<=0.25 && P(b) >= 0.75) = +1 ; 
  off(P(a)>=0.75 && P(b)<=0.25) = -1 ; 
  S = S+(P(a)-P(b)+off).^2;
end
S = sqrt(S); 

[Samuel Leeney]

Thanks, but it doesn't run at all like this:

"Operands to the || and && operators must be convertible to logical scalar values.

Error in Untitled (line 25)

off(pos(n,a)<=0.25 && pos(n,b) >= 0.75) = +1 ; "

[KSSV]

I have copied what you gave in the question...try using a single &.

[Samuel Leeney]

Thanks, but still no luck. It ran before (slowly) but doesn't run in the way you've sent.

I haven't included how 'pos' is setup, perhaps this is related?

N = 1000;
seed = 1234;
%seed=input('Random number seed = ');
rng(seed)
pos=rand(3,N);

[KSSV]

Edited the answer.....If it throws error...try using single &.

[Samuel Leeney]

Thanks a lot for your help KSSV, the code above you've sent works but runs about 50% slower than the code we started with.

[KSSV]

No way.....this should be definitely faster when compared to given loops.

[Samuel Leeney]

Its about twice as slow when N=20k

>> OptimisationPractice2 Your code

Elapsed time = 67.5016

Checked match

>> OptimisationPractice

Elapsed time = 33.2087 Initial code

Checked match

[KSSV]

How did you calculate the time? Check the below:

N = 1000;
seed = 1234;
%seed=input('Random number seed = ');
rng(seed)

pos=rand(3,N);
t1 = tic ; 
s = zeros(N,N);
for n = 1:3
    for a = 1:N
        count = 0;
        for b=1:N
            off=0;
            if (pos(n,a)<=0.25 && pos(n,b) >= 0.75)
                off=1;
            elseif pos(n,a)>=0.75 && pos(n,b)<=0.25
                off=-1;
            end
            s(a,b)=s(a,b)+(pos(n,a)-pos(n,b)+off)^2.0;
            if n == 3
                s(a,b)=sqrt(s(a,b));
            end
        end
    end
end
t1 = toc(t1) ; 

t2 = tic ;
a = 1:N ;
b = 1:N ;
[a,b] = meshgrid(a,b) ; 
S = zeros(N,N) ; 
for n = 1:3
  off = zeros(N,N) ;
  P =  pos(n,:) ; 
  off(P(a)<=0.25 && P(b) >= 0.75) = +1 ; 
  off(P(a)>=0.75 && P(b)<=0.25) = -1 ; 
  S = S+(P(a)-P(b)+off).^2;
end
S = sqrt(S); 
t2 = toc(t2) ; 

Now check t1, t2 for different size inputs. My bet is always the second code will be faster. I am comparing only the first part off the code.

[Samuel Leeney]

Thanks,

I calculated the tim using the toc function too.

I have also used the 'Rune and Time' option and the vast majority of time is spent computing the code in the

for n = 1:3

loop.

The second code is not faster. You can try it yourself.

Answer 3

[Bruno Luong]

I remove your hanling of offset (not sure what is the purpose), and this is a much faster method using delaunay triangulation:

clear
N=20000;
nd=3;

pos=rand(nd,N);

tic
s = zeros(N,N);
for n = 1:3
    for a = 1:N
        count = 0;
        for b=1:N
            s(a,b)=s(a,b)+(pos(n,a)-pos(n,b))^2.0;
            if n == 3
                s(a,b)=sqrt(s(a,b));
            end
        end
    end
end
match = zeros(1,N);
for a=1:N
    mindist=1e10;
    for c=1:N
        if  (a~=c)
            mindist=min(s(a,c),mindist);
            if (mindist==s(a,c))
                match(a)=c;
            end
        end
    end
end
toc % Elapsed time is 56.414344 seconds.

Find nearest neighbour of the same set of point in 2D or 3D

% INPUT: pos is array of size (nd x N), coordinates of N points in R^nd
tic

T = delaunay(pos.');
p = nchoosek(1:size(T,2),2);
P = T(:,p);
P = reshape(P,[],2);
P = unique(sort(P,2),'rows');
P1 = P(:,1);
P2 = P(:,2);
d2 = sum((pos(:,P2)-pos(:,P1)).^2,1);
A = [P1(:), d2(:), P2(:);
     P2(:), d2(:), P1(:)];
A = sortrows(A,[1 2]);
b = [true; diff(A(:,1),1)>0];
A = A(b,:);
nn = A(:,3).';      % index of nearest neighbour
d = sqrt(A(:,2)).'; % corresponding distance

toc % Elapsed time is 0.349317 seconds.

isequal(match,nn) % 1

Now it doesn't help you for your practice, but I still pot it here for future readers who seek for fast method.

[Samuel Leeney]

Thanks, this takes 56 seconds to run whereas mine takes 20 seconds.

I will post my answer when I figure it out!

[Bruno Luong]

Well I don't run on the same computer/data than your. Not sure if you have to bother.

[Samuel Leeney]

The code you've run is faster than mine but gives the wrong answer.

[Bruno Luong]

I run it 10 times and isequal(match,nn) return TRUE. So the answer match 10 times with random points.

> for k=1:10; benchnntest; end
Elapsed time is 49.127755 seconds.
Elapsed time is 0.516644 seconds.

ans =

  logical

   1

Elapsed time is 51.189459 seconds.
Elapsed time is 0.459055 seconds.

ans =

  logical

   1

Elapsed time is 50.431960 seconds.
Elapsed time is 0.465887 seconds.

ans =

  logical

   1

Elapsed time is 50.426246 seconds.
Elapsed time is 0.454885 seconds.

ans =

  logical

   1

Elapsed time is 50.651649 seconds.
Elapsed time is 0.567084 seconds.

ans =

  logical

   1

Elapsed time is 50.889422 seconds.
Elapsed time is 0.461514 seconds.

ans =

  logical

   1

Elapsed time is 50.678441 seconds.
Elapsed time is 0.491820 seconds.

ans =

  logical

   1

Elapsed time is 50.476219 seconds.
Elapsed time is 0.451430 seconds.

ans =

  logical

   1

Elapsed time is 52.659327 seconds.
Elapsed time is 0.443114 seconds.

ans =

  logical

   1

Elapsed time is 52.004992 seconds.
Elapsed time is 0.459873 seconds.

ans =

  logical

   1

>> 

Answer 4

[Image Analyst]

MATLAB is column major order, which means that the left most indexes go faster because they are adjacent in memory. MATLAB goes down rows first, then moves over to the next column and goes down its rows. So this slow code

for row = 1 : rows
    for col = 1 : columns
        s(row, col) = whatever; % Col iterates fastest
    end
end

will (or may be) be slower than this fast code

for col = 1 : columns
    for row = 1 : rows
        s(row, col) = whatever; % row iterates fastest
    end
end

Note that, in your code, n is your left most index of your arrays, yet you had the n loop as the outer loop, which is the slowest possible to do it. If possible, see if you can move n to an inner loop. I've had luck in the past getting nested loops to speed up doing that.

[Samuel Leeney]

Thanks, that helps a lot!

I have moved n inside the other two loops and that has gained me five seconds!

I think I need to vectorise some of the loops to gain more time, but can't seem to figure this out.

[Bruno Luong]

"I think I need to vectorise some of the loops to gain more time, but can't seem to figure this out."

The entire double for-loop

for a=1:N
    mindist=1e10;
    for c=1:N
        if  (a~=c)
            mindist=min(s(a,c),mindist);
            if (mindist==s(a,c))
                match(a)=c;
            end
        end
    end
end

can be replaced by

s(1:N+1:end) = Inf;
[~,match] = min(s,[],2);
match = match.'; % row vector

[Samuel Leeney]

How about the first loop? this seems to be using up the time. I know I need to vectorise it but really dont understand how.

[Bruno Luong]

First loop

pos_r = reshape(pos,[nd 1 N]);
s = zeros(N);
for b=1:N
    posb = pos_r(:,:,b);
    off = (pos<=0.25 & posb>=0.75) - ...
          (pos>=0.75 & posb<=0.25); 
    sb = sum((pos-posb+off).^2,1);
    s(:,b) = sb(:);
end
s = sqrt(s);

[Image Analyst]

We're not sure because you removed all your comments (which I'm sure you put in like all professional programmers do) before you posted your code so we don't really know what it's doing unless we spend a lot of time trying to understand it. Time that could be eliminated if you just put the comments back into the code to explain what it's doing. Looks like you're doing some kind of thresholding/masking but not sure.

[Samuel Leeney]

Image Analyst, what I am trying to do is generate N random co ordinates in 3 dimensions and then calculate the nearest neighbor of each point. An additional point to add is that if the nearest neighbor of a point on the far left of the 'square' would also be near to a point on the right of the 'square'. My code successfully does that however it needs to run in under 10 seconds.

Bruno has very helpfully shown me multiple ways to do this, however this is a practice excersise so I really want to understand what I'm doing so I can replicate it in the future.

Normally I don't ask so many questions however we have next to no support this year as teaching is online!

Thanks all!

[Samuel Leeney]

Bruno, this works peftectly, thank you!

[Image Analyst]

Samuel, going by your description, I'd try something like this to find the closest point.

N = 1000;
xyz = rand(N, 3); % Get N randomly located points in 3-D.
for k = 1 : N
	% Get the squared distance of point k to every other point in the array.
	distancesSquared = ((xyz(k, 1) - xyz(:, 1)) .^2 + ...
		(xyz(k, 2) - xyz(:, 2)) .^2 + ...
		(xyz(k, 3) - xyz(:, 3)) .^2);
	% We don't want to consider the distance of the point to itself, so set any zeros to infinity.
	distancesSquared(distancesSquared==0) = inf;
	% Find the min value and the index of that min for the other points.
	[minDist2, index] = min(distancesSquared);
	% Print it out.
	fprintf('Point %d at (%.2f, %.2f, %.2f) is closest to point %d at (%.2f, %.2f, %.2f).\n',...
		k, xyz(k, 1), xyz(k, 2), xyz(k, 3), index, xyz(index, 1), xyz(index, 2), xyz(index, 3));
end

It prints out stuff like:

Point 1 at (0.70, 0.37, 0.03) is closest to point 312 at (0.70, 0.37, 0.04).

Point 2 at (0.09, 0.71, 0.91) is closest to point 918 at (0.10, 0.72, 0.94).

Point 3 at (0.53, 0.95, 0.47) is closest to point 50 at (0.54, 0.93, 0.43).

etc.

Point 998 at (0.03, 0.67, 0.34) is closest to point 972 at (0.05, 0.59, 0.39).

Point 999 at (0.99, 0.87, 0.16) is closest to point 592 at (0.94, 0.93, 0.16).

Point 1000 at (0.54, 0.34, 0.42) is closest to point 540 at (0.55, 0.33, 0.42).

It doesn't do the stuff about the edges of the square though.