1 Introduction↩︎

Spherical range queries are a fundamental primitive for working with spatial data. Given a set of \(n\) points in \(d\)-dimensional Euclidean space, the objective is to retrieve all points that lie within the specified radius of a query point. In low-dimensional spaces, spatial indices significantly accelerate these queries by pruning large portions of the search space. However, in higher dimensions a brute-force approach that checks the distance to every point is often the only viable option due to the curse of dimensionality [1]. One application of spatial indices is the computation of graph embeddings [2]–[5], which play an important role in machine learning tasks such as community detection, link prediction, and node classification [6]–[10]. These methods map the vertices of a graph into a continuous space such that spatial distance reflects structural similarity [11]. This is often achieved by force-directed optimization, which applies attracting forces to adjacent vertices that are too distant from each other and repelling forces to non-adjacent vertices that are too close. The attracting forces can be computed in linear time \(O(m)\) by iterating over all \(m\) edges. However, computing repelling forces requires identifying all non-adjacent vertices within a specific distance of each other. A naive brute-force implementation would check all \(O(n^2)\) pairs of vertices, which is particularly prohibitive for large sparse graphs. Alternatively, the close vertex pairs can be identified by performing a spherical range query for each vertex, making the performance of the embedding framework heavily dependent on the efficiency of the underlying spatial index.

Beyond graph embedding, range queries also play a crucial role in numerous other domains such as particle-based simulations, database systems, and computer graphics [12]–[14]. For a comprehensive overview of applications, see the book by Samet [15]. To support these applications, numerous spatial indices have been developed; see 1 for a list. The theoretical worst-case bounds also shown in 1 give the impression that these spatial indices do not yield significant improvements over the brute-force approach. However, this is highly misleading and comes from the fact that the worst-case is often too pessimistic. In practice, the different spatial indices are often successful in pruning large chunks of the search space. In the following, we give a brief overview of the different spatial indices and broadly categorize them by how they organize and prune the search space.

The first category of spatial indices is based on partitioning the geometric space. Uniform grids partition space into fixed-size cells and insert the points into them. Queries are answered by a lookup of the corresponding cell and inspecting it and adjacent cells. These yield especially fast lookups in low dimensions when the points are evenly distributed. However, they lack adaptability and in datasets with heterogeneous density, dense cells often become a bottleneck. Orthtrees (also known as Quadtrees in 2D and Octrees in 3D) [16] address this by recursively dividing the space only where necessary. While more flexible than grids, they share the same limitation, where each split generates \(2^d\) children, which is infeasible in higher dimensions.

Unlike space partitioning methods, data partitioning methods divide the point set itself. KD-trees [17], [18] split the point set along axis-aligned hyperplanes, to construct a balanced tree. To answer a query, the tree is traversed recursively, pruning subtrees that cannot contain query results. R-trees [19] group points into bounding rectangles, VP-trees [20], [21] partition by distance to a vantage point, and Balltrees [22] use bounding hyperspheres.

The SNN structure of Chen and Güttel [23] avoids hierarchical traversal entirely. It uses the singular value decomposition (SVD) to identify the direction of maximum variance in the dataset. Points are projected onto this axis and sorted. The query algorithm uses a binary search to identify a contiguous range of candidate points that is then scanned exhaustively.

This shows that there are plenty of spatial indices to choose from. As the theoretical worst-case bounds do not provide strong guidance, choosing the right data structure for the application at hand, e.g., graph embedding, can be a daunting task. While there are some average-case results that give guarantees in certain situations [24], they are rather restricted and sometimes difficult to compare. Moreover, the performance depends on the data distribution of the specific application as well as implementation details. Consequently, an empirical evaluation is essential to choose the best data structure. The only empirical study in this direction we are aware of is by Lawson, Gropp, and Lofstead [25], whose comparison is tailored for the HPC-context, focusing on orthogonal range queries in 3-dimensional space in a distributed setting. Additionally, some libraries, e.g., the KD-tree library Kiddo [26], include benchmarks. However, to the best of our knowledge, there has been no extensive scientific comparison of spatial indices on a diverse set of inputs with varying dimensions yet.

1.1 Our Contribution↩︎

We have two main contributions. First, we address the need for an empirical comparison by providing a comprehensive evaluation of different state-of-the-art spatial indices. Secondly, we provide a new KD-tree variant that we call Sorted-Projection Radius KD-tree (SPRK-tree; pronounced “spark”), which outperforms the competition on almost all benchmarks. In the following, we discuss these two contributions in more detail; starting with the SPRK-tree.

1.1.0.0.1 SPRK-Tree.

