Searching for FixedRadius Near Neighbors with Cell Lists Algorithm in Julia Language
Table of Contents
Introduction
This article defines the fixedradius near neighbors search problem and why efficient algorithms are essential. We present two algorithms for solving the problem, Brute Force and Cell Lists. We explain both algorithms indepth and analyze their performance. We accompanied this article with a Julia implementation of the Cell Lists algorithm, located in CellLists.jl repository. We can use the algorithm in a simulation or as a reference for implementing it in another language. Cell Lists is applicable in simulations with finite range interaction potentials where we are required to compute all close pairs in each time step, such as particle simulations. For longrange interaction potentials, we recommend reading about BarnesHut simulation and Fast multipole method.
Fixedradius near neighbors is an old computer science problem closely related to sorting. We define the problem as finding all pairs of points within a certain distance apart, referred to as all close pairs. The author of ^{1} has surveyed different methods for solving the fixedradius near neighbors problem. These methods include brute force, projection that is, dimensionwise sorting and querying using binary search, cell techniques, and kd trees. The paper ^{2} analyzes the brute force, projection, and cell techniques. In the context of particle simulations, ^{3} discusses Verlet Lists, a method with similarities to Cell Lists. These methods have also been applied to large scale particle simulations using GPU by ^{4}.
Similar to sorting, the computational complexities of the different methods in respect to the number of particles $n$ are $O(n^2)$ for the Brute Force approach, $O(n\log n)$ for comparison sort based approaches such as projection and kd trees, and $O(n)$ for bucket sort based approaches such as cell techniques. However, the $O$notation hides the constant that depends exponentially on the dimension $d$ for comparison and bucketbased approaches. Nevertheless, in low dimensional spaces, these methods can be significantly more efficient than the Brute Force.
Regarding analyzing algorithms, we recommend Introduction to Algorithms ^{5} textbook. Foundations, sorting, hash tables, and computational geometry sections are relevant to this article. The computational geometry section contains a chapter (33.4) about finding the closest pair of points in two dimensions, which is related to finding all close pairs.
FixedRadius Near Neighbors Problem
Let $I=\{1,…,n\}$ denote the indices of $nââ$ points, $r>0$ denote a fixed radius, and $ð©_iââ^d$ denote the coordinates of the point $iâI$ in $dââ$ dimensional Euclidean space and norm as $\â \$. The near neighbors of point $i$ are the distinct points $jâI, jâ i$ within radius $r,$ denoted as
Using the definition of near neighbors, we define all close pairs as
Because the norm is symmetric, that is $\ð©_ið©_j\=\ð©_jð©_i\,$ the order of $i$ and $j$ does not matter. We can simplify the all close pairs to symmetric all close pairs where we use the set $\{i, j\}$ to denote the inclusion of either $(i, j)$ or $(j, i),$ defined as
The next sections discuss how to implement algorithms for fixedradius near neighbors search.
Brute Force Algorithm
The Brute Force is the simplest algorithm for computing the near neighbors. It serves as a benchmark and building block for the Cell Lists algorithm, and it is useful for testing the correctness of more complex algorithms.
The Brute Force algorithm directly follows the definition symmetric all close pairs. It iterates over the points pairs, checks the distance condition, and yields all pairs that satisfy the condition. The following Julia code demonstrates the Brute Force approach.
@inline function distance_condition(p1::Vector{T}, p2::Vector{T}, r::T) where T <: AbstractFloat
sum(@. (p1  p2)^2) â€ r^2
end
function brute_force(p::Array{T, 2}, r::T) where T <: AbstractFloat
ps = Vector{Tuple{Int, Int}}()
n, d = size(p)
for i in 1:(n1)
for j in (i+1):n
if @inbounds distance_condition(p[i, :], p[j, :], r)
push!(ps, (i, j))
end
end
end
return ps
end
n, d = 10, 2
r = 0.1
p = rand(n, d)
brute_force(p, r)
The Brute Force algorithm performs a total number of iterations:
$$B(n) = â_{i=1}^{n1} i = \frac{n(n  1)}{2} \in O(n^2). \tag{4}$$
Since Brute Force has quadratic complexity, it does not scale well for many points, motivating the search of more efficient algorithms.
Cell Lists Algorithm
Definition
The Cell Lists algorithm works by partitioning the space into a grid with cell size $r$ and then associates each cell to the points that belong to it. The partitioning guarantees that all points within distance $r$ are inside a cell or its neighboring cells. We refer to the set of all such point pairs as cell list pairs $C.$ It contains all close pairs by the relation $\hat{N}_r â C â \hat{N}_{R},$ where $R=2r\sqrt{d}$ is the maximum distance between two points in neighboring cells. Therefore, we can find all close pairs by iterating over cell list pairs and checking the distance condition.
We implement the Cell Lists algorithm in Julia language. The algorithm has two parts. The first one constructs the Cell List, and the second one iterates over it to produce the symmetric all close pairs.
Constructor
The CellList
constructor maps each cell to the indices of points in that cell.
struct CellList{d}
indices::Dict{CartesianIndex{d}, Vector{Int}}
end
function CellList(p::Array{T, 2}, r::T, offset::Int=0) where T <: AbstractFloat
@assert r > 0
n, d = size(p)
cells = @. Int(fld(p, r))
data = Dict{CartesianIndex{d}, Vector{Int}}()
for j in 1:n
@inbounds cell = CartesianIndex(cells[j, :]...)
if haskey(data, cell)
@inbounds push!(data[cell], j + offset)
else
data[cell] = [j + offset]
end
end
CellList{d}(data)
end
Near Neighbors
We define the indices to the neighboring cells as follows:
@inline function neighbors(d::Int)
n = CartesianIndices((fill(1:1, d)...,))
return n[1:fld(length(n), 2)]
end
We only include half of the neighboring cells to avoid double counting when iterating over the neighboring cells of all nonempty cells. We also define functions for bruteforcing over a set of points.
@inline function distance_condition(p1::Vector{T}, p2::Vector{T}, r::T) where T <: AbstractFloat
sum(@. (p1  p2)^2) â€ r^2
end
@inline function brute_force!(ps::Vector{Tuple{Int, Int}}, is::Vector{Int}, p::Array{T, 2}, r::T) where T <: AbstractFloat
for (k, i) in enumerate(is[1:(end1)])
for j in is[(k+1):end]
if @inbounds distance_condition(p[i, :], p[j, :], r)
push!(ps, (i, j))
end
end
end
end
@inline function brute_force!(ps::Vector{Tuple{Int, Int}}, is::Vector{Int}, js::Vector{Int}, p::Array{T, 2}, r::T) where T <: AbstractFloat
for i in is
for j in js
if @inbounds distance_condition(p[i, :], p[j, :], r)
push!(ps, (i, j))
end
end
end
end
Finding all cell list pairs works by iterating over all nonempty cells and finding the pairs within that cell using Brute Force. Iterating over its neighboring cells checks if the cell is nonempty and finds the pairs between the current and the nonempty neighboring cell using Brute Force.
function near_neighbors(c::CellList{d}, p::Array{T, 2}, r::T) where d where T <: AbstractFloat
ps = Vector{Tuple{Int, Int}}()
offsets = neighbors(d)
# Iterate over nonempty cells
for (cell, is) in c.data
# Pairs of points within the cell
brute_force!(ps, is, p, r)
# Pairs of points with nonempty neighboring cells
for offset in offsets
neigh_cell = cell + offset
if haskey(c.data, neigh_cell)
@inbounds js = c.data[neigh_cell]
brute_force!(ps, is, js, p, r)
end
end
end
return ps
end
Usage
We can use the Cell Lists by first constructing the CellList
and then query the symmetric all close pairs with near_neighbors
.
n, d = 10, 2
r = 0.1
p = rand(n, d)
c = CellList(p, r)
indices = near_neighbors(c, p, r)
We receive the same output as for brute force.
indices2 = brute_force(p, r)
@assert Set(Set.(indices)) == Set(Set.(indices2))
In the next section, we analyze the Cell Lists algorithm’s performance compared to Brute Force and theoretical limits and discuss the assumptions when Cell Lists is practical.
Analysis
Number of Iterations
We denote the number of nonempty cells using
Each nonempty cell has $n_1,n_2,…,n_{\hat{n}}$ points such that that the total number of points is $n=n_1+n_2+…+n_{\hat{n}}.$ Then, the number of iterations to find cell list pairs is
where $N(i)$ is half of the nonempty neighboring cells of cell $i$ where $N(\hat{n})=â .$ We can iterate over the neighboring cell pairs (highlighted in red) in two ways.
Method 1: For each nonempty cell $iâ\{1,…,\hat{n}1\},$ iterate over the remaining nonempty cells $j=\{i,…,\hat{n}\}$ and check if they are a neighbor, which takes $B(\hat{n})$ iterations.
Method 2: For each nonempty cell $iâ\{1,…,\hat{n}1\},$ iterate over neighboring cells and check if they are nonempty (as in near_neighbors
code), which takes $\hat{n}â
m_d$ iterations where $m_d=(3^d1)/2$ is the number of neighboring cells in dimension $d.$
We derive the dimensionality condition such that method 2 performs fewer iterations compared to method 1:
When true, finding neighboring cells grows linearly, necessary for the Cell Lists to be efficient compared to Brute Force. Note that higherdimensional spaces require exponentially higher numbers of nonempty cells $\hat{n}$ to satisfy the condition.
WorstCase Instances
If one cell contains all the points, we have $\hat{n}=1,$ and $n_{1}=n,$ the number of iterations equal to Brute Force
$$M(1)=\textcolor{darkblue}{B(n_1)}+\textcolor{darkred}{0}=B(n). \tag{8}$$
If all points are in unique cells and we do not satisfy the dimensionality condition, we have $\hat{n}=n$ and $n_1,n_2,…n_{\hat{n}}=1$ and $\hat{n}â€3^d,$ the number of iterations equals to Brute Force
$$M(n)=\textcolor{darkblue}{\hat{n}}+\textcolor{darkred}{B(\hat{n})}=B(n)+n. \tag{9}$$
Curse of Dimensionality
Cell Lists also suffers from the curse of dimensionality. As we increase the dimension, the number of possible cells grows exponentially. That is,
$$â_{i=1}^d \left(\frac{l_i}{r}\right) â¥ \left(\frac{l}{r}\right)^d, \tag{10}$$
where $l_1,l_2,…,l_{d}$ are the lengths of the minimum bounding box and $\min\{l_1,l_2,…,l_{d}\}=l>r$ is the minimum length. Therefore, if we draw two random cells from a nonexponential probability distribution where possible cells are its states, the probability that they are the same cell decreases exponentially as the dimension grows. Consequently, the number of nonempty cells approaches the number of points $\hat{n}ân,$ and the dimensionality condition requires larger $\hat{n}$ to satisfy. Thus the performance approaches the worstcase.
Practical Instances
In many practical instances, we can assume a limit to how many points can occupy a cell. For example, two and threedimensional physical simulations particles have finite size and nonoverlapping and satisfy the dimensionality condition.
For a problem instance that satisfies the dimensionality condition $\hat{n}>3^d,$ and the number of points per cell is limited by $1â€n_1,n_2,…,n_{\hat{n}}â€\bar{n}$ where $\bar{n}$ is independent of $n.$ By assuming that every cell is filled $n_1,n_2,…,n_{\hat{n}}=\bar{n}$, we have $n=\hat{n}â \bar{n}$, and we obtain an upper bound for the number of iterations
With fixed dimension and maximum density, the number of iterations is linearly proportional to the number of points $O(n)$. The constant term can be large; therefore, we should always measure the Cell Lists' performance against the Brute Force on real machines with realistic problem instances.
Benchmarks
Method
We performed each benchmark instance by setting fixed parameter values for size $n$, dimension $d$, and radius $r,$ and generating $n$ points $ð©=(p_1,…,p_d)$ where we draw the coordinate values from a uniform distribution, that is, $p_iâŒU(0, 1)$ for all $i=1,…,d.$ We measured the performance of each instance using BenchmarkTools.jl and selected the median time. We performed the benchmarks on the Aalto Triton cluster running on Intel(R) Xeon(R) CPU E52680 v4 @ 2.40GHz processors.
Cell Lists vs. Brute Force
The figures above show how Cell Lists compares to Brute Force as we increase dimension with a constant radius. As can be seen, Cell Lists performs better in low dimensions and beat Brute Force with a small number of points, but requires more points with a larger dimension.
Effects of Dimensionality
The figure above shows the relation between Cell Lists with increasing dimension and constant radius. As we can see, the increase in time scales roughly to $3^d.$
Conclusions
I first encountered the fixedradius near neighbors problem while developing the Crowd Dynamics simulation, which is mechanistically similar to molecular simulation with finite distance interactions (local interactions). The algorithms presented in this article can be used for research purposes or building small scale simulations. Large scale simulations may require additional improvements such as parallelism, which we discuss more in Multithreading in Julia Language Applied to Cell Lists Algorithm.
Contribute
If you enjoyed or found benefit from this article, it would help me share it with other people who might be interested. If you have feedback, questions, or ideas related to the article, you can write to my GitHub Discussions forum.
***
For more content, you can follow my YouTube channel and join my newsletter. Since creating content and opensource libraries take time and effort, consider supporting the effort by subscribing or giving a onetime donation.
References

Bentley, J. L. (1975). A Survey of techniques for fixed radius near neighbor searching, 37. ↩︎

Mount, D. (2005). Geometric Basics and FixedRadius Near Neighbors, 1â8. ↩︎

Yao, Z., Wang, J. S., Liu, G. R., & Cheng, M. (2004). Improved neighbor list algorithm in molecular simulations using cell decomposition and data sorting method. Computer Physics Communications, 161(1â2), 27â35. https://doi.org/10.1016/j.cpc.2004.04.004 ↩︎

Hoetzlein, R. (2014). Fast FixedRadius Nearest Neighbors: Interactive MillionParticle Fluids. GPU Technology Conference. ↩︎

Cormen, T. H., Leiserson, C. E., Rivest, R. L., & Stein, C. (2009). Introduction to Algorithms, Third Edition (3rd ed.). The MIT Press. ↩︎