Nearest neighbours problem

I was interested if somebody has any experience with quick finding closest neighbours for given points.

Lets say I have list of N points with x, y and z coordinates and want to find the closest neighbours for all points. The simplest way would be:

for(i=1;N;i++){
for(j=1;N;i++){
r(i,j)=(x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2
if(r(i,j)<rmin){
rmin=r(i,j)
x1=x(i)
x2=x(j)
y1=y(i)
.....
}
}
}

and so on... But the problem here is that time is increasing as O(N^2). I have found lots of talks over the web, that time can be reduced to O(N log N) what can be done by using various trees, sorting, etc. but I can't find anywhere programming code in any languages how this think can be treated. If anyone find or know something about this (not just Csharp but even in other languages, like C++, Fortran, etc.), it would be a great pleasure!

 

If I get you correctly, you have a group of points of which you want to determine which two are closest to each other?

If so, I would have a global list of the points. I would have each point access this list and compare points to itself. I would then report back each of these distances to a manager. The manager will then sort thru to find the smallest distance and which two points it involved.

Your approach seems kind if cluttered.

 

Hello
I think you could save some O by not checking every point against every point. I.e. if you have checked the first point against every other point there is no reason to check any more points against that first point. Then you have checked the second point against all other points and there is no reason to ..... etc. I changed your code. I did not check the O or the code.

Code (csharp):

1.  
2. for(i=1;N;i++){
3.  for(j=i;N;j++){ // Changed 1 to i here
4.  r(i,j)=(x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2
5.  if(r(i,j)<rmin){
6.  rmin=r(i,j)
7.  x1=x(i)
8.  x2=x(j)
9.  y1=y(i)
10.  .....
11.  }
12.  }
13.  }
14.  

Patrik

Edit: I did some calculations. I find that it cuts the number of calculations in half but is still O(N^2).
Edit Again: I see now that you wanted closest neighbour for each particle. A little different than the problem I solved. I still have some cool code in my next post that might at least give something to work off.

 

I was a little disappointed so I started matlab and invented something.
Here is the code for the O(N^2) algorithm followed by my invented algorithm.
For N = 100000 the O(N^2) code takes 162 seconds and the improved code takes 0.2 seconds.

Code (csharp):

1.  
2. clear
3. n = 100000;
4. x = rand(n,3);
5. minD2 = 1;
6. k = 0;
7. tic
8. for i = 1:n
9.     for j = i+1:n
10.         d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2;
11.         if d < minD2
12.             minD2 = d;
13.             r = [i,j];
14.         end
15.         k = k + 1;
16.     end
17. end
18. toc
19. minD2
20. r
21. k
22.  
23. minD2 = 1;
24. k = 0;
25. tic
26. x(1:n,4) = 1:n;
27. x = sortrows(x);
28. x(n+1,1) = 2;
29. for i = 1:n
30.     j = i + 1;
31.     while (x(i,1)-x(j,1))^2 < minD2
32.         d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2;
33.         if d < minD2
34.             minD2 = d;
35.             r = [x(i,4),x(j,4)];
36.        end
37.         k = k + 1;
38.         j = j + 1;
39.     end
40. end
41. toc
42. minD2
43. r
44. k
45.  

The output was

Code (csharp):

1.  
2. >> NearestNeighbour
3. Elapsed time is 161.684106 seconds.
4.  
5. minD2 =
6.  
7.    9.3041e-08
8.  
9.  
10. r =
11.  
12.        60246       96101
13.  
14.  
15. k =
16.  
17.    5.0000e+09
18.  
19. Elapsed time is 0.198043 seconds.
20.  
21. minD2 =
22.  
23.    9.3041e-08
24.  
25.  
26. r =
27.  
28.        60246       96101
29.  
30.  
31. k =
32.  
33.      3517366
34.  

It depends on you having an efficient sortrows algorithm. I did some numeric calculations on N and time and found something close to O(N*logN). Pretty cool.

Patrik

Edit: This code finds the shortest distance between any two particles. In my next post I find all closest neighbours.

 

This will find all closest neighbours of all particles.
For N = 100000 the
O(N^2) code is 293 seconds and
improved code is 13 seconds.

I shall do my own stuff now

Code (csharp):

1.  
2. clear
3. n = 10000;
4. x = rand(n,3);
5. k = 0;
6. tic
7. for i = 1:n
8.     minD2 = 2;
9.     for j = 1:n
10.         if i ~= j
11.             d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2;
12.             if d < minD2
13.                 minD2 = d;
14.                 r(i,1) = j;
15.                 r(i,2) = minD2;
16.             end
17.             k = k + 1;
18.         end
19.     end
20. end
21. toc
22. r1 = r;
23. k
24.  
25. clear r
26. k = 0;
27. tic
28. x(1:n,4) = 1:n;
29. x = sortrows(x);
30. x = [-4 0 0 0; x];
31. x(n+2,1) = 4;
32. for i = 2:n+1
33.     j = i + 1;
34.     minD2 = 3;
35.     while (x(i,1)-x(j,1))^2 < minD2
36.         d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2;
37.         if d < minD2
38.             minD2 = d;
39.             r(x(i,4),1) = x(j,4);
40.             r(x(i,4),2) = minD2;
41.        end
42.         k = k + 1;
43.         j = j + 1;
44.     end
45.     j = i - 1;
46.     while (x(i,1)-x(j,1))^2 < minD2
47.         d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2;
48.         if d < minD2
49.             minD2 = d;
50.             r(x(i,4),1) = x(j,4);
51.             r(x(i,4),2) = minD2;
52.        end
53.         k = k + 1;
54.         j = j - 1;
55.     end
56. end
57. toc
58. r2 = r;
59. k
60.  
61. % [r1 r2]
62.  

 

Thanks very much Kavorka for sharing your knowledge! It looks great. By the way, I was wondering if you have complete code of this script (with include, use and other statements - if there are such ones), which I could quickly copy and run to test with other numbers (lets say N=7000,8000,9000,10000,11000,12000) to check time increment profile over N? And which compiler do you use - that I could easily check it?

I am usually using heapsort, which I tested many times and it works as O(N log N). Examples can be used from here:
http://rosettacode.org/wiki/Sorting_algorithms/Heapsort#C.23

This sort is actually unstable, but it means that it can be not sorting values which are exactly the same. But as my coordinates are always a bit different, heapsort seems to be the best way. Which sorting is used by this sortrows(x) - I see it should be sorting by all three coordinates, right?

 

Hello
No it sorts only on the first coordinate.
That is the complete script. Matlab is an interpreted language so there is no compiler. The sortrows() is some built in function and I have no idea of its algorithm.
You told me I could use any language so I used the one I always use for calculations. Matlab has extended support for manipulating vectors like I do on rows 27-30. I am sure you can figure out a way to port it to your language of choise. I don't have to declare variables in matlab, they are created by assigning to them.

 

Oh, ok, I will try to translate it to my language and will check out how sortrows() actually works. So if this is correct the example of coordinates:

3 8 1
1 4 9
2 1 3

should be sorted like that:
1 4 9
2 1 3
3 8 1

is that correct?

 

Yes.

 

well, problem seems to be almost solved. By the way, as we are still fresh memory on this question, what you think would be the most efficient way in this code to find m-nearest neighbours (for example 2, 3-nearest neighbours for each point)?

 

I just translated this code to Fortran 90 and plotted time(N) law for Kavorka algorithm.

program main
integer :: nn,n,m,k
double precision, dimension), allocatable :: x,y,z,r
integer, dimension), allocatable :: index,ir