At its core, the SPRK-tree is a KD-tree variant that recursively partitions the data points until the number of points falls below a certain predefined bucket size for the leaves; see 2.1. In the leaves, we use a variant of SNN [23], i.e., we sort the remaining points with respect to a projection line; see 2.3. For the query, we utilize pruning techniques, which reduce the query radius incrementally while traversing the hierarchy [18], [27], [28]. This allows more aggressive pruning of subtrees, while still ensuring that no points in the original query ball are omitted; see 2.2. For high-dimensional data sets, we initially rotate the data based on a singular value decomposition of sampled points, such that the first splits done by the KD-tree separate along directions with high spread; see 2.4. Our Rust implementation features SIMD vectorization and cache-aware memory layout (see 2.5), and provides bindings for C and Python.

1.1.0.0.2 Experimental Evaluation.

Our evaluation uses a large and diverse set of benchmarks, including a newly created benchmark set. Our new set comes from the graph embedding application discussed above. It features over 240 embedding instances with up to 1 M points in dimensions \(\{2,3,\ldots,16,32\}\). Beyond that, we also consider points distributed uniformly at random and benchmarks from the literature, including data sets from clustering, point of interest queries, and a high-dimensional real-world data set with point sets of up to 960 dimensions. Besides our SPRK-tree, we evaluate 10 other implementations as well as a brute-force approach. The evaluated data structures are diverse (Quadtree, uniform grid, Balltree, VP-tree, R-tree, KD-tree, SNN) and include state-of-the-art implementations of KD-trees in multiple programming languages. Our findings can be summarized as follows.

• For the embedding benchmark, KD-trees are the method of choice. While the KD-tree libraries Kiddo [26] and Neighbourhood [29] already perform well, our SPRK-tree consistently outperforms them. For sufficiently high dimensions, brute-force eventually becomes the best option, yet our SPRK-tree still outperforms brute-force by an order of magnitude for 16 dimensions.

• For randomly distributed points, tree-based methods perform well for low dimensions but degrade for higher dimensions. This is contrasted by SNN and brute-force, which perform much worse on lower dimensions but outperform the tree-based methods for dimension 32. The exception to this is our SPRK-tree. While being the best tree-based method for low dimensions, it is also (slightly) faster than SNN and brute-force for dimension 32, outperforming the other tree-based methods by an order of magnitude.

• While we developed the SPRK-tree with the embedding application in mind, it still outperforms the competitors on most benchmarks from other applications. The only exception is the benchmark of fixed-radius search in (very) high dimensions, where the SPRK-tree is usually second after SNN or brute-force, but still stays competitive.

Beyond the extensive comparison with other data structures, we provide an empirical analysis of how our design and optimization choices contribute to the performance of the SPRK-tree.

2 The SPRK-tree↩︎

In this section we introduce the SPRK-tree, which improves upon the standard KD-tree through three primary mechanisms: aggressive pruning of subtrees using radius reduction, the use of SNN-inspired sorted leaf buckets, and an SVD-based rotation to improve performance on high-dimensional data. We begin with a brief description of the standard KD-tree, followed by a detailed explanation of our improvements. We conclude by describing implementation optimizations that further boost the performance of the SPRK-tree.

2.1 KD-tree↩︎

The KD-tree [17], [18] is a spatial index that organizes a point set \(P\subseteq\mathbb{R}^d\) for range queries. It chooses an axis-aligned hyperplane to split the whole space into two half-spaces such that both half-spaces contain the same number of points. This procedure is applied recursively, yielding a binary tree, where each internal node corresponds to a hyperplane and its two children correspond to the half-spaces. When the number of points falls below a predefined bucket size, the recursion stops and creates a leaf node. Each node is associated with a region, namely the intersection of all half-spaces on the path from the root to the node. With this, the leaves form a partition of the whole space \(\mathbb{R}^d\); for an example see 1.

A degree of freedom to fill in is the choice of the axis-aligned hyperplane in each step. We achieve this by cycling through the dimensions in round-robin order and placing the hyperplane at the median coordinate of the current points along that dimension. This halves the point set at every split and produces a tree of logarithmic depth.

A spherical range query with query point \(q\in\mathbb{R}^d\) and radius \(r\) asks for all points in \(P\) within distance \(r\) of \(q\). The KD-tree resolves this by recursively traversing the tree and pruning subtrees that cannot contain query results. Specifically, at an internal node, a child subtree is pruned if the query sphere does not intersect the child’s half-space. Upon reaching a leaf, every point in the leaf is checked against the query radius by an exhaustive search.

2.2 Radius Reduction↩︎

