# $k$ nearest neighbors for $n$ points with periodic boundary conditions in R^3

Take $n$ points ${\bf r} =\{ r_1 = (x_1,y_1,z_1), r_2=(x_2,y_2,z_2) \ldots r_n=(x_n,y_n,z_n) \}$ enclosed in a periodic box of length $L$, such that that the point $(0,0,0)=(L,0,0)$, $(0,0,0)=(L,L,L)$, etc... The difference between two points is the minimum distance of the first point and any of the images of the second point.

Given a point $p$ inside the box, how can I find the $k$ nearest neighbors without brute force? I'm not sure how to apply a standard space partition (like a kd-tree) as I don't see how I can "split" the space.

Do you have to compute the nearest neighbours of all points, or just of a single point? And do you have to compute it/them once, or several times?

If you cannot assume anything about your points, then a spatial Octree is probably your best bet, cutting the space into eight equal cubes recursively until you have a small number of points per cell. You can look for neighbours of a point by inspecting all points in the same cell and successively climbing up the tree and inspecting neighbouring cells as you go.

If you can assume a maximum distance $r_c$ within which all $k$ neighbours will lie, then you could use a Cell List, which consists of partitioning the space into cells of edge length of at least $r_c$ in all dimensions. The neighbours of any given point will then be in any of the 26 cell surrounding that point's cell, so these are the only candidate points you need to consider. This approach is of course only useful if $r_c$ is reasonably bounded and if you are trying to find the neighbours of more than one point.

Dealing with boundary conditions is not much of an added difficulty, since in both cases you are dealing with cells of points which can be wrapped periodically.

Update

If you're only considering a single target point, you may be better off just computing the $k$ nearest brute-force.

If you have to compute the $k$ neighbours continuously, and the points only move by a bit, you can get some improvement by using a Verlet list (a.k.a. Neighbour list): Instead of finding the $k$ nearest neighbours, find the $K>k$ nearest. If $r_0$ is your target point, compute $s = \|r_0-r_K\| - \|r_0-r_k\|$, i.e. the difference in distance from the $k$th and $K$th nearest neighbours.

By doing this, in every iteration you only have to re-check these $K$ nearest neighbours and pick out the $k$ nearest. That is, until any of your points has moved more than $s$ since you first looked for the $K$ neighbours. If that happens, just re-start the whole procedure.

Update 2

You may already know this, but if you're aiming for efficiency, you should try using the Median of Medians algorithm to find the $k$ or $K$ smallest point distances.

• I would like to compute the NN of a single point, but the underlying point cloud would be slowly changing (the points are coordinates to a Monte-Carlo simulation). – Hooked May 30 '12 at 21:06
• I have updated my reply for the special case of a single target point. – Pedro May 30 '12 at 21:52

There is a simple answer where you just replicate and translate the points into your acceleration structure. Requires 8x more ram, but 1/8th the code of a real "solution." For less than 1e7 points or so this works pretty well actully reasonably.

There is a nice KD-tree implementation in scipy.spatial that makes this easy.

I used it once to calculate the Renyi parking coefficient for high-dimensional spaces. If you are interested, I can see if I can dig up code.

For a single point query, there is no beating brute-force.

• Isn't it only 8 times (32-1) more ram for a 2D simulation and 26 times (33-1) for a 3D simulation? While an 8-fold increase in memory is possible, 26 times may not be. – Hooked May 31 '12 at 18:28
• I too believe that the factors in @Hooked's comments are the correct ones. – Wolfgang Bangerth May 31 '12 at 21:06