double precision :: d, minD2, minD2a


double precision :: t1, t2

integer :: useMod

useMod = 1

nn=200

open(1,file="rez.dat",status="replace")
open(2,file="rez2.dat",status="replace")

write(1,*) '#index ','jj ','x ','y ','z ','x2 ','y2 ','z2 ','ir ','r '

write(2,*) '#n ','t '

do jj=1,nn

print *, jj

n = 100*jj

print *, '1'

allocate(index(n+2))
allocate(x(n+2))
allocate(y(n+2))
allocate(z(n+2))

allocate(ir(n+2))
allocate(r(n+2))


print *, '2'

call cpu_time(t1)

index(1)=1
x(1)=-10.0
y(1)=0.0
z(1)=0.0

index(n+2)=n+2
x(n+2)=10.0
y(n+2)=0.0
z(n+2)=0.0

print *, '3'

do i =2,n+1

index(i)=i
x(i)=rand()
y(i)=rand()
z(i)=rand()



end do

print *, '4'

if(useMod.EQ.0) then

do i =2,n+1
minD2a = 10.0
do j =2,n+1
if(i.NE.j) then
da = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2

if(da.LT.minD2a) then
minD2a = da
ir(i) = j
r(i) = minD2a

end if

end if
end do
end do




end if

print *, '5'





if(useMod.EQ.1) then
! call sort(n,x)
print *, '6'
call heapsort(x,y,z,index)
print *, '7'
i=0
j=0
k=0

do i = 1,n+2
index(i)=i
end do

print *, '8'

do i =2,n+1

j = i+1
minD2 = 10.0

do while(((x(i)-x(j))**2).LT.minD2)
if(((x(i)-x(j))**2).LT.minD2) then
d = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2

if (d .LT. minD2) then
minD2 = d
ir(index(i)) = index(j)
r (index(i)) = minD2
end if

k = k + 1
j = j + 1
end if
end do



j=i-1
do while(((x(i)-x(j))**2).LT.minD2)
if(((x(i)-x(j))**2).LT.minD2) then
d = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2