Note that the default pruning of a KD-tree does visit unnecessary regions, e.g., Region \(2\) in 1 despite the fact that Region \(2\) does not intersect the query ball. This can be prevented by a more informed pruning procedure [18], [27], [28]. For this, consider the example in [fig:sprk1] (which is the same as in 1), and let \(p\) be the point of Region \(2\) closest to the query point \(q\). Let \(\delta_1\) and \(\delta_2\) be the distance between \(p\) and \(q\) in dimensions \(1\) and \(2\), respectively. If we had computed \(\delta_1\) and \(\delta_2\), we could have pruned early by observing that \(\|(\delta_1, \delta_2)\| > r\). To compute \(\delta_1\) and \(\delta_2\), note that they are also the distances of \(q\) to the hyperplanes \(a\) and \(b\), respectively, which both separate \(q\) from Region \(2\). This is not a coincidence: for any given region, we can compute the corresponding \(\delta_1\) and \(\delta_2\) by going up the tree; for each dimension taking the maximum distance between \(q\) and a hyperplane (of that dimension) separating \(q\) from the region.

Similar to [27], [28], we do not recompute all \(\delta\) values for each node of the tree individually but keep track of a vector \(\Delta = (\delta_1, \dots, \delta_d)\) while traversing the tree. This allows us to update \(\|\Delta\|\) using a constant number of arithmetic operations for the distance check. We start with \(\delta_i = 0\) for all \(i \in [d]\). Now let \(u\) be a node in the tree with corresponding hyperplane \(H_u\). Consider traversing from \(u\) to a child \(v\) such that \(H_u\) separates the query point \(q\) from the region corresponding to \(v\). Then we update \(\delta_i\) to the distance from the query point \(q\) to the hyperplane \(H_u\), where \(i\) is the dimension of \(H_u\).

2.3 SNN Buckets↩︎

In a regular KD-tree, the points in a leaf node are filtered by a linear scan. The SPRK-tree replaces this with an SNN data structure [23] for each leaf. The SNN structure is constructed by projecting all the points of the leaf onto an axis-aligned line and sorting along its dimension; see Region \(8\) in [fig:sprk2]. Instead of now scanning all points in the leaf, we can prune some points based on the projection. We can do this by projecting the query point \(q\) to the SNN line and then finding all points within distance \(r\) in both directions from \(q\) on the line. This can be done with a binary search followed by a linear scan along the order. Note that the resulting set of points has to be filtered for points actually in the query ball.

Observe in [fig:sprk2] that the intersection of the query ball with Region \(8\) gives rise to an improvement over this. In the direction of the SNN line, this intersection spans only a smaller range of size \(2\cdot r_{\text{SNN}}\) than the range of \(2\cdot r\) queried above. To compute this smaller radius \(r_{\text{SNN}}\), we can utilize the distance vector \(\Delta\) from the previous section. To make this precise, let the dimension of the SNN line be, without loss of generality, the first dimension. Then \(r_{\text{SNN}}\) is the maximum value such that \((r_{\text{SNN}}, \delta_2, \dots, \delta_d)\) has length at most \(r\).

Furthermore, to avoid the binary search, the SPRK-tree computes a lookup table for the sorted coordinates of the leaf nodes by discretizing the coordinate range into fixed-width bins. Pseudocode for the query procedure of a node is given in 5.

2.4 Singular Value Decomposition↩︎

In many scenarios, high-dimensional data often resides within a lower-dimensional subspace. The SPRK-tree exploits this property by rotating the data so that the \(d\) principal components with the largest variance align in descending order with the coordinate axes. Because the tree splits dimensions in a round-robin sequence, this rotation ensures that the earliest partitions occur along the directions of maximum spread, enabling more efficient pruning near the root. The rotation is computed via the SVD of a sampled subset of up to 10 k points of \(P\), which ensures that the SVD computation does not dominate the overall construction cost.

However, this transformation introduces a slight overhead during the query phase, as every subsequent query requires an initial rotation into the new coordinate space. Consequently, the SPRK-tree only applies this rotation for datasets with more than 16 dimensions. In these high-dimensional cases, the tree depth is often smaller than \(d\), meaning a standard traversal would never cycle through all dimensions. The SVD rotation guarantees that the limited number of splits are used to separate the data along the most informative dimensions.

2.5 Performance Engineering↩︎

Beyond algorithmic improvements, we engineered our implementation to utilize the full capabilities of modern hardware. In the following, we describe the most important optimizations.

2.5.0.0.1 Distance Computations.

During a query, the spatial index spends a large part of the time computing distances between the query point and candidate points in the leaf nodes. A naive SIMD (Single Instruction, Multiple Data [30]) implementation might attempt to accelerate this by computing the differences across the \(d\) dimensions of a single point simultaneously. However, this can be wasteful since the dimensionality \(d\) varies per dataset and does not necessarily align with the fixed width of hardware SIMD registers.

