How do I choose an optimal cell size when searching for close pairs of points, and using cells to implement this?

Question

Suppose that I have a set of $N$ points in $k$-dimensional space ($k>1$), such as in this question, and that I need to find all pairs with a distance¹ smaller than a certain threshold $t$. The brute-force method would require $N(N-1)$ distance calculations, which is not acceptable. I attack the problem by first sorting the cells in a grid¹, such as in this answer, followed by brute-force within each grid cell and a number of neighbours (which is easily calculated from the cell size $w*h$ and the maximum distance $t$).

My solution seems to work acceptably well for my purposes, and the results appear to be correct. However I'm neither a computer scientist nor a mathematician, and I'm not sure what tools I could use to calculate the optimal cell size. In fact, I developed the aforementioned possibly naive algorithm because it seemed like a reasonably okay method. I guess the optimal cell size depends in some way on $N$, $t$, on the cost of the distance function, and on the implementation of the sorting in cells, on the distribution of points, and on other things. How would I make a guess of the optimal values of $w$ and $h$, with or without a priori knowledge on the approximate number of pairs I expect to find?

Does the answer change if the N points are divided in two sets $S_1$ and $S_2$, and each pair shall consist of one element from each set?

---

¹^{Not necessarily euclidian. The points may, for example, be locations on a sphere, i.e. on Earth, with latitude and longitude.}

Answer 1

There's an enormous amount of work on data structures and algorithms for this sort of problem. I suggest you start by reading the general literature on this problem.

Start by reading about algorithms for nearest neighbor search, including all nearest neighbors and the fixed-radius nearest neighbors problem. Those techniques are applicable to your problem.

Then, read about quadtrees, octrees, k-d trees, VP trees, and the general category of BSP trees.

Answer 2

Let's start with a special case, $k=2$, and a bunch of assumptions. Let's assume that your grid is a dense $k$ dimensional square array of cells, and each cell contains some relatively simple implementation of a set of points (for example, a linked list of points.) Let's also stick with Euclidean distances, assume you don't care about the edge conditions (points that move outside the grid), and that the points are independently uniformly randomly distributed.

I'm also going to assume you know the average density (in terms of expected number of points within a sphere of radius $t$).

Now let's look at the two extreme cases. If you make your cells too small then most of the cells will be empty, and you'll spend most of your time iterating through empty cells. If you make your cells too large then you will do a larger number of extra distance computations for points that are farther than $t$ away.

I can't think of any reason you would want your cells to be larger than $t \times t$. With $k=2$ and cells of this size you need search 9 cells. That's a search area of $(3t)^k$. With cells of size $t/2 \times t/2$ you need to search 25 cells. That's a search area of $\frac{25}{4}t^2$. With cells of size $t/M \times t/M$ you need to search approximately $(2M+1)^k$ cells, and the area you search is $(\frac{2M+1}{M}t)^k$, so in the limit you are getting arbitrarily close to $(2t)^k$. (When $M$ gets sufficiently large you may be able to avoid searching some of the cells in the corners, but I'm ignoring that too.)

This is where the density, $\delta$, (the average number of points in a volume of $t^k$) comes in. If the cost of iterating over an empty cell is similar to the cost of doing a distance calculation, then you are trying to minimize $$(2M+1)^k + \delta (\frac{2M+1}{M})^k$$. Take the derivative with respect to $M$, set it equal to 0, and solve for $M$ and you get something like $\hat{M}^{k+1} = \delta/2.$