Repeated nearest neighbor calculation for millions of data points too slow

Question

I have a dataset running into millions of data points in 3D. For the calculation I am doing, I need to calculate neighbor (range search) to each data point in a radius, try to fit a function, calculate the error for the fit, repeat this for next datapoint and so on. My code works properly but it is taking very long time to run, around 1 second per datapoint! It is probably because for each point, it has to search in whole dataset. Is there a way I can make the process fast. I have an idea that if I can somehow establish some adjacency relationship between first neighbors, then this can be less slow. If it helps, I am trying to find optimum Parzen window width in 3D.

Answer 1

I'd suggesting googling for bounding volume hierarchies (BSP tree in particular). Given your point cloud, you can find a plane that splits it into two equal subclouds. Then when you need to find the collection of points that are within some radius R of a test point, you can first compare your test point to that plane, and if it's height above it is more than R, then the entire subcloud below the plane must also be farther away than R too (so you don't need to check any of those points). You can apply this idea recursively too, ultimately yielding n log n type complexities instead of n-squared. (This is BSP / binary space partitioning, but the other kinds of bounding volume hierarchies are conceptually similar - recursively clustering your points in such a way that you can quickly cull away large subcollections with a single test of your geometrical predicate).

Answer 2

There are several data structure for storing data that preserves information about position and proximity; there by allowing fast nearest neighbor(s) determination.

In particular R-trees (and specialized forms like R*-trees) and X-trees. Lots of choices that are optimized for slightly different uses.

Choosing a R*-tree rather than a naive nearest neighbor look-up was a big part of my getting a factor of 10000 speedup out of a particular code. (OK, maybe a few hundred of that was the R*-tree, most of the rest was because the naive look-up had been badly coded so that it smashed the cache. ::sigh::)

These structures have typical $O(N \log N)$ ($N$ the number of point stored) insertion performance and storage requirement and $O(\log N)$ look-up performance, so they work well if you will be doing many look ups (say one for each point as in DBSCAN); however some of them have very bad worst-case performance.

Answer 3

This is very similar to one of the biggest challenges in the field of molecular dynamics—computing all of the pairwise interactions between nonbonded particles.

There, we use cell lists (or neighbor lists) to help us figure out what's nearby; for this application, the cell list is probably the easier algorithm to use:

• Divide the box into a series of cells.
• For each particle, determine to which cell it should be assigned (O(1) per particle).
• Then, for each particle, check the "own" cell plus the neighbor cells; if any of these are occupied, then no further search is necessary.
• If all nearest neighbors are empty, then expand to next-nearest neighbors, and so on, until a particle is found.

If your system has a more or less uniform distribution of particles, this will greatly reduce the cost of your algorithm, according to the coarseness of the grid. However, some fine tuning is necessary: too coarse a grid and you won't save much time; too fine, and you'll spend a lot of time cycling over empty grid cells!

Answer 4

You should definitely check K-D trees and octrees which are the methods of choice for point sets (while BSPs are for general objects, and grids for more or less uniform densities). They can be very compact and fast, minimizing overhead in both memory and computation, and are simple to implement.

When your points are more or less uniformly distributed (even with empty areas, but there must be no density singularity or other high concentration) check sphere packings if you want to try a grid-like non-hierarchical space subdivision.

Answer 5

You should probably consider building the Delaunay triangulation (well, its 3D analogue). In 2D, that's a special triangulation of the data points that always contains the nearest neighbor. The same holds in 3D, but with tetrahedra.

You can build once and for all the triangulation, and then search for the nearest neighbor directly in the triangulation. I think that there are some good algorithms to build the triangulation: in 2D, the construction of the triangulation is in $n \log(n)$ and the later searches for the nearest neighbor is linear in the number of data points.

Hope it helps!