To resolve this, we vectorize across points rather than dimensions. The leaf buckets already store points in a contiguous array, and we group them into blocks of fixed width \(W\). We represent each block as a matrix \(M \in \mathbb{R}^{d\times W}\) in row-major order. By choosing \(W = 8\), the matrix rows fit exactly into one 256-bit SIMD register using 32-bit floating-point arithmetic.

Let \(q\) be the query point. To compute the squared distance to all \(W\) points simultaneously, we broadcast its \(i\)-th coordinate \(q_i\) into a SIMD register \(Q_i\), creating \(W\) copies of that value. Second, we calculate the coordinate-wise differences \(D_i = Q_i - M_i\) using one SIMD subtraction per dimension, where \(M_i\) is the \(i\)-th row of the matrix. We then square and sum these differences element-wise into an accumulator vector \(\mathrm{acc}\) using one fused multiply-add (FMA) instruction per dimension, effectively computing \(\mathrm{acc} = \sum_{i=1}^d D_i^2\) element-wise. This requires \(2d\) SIMD instructions (\(d\) subtractions and \(d\) FMAs) to evaluate the distance to \(W\) vectors; see 3.

For higher dimensions, we reduce the instruction count further by rewriting the squared distance \(\|p - q\|^2\) between a candidate point \(p\) and the query point \(q\) as \((p^\top p + q^\top q) - 2p^\top q\) and precomputing the squared norms \(p^\top p\), as suggested by Chen and Güttel [23]. By storing the halved squared norms, we also avoid a multiplication by \(2\) during the accumulation phase: \[\|p - q\|^2 = (p^\top p + q^\top q) - 2p^\top q \le r^2 \Longleftrightarrow \left(\frac{p^\top p}{2} + \frac{q^\top q}{2}\right) - p^\top q \le \frac{r^2}{2}.\] Computing the left side of the inequality for all \(W\) points simultaneously requires only \(d\) FMA instructions. However, it requires loading the precomputed values \(\frac{p^\top p}{2}\) from memory. Due to this overhead, it only becomes faster in practice for \(d \geq 6\).

To break data dependencies and maximize instruction throughput, we split the computation across 2 independent accumulators. We initialize them to the SIMD vectors containing \(\frac{p^\top p}{2}\) and \(\frac{q^\top q}{2}\) respectively, and alternate the subtraction of the coordinate-wise products of \(p^\top q\) using the previously allocated SIMD registers \(M_i\) and \(Q_i\). For \(d > 30\), utilizing 4 accumulation registers provides additional performance gains by further hiding instruction latency; see the detailed analysis in 3.3. Note that the dot-product reformulation potentially introduces numerical instabilities. We therefore add a small fixed \(\varepsilon=10^{-4}\) to the query radius in our implementation. Across all our test cases, this returned every point found by a naive reference implementation.

After computing all candidate distances, the indices corresponding to points fulfilling \(\|p - q\| \le r\) must be compacted into the result buffer. On CPUs with AVX-512F [31], [32] support, we use the vcompressps [33] instruction to gather the indices in a single step.

2.5.0.0.2 Memory Layout.

In addition to vectorized arithmetic, we optimize for cache-efficiency. Inner nodes are extremely lightweight, storing only the 32-bit floating-point split value of the dividing hyperplane. In contrast, leaf nodes are heavy, containing the SNN bucket with metadata and point range indices. Therefore, we keep leaves and inner nodes in separate arrays, and use for the inner nodes an implicit tree layout analogous to binary heaps.

Beyond the memory layout, we also optimize memory access patterns. During a query, a traditional KD-tree implementation frequently switches between traversing the tree and processing distance computations in the leaf nodes. This causes cache misses, as the CPU evicts the tree structure from the cache to load point coordinates. To resolve this, we decouple the query into a two-phase algorithm. The first phase traverses the tree in depth-first search order to compute the candidate ranges for the leaves. It identifies all ranges that intersect the query ball and collects them into an intermediate buffer. The second phase then computes the actual distances. Because the coordinate data is stored contiguously, iterating over the collected ranges results in a mostly linear access pattern (with skips between ranges), thereby improving the memory prefetcher performance.¹

3 Experiments↩︎

In this section, we present our experiments comparing the SPRK-tree against state-of-the-art spatial data structures. 3.1 outlines the experimental setup, introduces the baseline data structures, and details our new graph embedding benchmark. Following this, we compare query performance across a diverse set of benchmarks in 3.2, showing that the SPRK-tree consistently achieves the fastest query times across most evaluated datasets. Finally, we identify which optimizations contribute most to its efficiency in 3.3.