if (d .LT. minD2) then
minD2 = d
ir(index(i)) = index(j)
r (index(i)) = minD2
end if

k = k + 1
j = j - 1
end if
end do

end do

print *, '9'



end if




print *, '10'

do i =2,n+1
! write(1,*) index(i),jj,x(i), y(i), z(i),x(ir(i)), y(ir(i)), z(ir(i)), ir(i), r(i)
end do
print *, '11'
call cpu_time(t2)


write(2,*) n, t2-t1
print *, '12'
deallocate(index)
deallocate(x)
deallocate(y)
deallocate(z)

deallocate(ir)
deallocate(r)

print *, '13'

end do

close(1)
close(2)













contains

! sorting along double precision with integer is in addition
subroutine heapsort(a,b,c,ia)

double precision, intent(in out) :: a(0
double precision, intent(in out) :: b(0
double precision, intent(in out) :: c(0
integer, intent(in out) :: ia(0
integer :: start, n, bottom
double precision :: atemp,btemp,ctemp
integer :: itemp

n = size(a)
do start = (n - 2) / 2, 0, -1
call siftdown(a,b,c, ia, start, n);
end do

do bottom = n - 1, 1, -1
atemp = a(0)
btemp = b(0)
ctemp = c(0)
itemp = ia(0)

a(0) = a(bottom)
b(0) = b(bottom)
c(0) = c(bottom)
ia(0) = ia(bottom)

a(bottom) = atemp;
b(bottom) = btemp;
c(bottom) = ctemp;
ia(bottom) = itemp;
call siftdown(a,b,c, ia, 0, bottom)
end do

end subroutine heapsort

subroutine siftdown(a,b,c, ia, start, bottom)

double precision, intent(in out) :: a(0
double precision, intent(in out) :: b(0
double precision, intent(in out) :: c(0
integer, intent(in out) :: ia(0
integer, intent(in) :: start, bottom
integer :: child, root
double precision :: atemp,btemp,ctemp
integer :: itemp

root = start
do while(root*2 + 1 < bottom)
child = root * 2 + 1

if (child + 1 < bottom) then
if (a(child) < a(child+1)) child = child + 1
end if

if (a(root) < a(child)) then
atemp = a(child)
btemp = b(child)
ctemp = c(child)
itemp = ia(child)

a(child) = a (root)
b(child) = b (root)
c(child) = c (root)
ia(child) = ia (root)

a(root) = atemp
b(root) = btemp
c(root) = ctemp
ia(root) = itemp

root = child
else
return
end if
end do

end subroutine siftdown








end program main

 

Wow, I am a little impressed. I wrote the original code and I still have a hard time reading that Fortan. I was a little worried that the part where I extend the size of the vector would give you some trouble, but you apparently had none. Result is not quite O(NlogN), is it? But at least it is a lot better than O(N^2), see picture. I can't judge, are you somewhat pleased or disappointed.

 

If the algorithm is not the fastest available I excuse myself with the fact that I developed it in an hour. I noticed you changed some things in the code to suit your ensemble, that really tells me you understand what you are doing.

This will find the m nearest neighbours of all n particles and their relative distances.

Code (csharp):

1.  
2. clear
3. n = 10000;
4. m = 10; % m < n
5. x = rand(n,3);
6.  
7. clear r
8. k = 0;
9. tic
10. x(1:n,4) = 1:n;
11. x = sortrows(x);
12. x = [-4 0 0 0; x];
13. x(n+2,1) = 4;
14. for i = 2:n+1
15.     j = i + 1;
16.     rtemp = ones(m,2) * 3;
17.     b = 3;
18.     while (x(i,1)-x(j,1))^2 < b
19.         d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2;
20.         if d < b
21.             rtemp = [rtemp(1:m,1:2); d x(j,4)];
22.             rtemp = sortrows(rtemp);
23.             b = max(rtemp(1:m,1));
24.         end
25.         k = k + 1;
26.         j = j + 1;
27.     end
28.     r(x(i,4),1:m,1:2) = rtemp(1:m,1:2);
29.     j = i - 1;
30.     while (x(i,1)-x(j,1))^2 < b
31.         d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2;
32.         if d < b
33.             rtemp = [rtemp(1:m,1:2); d x(j,4)];
34.             rtemp = sortrows(rtemp);
35.             b = max(rtemp(1:m,1));
36.         end
37.         k = k + 1;
38.         j = j - 1;
39.     end
40.     r(x(i,4),1:m,1:2) = rtemp(1:m,1:2);
41. end
42. toc
43. % r
44. k % k is not essential to the script. It is number of particles checked.
45.  

Checking for 5 nearest neighbours of 6 particles yields:

Code (csharp):

1.  
2. >> NearestNeighbours
3. Elapsed time is 0.007401 seconds.
4.  
5. r(:,:,1) =
6.  
7.     0.2305    0.4006    0.4595    0.7971    1.2752
8.     0.2305    0.4034    0.4647    0.9743    1.0744
9.     0.2411    0.3653    0.4034    0.7971    1.0028
10.     0.3653    0.5359    0.7394    1.0744    1.2752
11.     0.2411    0.4006    0.4127    0.4647    0.5359
12.     0.4127    0.4595    0.7394    0.9743    1.0028
13.  
14.  
15. r(:,:,2) =
16.  
17.      2     5     6     3     4
18.      1     3     5     6     4
19.      5     4     2     1     6
20.      3     5     6     2     1
21.      3     1     6     2     4
22.      5     1     4     2     3
23.  
24.  
25. k =
26.  
27.     30
28.  

Patrik

 

I calculated the time law (analytically), it is O(N^1.67).

 

I am pretty happy with results and here I have a plot which explains everything. Unimproved algorithm seems to be working on O(N^2), while improved algorithm could be working on O(NlogN), while sorting algorithm could be working on O(N).

If your law would be correct they could be working on O(N^2), O(N^1.67) and O(NlogN), so which version is correct?

I am also curious what law we will get if m is getting closer to n? Do we come back to O(N^2) or not?

 

well, it seems like you are right and this algorithm works with O(N^1.67), as I plotted const*x*log10(x) together with sort law on this plot.

 

Intreresting work. What will you use this for? NlogN would have been nice but this works, and is fairly easy to implement. Maybe you could improve it furter or find the NlogN algorithm. I have assumed that m << n but it works also for m = n-1 but I don't know at what O. I guess it is not good for large m.

 

I found some stunning results - performance in this algorithm depends on how many space dimensions we use (I just commented z component when calculating d in 2D space (x,y plane) and y component in 1D space). It seems like laws are the following:

3D O(N^1.67)
2D O(N^1.33(?))
1D O(N^1.00(?))

so in 1D we really have O(NlogN) (not O(N)) if we account sorting. So if we repeat 1D case for y and z coordinates separately, maybe we can get something even more close to O(NlogN). But I am not sure if it would be possible to combine these results for final r. Do you have any ideas about that?

Here is that updated plot again for 2D and 1D.




Oh, I am creating RTS Battle system, which looks for nearest neighbour (NN) to attack:
http://forum.unity3d.com/threads/223334-3-9k-units-battle-system
it currently involves the 'inefficient' i-j looping to get NN. I also included skipping j loops frequency based i-NN from previous r calculations. But everything got dependent on NN problem. So if we would manage to get O(NlogN) - it would be endless number of units used in my system compared to what it is now.

 

Cool. I was just thinking that if my calculations are correct it would be O(N^1.5) in 2D. I calculated (again not by book but by thinking about it) that it should be O(N^(2-1/d)) where d is number of dimensions.

I cannot think of an easy way to combine three 1D runs. But then again nothing I say here have been published in Nature so consider everuthing I say as potentially wrong.

RTS is fairly flat? Then 2D could be good enough? I am also feeling something about a flat 3D problem maybe could be better than O(N^1.67).

Post if you find more interesting results and/or make improvements.

P.S. In 1D you just need sort plus checking two scalars, right? Did you try it for 1D? Is that actual times?

 

Saw your system, cool. Now that I have seen the application I may have something. Since you have more information than I assumed, how about this:
Get the m closest once.
Next time look for closest neighbours among your closest neihbours and their closest neighbours.
Statistically, with a well chosen m, it will be correct most of the time.
This, I beleive could be about O(NM^2). Which for fixed m is O(N).
That would be something! Am I right, is it doable?

 

Edit: Back to thinking it will work and be fantastic.

 

well, I still have a hope to try here so called k-d tree which is binary search tree and is searching in k-dimensional space:
http://en.wikipedia.org/wiki/Kd-tree

quite good implementation seems to be implemented in C++ here:
https://code.google.com/p/kdtree/downloads/list

I still didn't find out how to run this code and what it gives. Another problem is that it could be not so easy to translate it from C++ to C#, as Unity uses C#.

Yes, RTS is pretty flat and I can use 2D in case if there would be no other options. However as I am building 3D game, I would start experiencing problems in mountain areas. I am also thinking about using these algorithms in collision detecting problems, which is now working I guess as O(N^2). I don't like to apply current collision detection also because it applies force along r if r is close enough for NN instead of keeping NN at the constant r if critical value of r is reached. So there are several cases where I would apply NN problem.

well, I assuming that in 1D the only change in your code would be that instead of
d = (x(i,1)-x(j,1))^2 + (x(i,2)-x(j,2))^2 + (x(i,3)-x(j,3))^2;
you put
d = (x(i,1)-x(j,1))^2;
everything else remaining the same.