k-d Trees for Range Search and Nearest Neighbours

A k-d tree is a binary search tree over points in k dimensions. Each internal node splits on one axis at a threshold, sending points below the threshold to the left subtree and the rest to the right, and the axis cycles through the dimensions with depth. The structure is a generalisation of the one-dimensional search tree to the case where there is no total order on the data, and it is a common choice of spatial index because it is simple, compact, and needs no parameters beyond a leaf size.

The two queries it answers are range search, which reports every point inside an axis-aligned box, and nearest-neighbour search, which finds the closest point to a query. Both use the same recursive descent with a different pruning rule, and both have an expected cost far below the linear scan on typical data, with a worst case that is exactly the linear scan.

Construction: median split, not midpoint split

The naive split places the threshold at the middle of the coordinate range, which produces a balanced-looking tree whose subtrees can be badly unbalanced in point count, because the points may all lie on one side of the midpoint. The standard construction splits at the median coordinate value, which guarantees that the two subtrees differ in size by at most one.

Splitting at the median requires either sorting the points on the split axis, which costs O(n log n) per level and O(n log² n) overall, or a linear-time median selection such as quickselect, which gives O(n log n) overall. The second is what serious implementations use. A common shortcut is to sort the points once on the first axis and then, at each level, partition the current range in place by the median of the other axis; the bookkeeping is fiddly but the total stays O(n log n).

def build(points, depth=0):
    if not points:
        return None
    axis = depth % len(points[0])
    points.sort(key=lambda p: p[axis])          # or nth_element / quickselect
    mid = len(points) // 2
    return Node(point=points[mid], axis=axis, value=points[mid][axis],
                left=build(points[:mid], depth + 1),
                right=build(points[mid + 1:], depth + 1))

The node stores the split point, the axis, and the split value. Some implementations store only the split value and keep points in the leaves, which saves space when the leaf size is larger than one; a leaf holding 8 to 32 points is the usual compromise, because the last few levels of a point-per-node tree are almost entirely pointer chasing while a small array is scanned contiguously.

Range search and the pruning argument

A range query descends with the box. At a node splitting on axis a at value v, the left subtree contains only points with coordinate a less than v and the right subtree only points with coordinate a at least v. If the query box’s upper bound on axis a is below v, the entire right subtree can be skipped; if the lower bound is above v, the left subtree can be skipped. Otherwise both are visited.

def query(node, box, out):
    if node is None:
        return
    if contains(box, node.point):
        out.append(node.point)
    lo, hi = box[node.axis]
    if lo <= node.value:                 # left subtree may intersect
        query(node.left, box, out)
    if hi >= node.value:                 # right subtree may intersect
        query(node.right, box, out)

The pruning is what separates the k-d tree from a linear scan, and its effectiveness depends on the query box’s shape relative to the tree’s splits. A box that is small in all dimensions follows a single root-to-leaf path plus a few side branches, which is O(log n). A box that is thin in one dimension and long in another is the awkward case: the splits on the long axis cannot prune anything, and the query degenerates toward a scan of the whole tree.

The classical analysis gives O(√n + k) for an axis-aligned range query in two dimensions on a balanced tree, where k is the number of reported points. The √n comes from the boundary of the query box crossing the splits: a box perimeter intersects O(√n) cells of the partition, and each contributes work. The bound assumes the tree is balanced and the query is a rectangle; for k-dimensional trees the analogous bound is O(n^(1−1/k) + k), which for k = 2 is the same √n and for k = 10 is nearly n.

Nearest neighbour and the shrinking radius

Nearest-neighbour search uses the same descent but prunes with a distance. The algorithm keeps the best candidate found so far and a radius equal to its distance from the query. At a node, the subtree on the far side of the split can be skipped if the distance from the query to the splitting plane exceeds the current radius, because every point in that subtree is at least that far away.

The radius shrinks as better candidates are found, which makes the pruning progressively stronger. The order of traversal matters for this reason: visiting the near subtree first gives a good candidate early and prunes more of the far side. Visiting the far side first is still correct but explores far more nodes.

def nearest(node, q, best):
    if node is None:
        return best
    d = dist2(node.point, q)
    if d < best.d2:
        best = Candidate(node.point, d)
    delta = q[node.axis] - node.value
    near, far = (node.left, node.right) if delta < 0 else (node.right, node.left)
    best = nearest(near, q, best)
    if delta * delta < best.d2:          # plane closer than current best
        best = nearest(far, q, best)
    return best

The expected cost for points drawn uniformly from a bounded region is O(log n), with a constant that grows with dimension. The worst case is O(n), and it is easy to construct: points on a circle, with queries at the centre, give a split plane that never prunes, because every subtree is at distance roughly equal to the radius.

Why high dimensions defeat it

The pruning argument relies on a box or ball excluding most of a subtree. In k dimensions the volume of a ball of radius r relative to a cube of side 2r falls off as the dimension grows: the ratio is π^(k/2)/(2^k Γ(k/2+1)), which is already small at k = 10 and negligible at k = 30. A nearest-neighbour ball of the typical radius therefore covers a large fraction of the bounding box, and almost no subtree can be excluded.

This is the curse of dimensionality in its index form. A k-d tree over 100-dimensional vectors degenerates toward a linear scan, and no amount of balancing changes that, because the geometry itself offers no separation. The alternatives for high dimensions are approximate methods with a recall or error guarantee, or dimensionality reduction before indexing. For genuinely low-dimensional data — 2D map coordinates, 3D mesh vertices, 4D space-time — the k-d tree remains the right default.

Variants and the practical knobs

Splitting at the median of the widest dimension, rather than cycling axes, adapts better to anisotropic data: a point cloud stretched along x gets more x-splits, and the resulting cells are closer to cubes. The cost is that each node must know its points’ extent, which requires either a pre-pass or storing bounding boxes.

Storing a bounding box per node supports a different pruning rule: skip any subtree whose box does not intersect the query box. It is strictly stronger than the split-plane test for range queries and about equally strong for nearest neighbour, at the cost of 2k extra numbers per node. With k = 2 or 3 that is cheap; with k = 50 the bounding box dwarfs the rest of the node.

Leaf buckets of size b reduce the tree depth by log₂ b and replace the bottom levels with a linear scan. For b around 16 the bucket scan touches a few 64-byte cache lines and replaces several levels of pointer dereference. The trade is that a range query must scan the whole bucket even when only one point qualifies, so the optimum depends on the query’s typical selectivity.

Dynamic updates are the weak point. Insertion into a balanced k-d tree requires either rebuilding subtrees or accepting imbalance, and deletion requires a rebuild of the affected subtree. Workloads with frequent updates usually prefer a structure designed for them, or a log-structured approach that rebuilds periodically. The static construction cost of O(n log n) is what makes the tree attractive for read-mostly data, and that is the assumption most spatial indexes in practice are built on.

Related structures split the same job: the trapezoidal map point location answers containment queries in O(log n) with a linear-size structure but only for planar subdivisions, and the Voronoi diagram answers nearest-neighbour queries in O(log n) after an O(n log n) construction. A Delaunay triangulation also supports nearest-neighbour queries, by walking from a known triangle toward the query and choosing the edge that decreases distance; that walk visits O(√n) triangles in expectation for uniformly distributed points, and covers only a few triangles when successive queries are close together. The k-d tree’s advantage is that it answers both box and ball queries with one structure and no preprocessing beyond the tree itself.