3.1 Experiment Setup↩︎

All benchmarks are performed single-threaded on a dual-socket Intel Xeon Gold 6144 CPU @ 3.50 GHz with 8 cores per socket, 32   L1 and 1   L2 cache per core, and 24.75   L3 cache per socket, supporting AVX-512. The system has 192   DDR4 RAM and runs Ubuntu 24.04 LTS. Rust code was compiled with rustc version 1.94.0 in release mode with RUSTFLAGS = -Ctarget-cpu=native. C code was compiled with gcc version 14.2 using -O3 and -march=native. For the Python libraries, CPython v3.12.3 was used.

After a 3-second warm-up, each benchmark runs for at least 20 seconds and at least 10 repetitions. We additionally verified that the measured times are consistent with hardware performance counters (cpu_cycles, cpu_ref_cycles, and instructions).

3.1.0.0.1 Evaluated Data Structures.

We compare the SPRK-tree against a broad selection of state-of-the-art spatial index implementations. The evaluated implementations span the scikit-learn Balltree [36], a VP-tree implemented in Rust [37], the Boost R-tree [38], and SNN [23]. Since KD-trees proved to be the strongest category overall, we include KD-tree implementations across multiple languages: the Rust crates Kiddo [26], nabo [39] and Neighbourhood [29], the C libraries nanoflann [40] and CGAL [28], and the scikit-learn KD-tree [36].

In addition to our SPRK-tree, we provide a brute-force baseline and our own Orthtree and uniform grid implementations. Notably, our grid implementation requires the query radius to be specified at construction time, which may be prohibitive for some applications. For datasets with varying query radius, we build the grid using the average query radius. To exclude performance differences resulting from the choice of implementation language, we also provide our own Rust reimplementation of SNN, evaluated in 11.

3.1.0.0.2 Graph Embedding Benchmark.

As introduced in 1, we provide a new benchmark set based on graph embeddings. We use geometric inhomogeneous random graphs (GIRGs) [41] generated using the method in [42]. GIRGs are a widely used model for real-world networks with heterogeneous degree distributions and community structure. For generation, we set the average degree to 15, the power-law exponent to 2.5, the dimension to 4, and the temperature to 0 (corresponding to \(\alpha = \infty\)). We vary the number of vertices \(n\) such that \(n\in\{10^4, 10^{4.5}, 10^5, 10^{5.5}, 10^6 \}\) and use 3 different seeds for each parameter configuration.

To obtain a data set from these graphs, we compute embeddings using the WEmbed algorithm [3]. The algorithm iteratively adjusts vertex positions starting from a random initial state. In each iteration, the algorithm performs \(n\) spherical range queries, finding the set of vertices that are within a certain distance for each vertex. The query radius depends on the vertex degree, with high-degree vertices having a larger query radius. To obtain data sets with varying dimension, we set the embedding dimension \(d\) from 2–16 and additionally include 32. This results in \(5 \cdot 3 \cdot 16 = 240\) different embeddings across all graph sizes.

The embedding algorithm updates the vertex positions for \(t\) iterations until convergence (averaging \(t\approx 530\)), yielding \(t\) different position and query sets per embedding. For each run, we choose the positions from the latest iteration divisible by 100 to include in our benchmark, as later iterations feature less random and more interesting vertex positions and are usually more expensive. To normalize the number of queries, we choose at most 10 k queries per embedding, sorting the vertices by degree and then choosing every \(k\)-th vertex for \(k = n / 10^4\).

We also tested other parameter configurations during generation, but observed similar performance trends across all configurations; see 6. For readability, we only show results for \(n\in\{10^4, 10^5, 10^6\}\) in our evaluation. For benchmarks spanning all vertex counts see 7.

3.2 Performance Comparison↩︎

In this section, we compare the query performance of the evaluated data structures. We begin with the graph embedding benchmark, which motivated the design of the SPRK-tree. To evaluate the SPRK-tree’s applicability beyond graph embedding and obtain a more robust comparison, we consider four additional categories: uniform point distributions, high-dimensional real-world data, clustering workloads, and Point of Interest (POI) queries. For these, we also provide unaggregated results in 10.

3.2.0.0.1 Graph Embedding.

4 illustrates query times across embedding dimensions and vertex counts for a subset of data structures, including the best performing implementation for each algorithm. The SPRK-tree achieves the best query time in every tested configuration, across all vertex counts and dimensions. For \(n=\SI{10}{k}\) and \(d \ge 12\), all structures except the SPRK-tree are outperformed by brute-force, an expected consequence of the curse of dimensionality. Here, the SPRK-tree likely prevails due to better-optimized constant factors. SNN generally exhibits similar performance to brute-force for high dimensions. We observe that SNN performs worse in dimensions not divisible by \(4\), which is likely due to the BLAS [43] implementation using a \(4\times4\) kernel for the sgemv subroutine. Furthermore, at 1 M vertices several structures (Kiddo, Boost R-tree, Neighbourhood, and scikit-learn Balltree) are slower in 2 dimensions than in 4. This might be an artifact of the embedding benchmark, where lower dimensions often have a high density of points just outside the query radius, thereby slowing structures that are unable to efficiently prune these points.

2 lists the comprehensive results for all evaluated structures. Most competitors not shown in 4 are clearly outperformed by the remaining structures. In higher dimensions the grid and Orthtree ran out of memory (oom). For construction times see 8.

3.2.0.0.2 Uniform Distributions.

Next, we evaluate a benchmark of uniformly distributed points within a unit hypercube. We vary the number of points from 10 k to 10 M and the dimension from 2 to 32. For each point set, we perform 10 k queries with the query radius set such that the expected number of returned points is 15. 5 shows the results, restricted to the best-performing implementations for readability. The SPRK-tree is the fastest structure in every tested configuration. Its advantage is most pronounced at low point counts and low dimensions, where it outperforms the next-best competitor by up to an order of magnitude. This gap narrows at high dimensions with many points, maintaining a factor of roughly \(\num{1.6}\times\) over brute-force at dimension 32 with 10 M points. Most tree-based methods scale similarly to the SPRK-tree, but for low point counts the SPRK-tree is significantly faster, likely due to better cache performance; see 9 for additional details.

The relative ranking of the competitors shifts with dimensionality. In low dimensions (2 and 8), Kiddo and the Boost R-Tree are close to the SPRK-tree, while SNN is significantly slower. For larger dimensions, the more sophisticated tree-based structures tend to be increasingly outperformed by brute-force, and SNN converges toward brute-force performance.

3.2.0.0.3 High-Dimensional Real-World Datasets.

We further evaluate data structures on high-dimensional real-world datasets from the benchmark collection of Aumüller, Bernhardsson and Faithfull [44]. This collection features datasets with up to 10 M points and up to 960 dimensions; see 3 (first section). Some datasets exhibit very high dimension-to-point ratios; the small SIFT dataset, for instance, has 128 dimensions but only 25 k points.

Based on our previous benchmark results, one would expect tree-based methods to struggle heavily in this regime. While other tree-based implementations indeed fall behind by an order of magnitude or more, the SPRK-tree remains competitive, converging toward brute-force and SNN performance levels. SNN is the fastest method on four of the six nearest neighbor datasets, while the SPRK-tree consistently places second. We attribute the SPRK-tree’s strong performance mostly to the initial SVD rotation and optimized distance computation. The SVD is also applied by SNN, explaining the similar performance.

3.2.0.0.4 Clustering.

We also evaluate the clustering benchmarks introduced in the SNN paper [23]. These benchmarks use DBSCAN [47], a density-based clustering algorithm whose core operation is a fixed-radius range query for every point in the dataset. We isolate the all-to-all query component to measure spatial index performance independently of the clustering logic (note that [23] reports the end-to-end runtime). The datasets are small, containing 200–4500 points, but span dimensions from 4–256.

As shown in the second section of 3, the SPRK-tree achieves the best query time on all clustering datasets, often by a wide margin. Brute-force and Kiddo alternate as the second-best method depending on the dimensionality: Kiddo places second on the lower-dimensional datasets (Banknote, Ecoli), while brute-force takes second on the higher-dimensional ones (Dermatology, Wine). SNN performs poorly across all clustering datasets, likely due to constant overheads that dominate at these small point counts.

3.2.0.0.5 POI Queries.

To evaluate performance on geographic data, we benchmark fixed-radius POI queries, e.g., finding all bus stops within 500 m of a given location. These datasets are two-dimensional but exhibit highly heterogeneous density, with dense clusters in cities and sparse coverage in rural areas. We use data from the OpenStreetMap project for Germany [46], with five categories: ATMs (\(\approx \SI{12}{k}\)), bakeries (\(\approx \SI{33}{k}\)), parking facilities (\(\approx \SI{769}{k}\)), bus stops (\(\approx \SI{761}{k}\)), and restaurants (\(\approx \SI{103}{k}\)). Each query set is constructed by using POIs of one category as dataset and POIs of a different category as query points. Queries are run at radii of 500 m, 1000 m, 2000 m, and 5000 m, and times are averaged across all radii and queries.

As shown in the POI section of 3, the SPRK-tree achieves the best query time across all five datasets. SNN is 2 to 3 orders of magnitude slower in this setting. The reason is likely that in 2 dimensions, the pruning of tree-based methods is more effective than projecting onto a single axis. Among tree-based methods, the Orthtree is the strongest competitor, yet the SPRK-tree outperforms it by a factor of \(2.8 \times\) on the smallest dataset (ATM, 12 k points) and up to \(3.3 \times\) on larger datasets (Bus Stop, 761 k points). We attribute this to spherical pruning, which reduces the distance checks by roughly 50 % in 2 dimensions.

3.3 Breakdown of Performance Gains↩︎

In this section, we analyze the individual contributions of the different optimizations (both algorithmic and low-level) to the overall performance of the SPRK-tree. We first measure the reduction in distance checks achieved by radius reduction and SNN buckets, followed by evaluating the SIMD and memory layout optimizations.

3.3.0.0.1 Algorithmic Optimizations.

Recall that the SPRK-tree is a KD-tree with radius reduction and SNN in the leaves (Section 2). In 6, we present an ablation study for these optimizations on the 100 k embedding dataset. We measure the number of distance computations performed by the data structure, which is independent of low-level optimizations.

The baseline is our KD-tree implementation without radius reduction or SNN (KD-tree) as described in 2.1. Using only radius reduction (Radius Reduction) cuts the number of distance checks significantly, achieving up to 50 % fewer checks in higher dimensions. On the other hand, SNN buckets in isolation (SNN Buckets) are most effective in lower dimensions, where they achieve a similar reduction of 50 %. The full benefit is realized by combining both techniques (Radius + SNN) and also using the reduced radius to further narrow the SNN candidate range (representing the full SPRK-Tree). This brings the total number of distance checks down to about 25 % of a regular KD-tree for \(d \ge 6\).

3.3.0.0.2 Distance Computation.

To evaluate the low-level optimizations of the distance computation introduced in 2.5, we conduct a microbenchmark. We measure the computation of 128 distances to a single query point across varying dimensions. We compare the direct distance computation without optimizations (naive), the approach using fused multiply-add instructions (FMA), the inner product reformulation (SNN Acc), and the variants using 2 and 4 accumulation registers (SNN 2 Acc and SNN 4 Acc). The results are shown in 7.

The fused multiply-add approach provides moderate improvement over the naive baseline, reducing runtime by roughly 10 % on average. The inner product reformulation initially performs worse, but becomes faster than FMA at 6 dimensions. In low dimensions, loading the precomputed norms from memory likely outweighs the benefits of reduced arithmetic operations. The variants with multiple accumulators provide additional improvements in higher dimensions. For the dimensions targeted by our main benchmarks (\(d \leq 16\)), 2 accumulators offer a good trade-off, with a speedup of up to 35 % over one accumulator. This is likely due to fewer data dependencies, allowing the CPU to better hide instruction latency.

Based on these insights, the SPRK-tree uses a dimension-dependent case distinction. For \(d < 6\) the SPRK-tree uses the FMA-based distance computation. For dimensions up to \(d = 32\), it switches to the inner product method with 2 accumulators. Beyond that, it uses 4 accumulators. We note that the optimal number of accumulators is architecture-dependent, as it interacts with register pressure and the compiler’s auto-vectorization logic.

3.3.0.0.3 Cumulative Optimization Effects.

In 8, we compare the full SPRK-tree against a naive variant that retains all algorithmic optimizations but disables the low-level optimizations from 2.5, including the memory layout and access pattern optimizations. The gains are most pronounced for low point counts, where for example even writing results to the output buffer constitutes a significant fraction of the total query time. Notably, the naive SPRK-tree already performs well for small point counts and high dimensions. This demonstrates the effectiveness of the algorithmic optimizations alone, with the low-level engineering further improving the performance by up to an order of magnitude at small point counts.

4 Conclusion↩︎

While developed primarily for repelling-force computation in graph embeddings, the SPRK-tree proves to be an efficient, general-purpose spatial index across a wide range of applications. Even in extreme scenarios, such as high-dimensional real-world datasets, the SPRK-tree remains competitive with the best-performing methods. This is achieved with a combination of algorithmic optimizations, such as the spherical pruning and SNN buckets, and low-level engineering that maximizes hardware utilization. Integrating the SPRK-tree should allow to substantially improve the performance of applications across various domains.

Another common strategy for high-dimensional regimes is to resort to approximate queries. Our initial attempts in this direction were rather unpromising, due to a significant drop in embedding quality for only small performance gain. However, additional work is warranted to fully evaluate the potential of this direction.

5 SPRK-Tree Query Procedure↩︎

In 9 and 10, we provide the pseudocode for the SPRK-tree query procedure, which is described in 2. Since the SPRK-tree is balanced and every branch has the same depth, whether a node is a leaf or an internal node can be determined by its depth.

6 Graph Embedding Parameters↩︎

In 3, we described our new graph embedding benchmark set. In 11, we provide our justification for focusing on a single parameter configuration in the main benchmark. We benchmarked the effect of varying the GIRG generation parameters (average degree, power-law exponent, temperature, and sampling dimension) and observed similar performance trends across all configurations. Most notably, the change in relative performance of the different data structures remains largely within the margin of error across all configurations.

7 Additional Graph Embedding Results↩︎

In 4 we evaluated the structures on a subset of the vertex counts of our benchmarking set to highlight the interesting performance characteristics. In 12 we show the results for all vertex counts for dimensions 2, 8, and 16. The analyzed trends remain consistent with the subset shown in 4, with the SPRK-tree achieving the best query time across all configurations.

8 Index Construction Time↩︎

In the context of graph embedding, the index needs to be reconstructed from scratch in every iteration of the embedding algorithm, making the index construction time a relevant performance metric. The index construction time of the benchmarked data structures, on the other hand, is negligible compared to the query time of the \(n\) queries performed in each iteration, so we do not include it in the main benchmark. 13 shows the index construction time of all evaluated data structures for the embedding benchmark. The implementations were tested in their default configurations, some of which use multithreading by default. The SPRK-tree has a competitive construction time, staying within the same order of magnitude as the other fast-constructing, tree-based methods. The Orthtree exhibits a noteworthy construction time with a steep increase with the number of dimensions but periodic large drops. The steepness naturally correlates to the exponential increase in the number of cells with increasing dimension. The periodic drops on the other hand are likely due to the fact that the depth of the tree decreases. Both the Orthtree and the uniform grid run out of memory for larger point counts at higher dimensions.

9 Effect of Low-Level Optimizations for the uniform distribution benchmark↩︎

In 3 we highlighted the performance gains of the SPRK-tree on the uniform distribution benchmark, especially at low point counts and low dimensions. We attribute this to better cache performance, which is a result of the low-level optimizations described in 2.5. In 14 we demonstrate the effect of these optimizations for this benchmark by comparing the fully optimized SPRK-tree to a naive variant that retains all algorithmic optimizations but disables the low-level optimizations from 2.5. We observe that the optimizations provide a significant performance boost for lower point counts in 2D, emphasizing that the slope of the fully optimized SPRK-tree is due to better cache performance rather than worse scaling.

This can also be explained mathematically by considering the cache occupancy of the dataset. For example, in dimension 8, each point occupies 8 times 32 bits, so 10 k points (312  ) fit comfortably in the L2 cache (1  ). From 10 k to 100 k points the data starts to outgrow L2, causing a steep initial rise, exceeding its capacity at 100 k points (3.1  ). Beyond this range the scaling flattens again: although the full dataset surpasses the 24.75   L3 cache, the tree-traversal metadata remains substantially smaller and stays resident in L3. Because cache occupancy also depends on the dimension, this transition shifts to lower point counts as dimension increases, producing steeper initial slopes at higher dimensions.

10 Unaggregated Results for the Real-World Benchmarks.↩︎

In 3 we reported the average query time across all queries and multiple radii for the high-dimensional real-world datasets, clustering, and POI benchmarks. In 4, 5, and 6 we provide the unaggregated results. In contrast to the aggregated results in 3, 6 includes the uniform grid instead of the ball tree since it provides additional interesting insights. Keep in mind that the uniform grid profits strongly from the fact that the radius is fixed within each benchmark, providing a somewhat optimistic lower bound on its performance in realistic scenarios. Nonetheless, it is noteworthy that the uniform grid is still beaten by the SPRK-tree in almost all POI benchmark configurations, only surpassing it for small radii on the smaller datasets.

11 Comparison of SNN Implementations↩︎

Throughout the paper we used the Python implementation of SNN provided by the original paper [23]. In 7 we compare this implementation to our own Rust implementation of SNN, which utilizes some of the low-level optimizations described in 2.5 such as the accumulation register optimization from 7. We observe that our implementation is significantly faster than the original Python implementation for low dimensions but is beaten by the original implementation in very high dimensions. This is likely due to the original implementation using a highly optimized BLAS library for the matrix-vector multiplications. Additionally, our implementation always processes at least 8 points per iteration, while the vectorized approach of the original can be more efficient when the candidate set is very small.

References↩︎

---

1. As our access pattern is already very efficient, more complex inner node layouts yielded no further benefits in our experiments. For example, we observed no improvement with a van Emde Boas layout [34], [35].↩︎