Approximate Nearest Neighbor Search

Approximate nearest neighbor search (ANN search, or ANNS) finds database vectors that are close to a query vector without comparing the query with all of them, at the price of sometimes returning a neighbor that is not the closest. Matching SIFT descriptors between images, assigning descriptors to visual words, and ranking millions of image embeddings by their distance to a query all reduce to nearest-neighbor search over vectors of 64 to a few thousand dimensions. Exact search compares the query with every stored vector, and in these dimensions no exact index reliably does much better. Approximate indexes, built from trees, hash tables, quantizers or graphs, return the true nearest neighbor for most queries while touching only a small fraction of the data, and some also compress the vectors so that a billion of them fit in memory.

Definition

Given a set of $n$ vectors and a distance, the nearest-neighbor problem is to preprocess the set so that, for any query, the stored vector closest to it can be found quickly; the $k$-nearest-neighbor problem asks for the $k$ closest. An approximate method may return other vectors: either any vector whose distance is within a factor of the true nearest distance (an error-bounded guarantee), or whatever the best candidate is after a fixed amount of work (a time-bounded search), judged empirically by how often it is right. An index is the data structure built in the preprocessing step. The distance is usually the Euclidean distance or the inner product; for $\ell_2$-normalized vectors, cosine similarity, inner product and Euclidean distance give the same ranking, as the similarity and distance measures article shows.

Intuition

A good index rules out most of the database without looking at it. In two dimensions this is easy: divide the plane into cells, find the query’s cell, and check that cell and its few neighbors. In 128 dimensions the same idea fails in a specific way. A query’s nearest neighbor is no longer much closer than everything else, so the ball around the query that must be searched to be sure is large, and it crosses the boundaries of a great many cells. Proving that a candidate is the nearest means visiting almost all of them.

Approximate methods give up the proof. A tree visits only the most promising few cells; a hash table looks only in the buckets the query falls into; a quantizer scans only the clusters nearest the query and compares compressed codes; a graph walks from neighbor to neighbor toward the query. Each returns the best vector it saw. Because real descriptors and embeddings lie near low-dimensional structure rather than filling their space, in practice that vector is usually the true nearest neighbor, and when it is not, it is usually almost as close.

Formal Definition

Let $X = \{x_1, \dots, x_n\} \subset \mathbb{R}^d$ be the database, $q \in \mathbb{R}^d$ a query and $D$ a distance. The exact nearest neighbor is

$$\mathrm{NN}(q) = \arg\min_{x \in X} D(q, x).$$

Two approximate versions are standard in the theory. A $(1 + \varepsilon)$-approximate nearest neighbor is any $p \in X$ with

$$D(q, p) \le (1 + \varepsilon)\, D(q, \mathrm{NN}(q)),$$

the criterion of Arya, Mount, Netanyahu, Silverman and Wu (1998). The $(c, r)$-approximate near-neighbor problem, used for locality-sensitive hashing, is a decision version: if some point lies within distance $r$ of $q$, return, with a fixed probability of success, any point within distance $cr$ (Har-Peled, Indyk & Motwani, 2012, the journal version of Indyk & Motwani, 1998). Answering it for a sequence of radii yields an approximate nearest neighbor.

Practical systems are mostly time-bounded and are scored by recall. In the definition used by ANN-Benchmarks (Aumüller, Bernhardsson & Faithfull, 2018 arXiv version), the recall of a returned list $\pi$ of $k$ vectors is the fraction that lie no farther from the query than its true $k$-th nearest neighbor $p^*_k$,

$$\mathrm{recall}(\pi) = \frac{\left|\{ p \in \pi : D(q, p) \le D(q, p^*_k) \}\right|}{k},$$

averaged over queries. The distance threshold makes ties harmless; with distinct distances it is the fraction of the true $k$ nearest neighbors found. Jégou, Douze and Schmid (2011) measured instead recall@$R$, the fraction of queries whose single true nearest neighbor is ranked within the first $R$ results, which tells how long a shortlist must be re-ranked with exact distances to recover it.

Neither is the Recall@K of retrieval evaluation. There the relevant items are images of the same object, person or class, judged by people; here the relevant item is whatever exact search would have returned, so ANN recall measures how faithfully the index reproduces brute force, not whether the embedding is any good. A perfect index over a poor embedding retrieves the wrong images perfectly.

Why Exact Search Degrades in High Dimensions

Tree indexes and the sequential scan

The k-d tree (Bentley, 1975) splits the data recursively on one coordinate at a time, cycling through the coordinates level by level. Friedman, Bentley and Finkel (1977) split instead at the median of the coordinate with the largest spread and stop at small buckets of points, so that the leaves partition space into boxes whose size adapts to the local density, and gave the branch-and-bound search, with logarithmic expected time: descend to the query’s leaf, then backtrack into every box that could hold a closer point than the best found so far. In higher dimensions, as Beis and Lowe (1997) summarized earlier experiments, the backtracking often examines a large fraction of the stored points before the nearest one is certain. Lowe (2004, §7.2) put it plainly for his 128-dimensional descriptors: the best exact algorithms, such as the k-d tree, give no speedup over exhaustive search beyond about ten dimensions.

The quantitative argument

Weber, Schek and Blott (1998) analyzed this “dimensional curse” for indexes that partition the data space or cluster the data, under the assumption of uniformly and independently distributed coordinates. They counted the disk blocks an optimal search must visit (those whose bounding regions intersect the ball around the query whose radius is the expected nearest-neighbor distance) and compared with a sequential scan, which they credited with a factor-of-five advantage from sequential reading, so that an index “fails” once it visits more than 20% of blocks. Under these assumptions the standard tree structures were outperformed by a scan from roughly ten dimensions, and for bounding spheres the 20% threshold was crossed at 26 dimensions. A cruder bound for any partitioning or clustering scheme gave 610 dimensions, which the authors stressed would in practice be well below that. Their general conclusions were that every such method has a dimension beyond which a scan wins, and that its cost tends to $O(n)$ as the dimension grows.

In their model the expected nearest-neighbor distance grows with the dimension (their Figure 3), so the ball that must be searched reaches into ever more cells; in our reading, this is related to the distance concentration described in the embeddings article. The assumption of uniform data matters, however. Muja and Lowe (2014, §4.1.1) found that on 100,000 random uniform vectors their approximate search was no faster than a linear scan beyond about 800 dimensions at 68% precision, whereas on image patches resampled to as many as 4,096 dimensions the speedup did not fall with dimension, because neighboring pixels are strongly correlated. What limits an index, in other words, is the data’s intrinsic dimension rather than the length of the vector.

Exact search remains the reference, and for small databases the right choice: a batched matrix product computes all query–database distances with optimized linear algebra. Speedups of approximate methods are quoted against it.

Families of Indexes

Trees with best-bin-first search

Best-bin-first (BBF). Beis and Lowe (1997) observed that standard backtracking visits leaves in the order of the tree, not in the order of their distance from the query. BBF keeps a priority queue of the branches not taken, keyed by the distance from the query to each branch’s region, always resumes at the closest, and stops after a fixed number $E_{\max}$ of leaves, returning the best point seen. Arya and Mount studied the same order independently as priority search, which Arya et al. (1998) analyzed for their error-bounded method. In Beis and Lowe’s experiment on 100,000 uniformly distributed points in 12 dimensions, BBF limited to 200 leaves found the exact nearest neighbor for 94% of queries, against 59% for backtracking cut off after a comparable time, while exact search examined over 2,400 leaves; up to 20 dimensions the neighbors it returned were on average only about 2% farther than the true ones. Lowe (2004) used BBF to match SIFT descriptors, checking 200 candidates against a database of 100,000 keypoints: about two orders of magnitude faster than exact search, with less than 5% loss in correct matches. He noted that BBF suits this use because the ratio test discards ambiguous matches anyway, so the hardest queries, where many neighbors are at nearly equal distances, need no exact answer.

Randomized k-d forests. A single tree misses a neighbor that lies just across a splitting plane from the query. Silpa-Anan and Hartley (2008) built several k-d trees over the same data and searched them simultaneously, and also aligned the data’s principal axes with the coordinate axes by principal component analysis before building, which improved search further. In Muja and Lowe’s implementation (2014, §3.1) each tree splits on a coordinate drawn at random from the five of highest variance, one priority queue spans all trees, and a point already compared in one tree is not compared again; on 100,000 SIFT features the speedup grew with the number of trees up to about 20 and not beyond, while memory grows linearly with it. Different random decompositions make it likely that in one of them the query and its neighbor share a cell.

Priority search k-means trees. Muja and Lowe (2009, 2014) also cluster the data recursively with k-means, using all dimensions at each node instead of one coordinate, and search the resulting tree in the same best-bin-first order. On the data sets they evaluated, one of the two, the randomized forest or the k-means tree, was always best, which one depending on the data. Their FLANN library therefore chooses the algorithm and its parameters automatically for a target precision, minimizing a cost that weighs search time against build time and memory with user-set weights (Muja & Lowe, 2009). On a million SIFT features they reported a speedup over linear search of three orders of magnitude at 60% precision with the randomized forest, and less than two orders above 90% with the k-means tree (2014, §4.1.2). The bag of visual words article describes how the same trees speed up vocabulary construction.

Locality-sensitive hashing

Locality-sensitive hashing (LSH) hashes vectors so that close ones collide more often than distant ones. Formally, a family of hash functions $\mathcal{H}$ is $(r_1, r_2, p_1, p_2)$-sensitive if, for $h$ drawn from $\mathcal{H}$,

$$D(p, q) \le r_1 \Rightarrow \Pr[h(p) = h(q)] \ge p_1, \qquad D(p, q) > r_2 \Rightarrow \Pr[h(p) = h(q)] \le p_2,$$

with $r_1 < r_2$ and $p_1 > p_2$ (Har-Peled et al., 2012, Definition 3.3). One function separates near from far only weakly, so the construction concatenates $k$ functions into a key $g(x) = (h_1(x), \dots, h_k(x))$, which makes far points rarely share a bucket, and builds $L$ independent tables, which makes a near point likely to share a bucket in at least one. A query hashes into its $L$ buckets and compares itself only with the points found there. With $\rho = \log(1/p_1) / \log(1/p_2)$, the $(c, r)$ problem is solved with $O(n^\rho)$ distance computations per query and $O(n^{1+\rho})$ space beyond the data (their Theorem 3.4); for the Hamming distance, hashing by sampling one bit gives $\rho \le 1/c$. This is the sense in which Indyk and Motwani (1998) proposed to remove the curse of dimensionality: the query time is sublinear in $n$ and polynomial, not exponential, in $d$.

Hash families exist for the common measures:

• Random hyperplanes for the angle between vectors (Charikar, 2002): $h(x) = \mathrm{sign}(a^\top x)$ with $a$ a random direction collides with probability $1 - \theta(x, y)/\pi$, where $\theta$ is the angle. Each function yields one bit, so $k$ of them form a $k$-bit code, and the fraction of differing bits estimates the angle.
• p-stable projections for $\ell_p$ distances (Datar, Immorlica, Indyk & Mirrokni, 2004): $h_{a,b}(x) = \lfloor (a^\top x + b) / r \rfloor$, where the entries of $a$ are drawn from a $p$-stable distribution (Gaussian for $\ell_2$, Cauchy for $\ell_1$) and $b$ is uniform in $[0, r]$, cuts a random projection of the data into intervals of width $r$. It works directly on Euclidean vectors, and on synthetic data the authors’ implementation was up to 40 times faster than an approximate k-d tree search.

LSH gives guarantees under few assumptions, but its fixed, data-independent partitions do not adapt to where the data lie, and the $L$ tables can take more memory than the vectors themselves (Jégou et al., 2011, §I). Muja and Lowe (2014, §4.1.2) compared the E2LSH implementation, with parameters from its own tuning tool, on 100,000 SIFT features and found their priority search k-means tree about an order of magnitude faster. Learned binary codes, such as spectral hashing, followed; Jégou et al. (2011) found product quantization more accurate than spectral hashing at the same 64-bit code length.

Inverted files and product quantization

Inverted file (IVF). A coarse quantizer, $k'$ centroids learned by k-means, assigns every database vector to its nearest centroid and stores it in that centroid’s list. A query is compared with the centroids and then with the vectors in the lists of its $w$ nearest centroids only, so about $n w / k'$ vectors are visited if the lists are balanced (Jégou et al., 2011, §IV). The neighbor is lost when it lies in a cell that was not probed, which is why $w > 1$ is needed near cell boundaries. The structure, k-means and an inverted list per cluster, is the one the bag of visual words borrowed from text search, but here the lists hold vectors to be compared rather than image identifiers to be counted.

Product quantization (PQ). To compress the vectors themselves, Jégou, Douze and Schmid (2011) split each $d$-dimensional vector into $m$ subvectors of $d/m$ dimensions and quantize each with its own k-means codebook of $k^*$ centroids. A vector becomes $m$ centroid indices, a code of $m \log_2 k^*$ bits, while the number of representable reconstructions is $(k^*)^m$. A 64-bit code for a 128-dimensional SIFT descriptor would need $2^{64}$ centroids from a single quantizer, impossible to learn or store; eight sub-quantizers of 256 centroids each give the same code length from $8 \times 256$ centroids, a choice the authors called often reasonable. Distances are estimated without decompressing:

• Asymmetric distance computation (ADC) keeps the query exact and approximates each database vector $y$ by its reconstruction $q(y)$:

$$\tilde{D}(x, y)^2 = \lVert x - q(y) \rVert^2 = \sum_{j=1}^{m} \lVert u_j(x) - q_j(u_j(y)) \rVert^2 ,$$

where $u_j$ extracts the $j$-th subvector. The $m \times k^*$ squared distances between the query’s subvectors and all sub-centroids are computed once per query, after which each database code costs $m$ table lookups and additions.

• Symmetric distance computation (SDC) quantizes the query too. It was clearly less accurate: on their SIFT data with $m = 8$, ADC with 64 centroids per subspace matched SDC with 256.

PQ recall is limited by quantization error, so a shortlist of the best $R$ codes is usually re-ranked with exact distances when the original vectors are available. IVFADC combines both ideas: an inverted file whose lists store the PQ code of each vector’s residual from its coarse centroid, which is more precise than encoding the vector itself because the residual has less energy. The authors reported coarse codebooks of $k' = 1{,}000$ to $1{,}000{,}000$ centroids for SIFT. Their Table V on a million 960-dimensional GIST descriptors (500 queries, 64-bit codes) shows the trade-off: exhaustive ADC took 17.2 ms per query for recall@100 of 0.652, while IVFADC with $k' = 1{,}024$ and $w = 8$ compared about 28,000 codes in 8.8 ms for 0.682. Against FLANN on SIFT, IVFADC with re-ranking gave better precision–time trade-offs over a large range of operating points with an index of under 25 MB against more than 250 MB, and the method was tested on about two billion SIFT vectors.

Proximity graphs and HNSW

A proximity graph links each vector to some of its near neighbors; a query is answered by greedy routing, moving from an entry vertex to whichever neighbor is closer to the query until no neighbor is closer, while keeping a list of the best vertices seen. Navigable small world (NSW) graphs (Malkov, Ponomarenko, Logvinov & Krylov, 2014) are built, as Malkov and Yashunin later described them, by inserting vectors in random order and linking each, in both directions, to the $M$ closest vectors already inserted; links made early, when the graph was sparse, become long-range shortcuts. Malkov and Yashunin also identified the weakness: both the number of greedy steps and the degree of the nodes passed grow with $n$, so search cost grows polylogarithmically, and on low-dimensional data NSW could lose to trees by orders of magnitude.

Hierarchical NSW (HNSW) (Malkov & Yashunin, 2020; the details below are from the arXiv v4 version) separates the links by length into layers, like a skip list. Each inserted vector draws a maximum layer $l = \lfloor -\ln(\mathrm{unif}(0, 1)) \cdot m_L \rfloor$, so each layer holds an exponentially smaller subset of the data than the one below; the authors recommend $m_L = 1/\ln M$. A search starts at the single entry point of the top layer, descends greedily layer by layer with a candidate list of size one, and on the bottom layer, which holds every vector, runs a best-first search with a candidate list of size $ef$, the parameter that trades speed for recall. Three construction choices matter:

• $M$, the number of links added per insertion, with at most $M_{\max}$ links per vector on the upper layers and $M_{\max 0} = 2M$ on the bottom layer; the authors found $M$ from 5 to 48 reasonable, smaller values better at low recall or low dimension and larger values at high recall or high dimension.
• $ef_{\mathrm{construction}}$, the candidate list size used while inserting, set large enough that the neighbors found during construction have a recall near 0.95 or more.
• A neighbor-selection heuristic that links a candidate only if it is closer to the new vector than to any neighbor already linked. This keeps links in several directions and preserves connectivity between clusters; the authors found it at least as good as linking the nearest candidates in every case, and most helpful for clustered or low-dimensional data.

The index costs about $(M_{\max 0} + m_L M_{\max})$ links per vector, which with 4-byte links and $M$ between 6 and 48 is about 60 to 450 bytes per vector on top of the vectors themselves. Graph construction and search are incremental, but the authors listed deletion and updates as future work; implementations added them later.

Measuring an Index

An ANN index has no single score. It is described by curves of recall against speed, usually in queries per second (QPS) or as a speedup over linear search, together with the memory used per vector and the time to build the index; every index has parameters (leaves checked, tables, probes $w$, shortlist $R$, $ef$) that move it along its curve. Comparisons are meaningful only at equal recall on the same data.

ANN-Benchmarks (Aumüller, Bernhardsson & Faithfull, 2020; the details below are from the 2018 arXiv full version) standardizes this. Each implementation runs in its own container, by default on a single CPU core in a single thread, over a grid of its parameters, on fixed data sets with precomputed exact neighbors, among them a million 128-dimensional SIFT vectors with Euclidean distance and 1.18 million 100-dimensional GloVe word vectors with angular distance (Table 3). Results are reported as Pareto frontiers of recall against QPS, and also against index size and build time. In that version’s experiments the graph-based methods were fastest on both GloVe and SIFT, HNSW at every recall level on GloVe, with KGraph close at high recall and FAISS’s inverted file third; FLANN’s automatic tuning did not finish at high recall within the five-hour limit on SIFT; and only three implementations accepted a target recall as a parameter. The benchmark’s code (commit 2e081ad, July 2026) counts a returned neighbor as correct if its distance is at most the true $k$-th distance plus $10^{-3}$, and computes QPS from the best of several timed runs (five by default).

Single-thread QPS on a fixed data set says little about batched GPU search, memory-bound billion-scale settings or data unlike the benchmark’s. Recall of true neighbors is also not task accuracy: Muja and Lowe (2014, §1) note that nearest-neighbor search is often one approximation inside a larger system that loses very little from approximate neighbors. Definitions differ between papers, too: Muja and Lowe’s “precision”, the percentage of queries for which the correct nearest neighbor is found (2009, §3), is what Jégou et al. call recall@1.

Libraries

As of October 2026:

• FAISS (Johnson, Douze & Jégou, 2021; release v1.15.1, September 2026, MIT licence) is a C++ library with Python wrappers for search with Euclidean distance or inner product, cosine being the inner product of normalized vectors. It implements exact search, inverted files, product quantization and its variants, and graph indexes (HNSW and NSG), some of them on GPUs. The accompanying paper built a $k$-nearest-neighbor graph on one billion vectors in under 12 hours on four GPUs (2017 arXiv version).
• hnswlib (latest stable release v0.9.0, March 2026, Apache 2.0 licence; v0.10.0 release candidates appeared in September 2026) is a header-only C++ implementation of HNSW from the nmslib project, with Python bindings, the spaces l2, ip and cosine, defaults $M = 16$ and $ef_{\mathrm{construction}} = 200$, and incremental insertion and deletion by marking.
• FLANN (Muja & Lowe; latest tagged release 1.9.2, April 2021, BSD licence) provides randomized k-d forests, priority search k-means trees, hierarchical clustering trees for binary descriptors and automatic configuration. Muja and Lowe (2014) note that it was incorporated into OpenCV; the repository’s last code change was in October 2022.

Examples

The script below builds a database of 10,000 local descriptors from 20 synthetic $256 \times 256$ grayscale images (intensities in $[0, 1]$) of blurred disks, rectangles and striped patches. From each image it keeps the 500 most textured $16 \times 16$ windows on a 6-pixel grid and describes each with a 128-dimensional SIFT-like vector: gradient magnitude histogrammed into 8 orientations in each of $4 \times 4$ cells, $\ell_2$-normalized, clipped at 0.2 and renormalized as in Lowe (2004), without SIFT’s detector, Gaussian weighting or interpolation. Each of 500 queries is the same window in a copy of its image with a random contrast and brightness change and fresh noise. Exact search gives the ground truth. The script first runs an exact k-d tree for 200 queries on the leading principal components of the descriptors and on the original coordinates, then compares best-bin-first search in one k-d tree and in four randomized trees, random-hyperplane LSH, an inverted file and 8-byte product-quantization codes, all written in NumPy. It reports recall@1 of the true nearest neighbor for all queries and for those that pass Lowe’s ratio test at 0.8, the mean ratio of the returned neighbor’s distance to the true nearest distance, and a cost: the number of 128-dimensional distances or dot products per query (database vectors, plus centroids, hyperplanes or the PQ table) divided by the database size. Timings in Python would mostly measure the interpreter, so the cost is counted instead. It uses NumPy 2.3 and SciPy 1.16, runs in about 20 seconds, and two runs printed identical output.

import os
os.environ["OMP_NUM_THREADS"] = os.environ["OPENBLAS_NUM_THREADS"] = "1"   # before NumPy loads
import heapq
import numpy as np
from scipy.ndimage import gaussian_filter

# ---------- data: SIFT-like descriptors of patches from synthetic images ----------
def synthetic_image(rng, size=256):
    """Grayscale image in [0, 1]: a smooth ramp plus 40 random disks, rectangles and
    stripe patches of random intensity, blurred and with pixel noise."""
    yy, xx = np.mgrid[0:size, 0:size].astype(float)
    img = 0.3 + 0.2 * (xx * rng.uniform(-1, 1) + yy * rng.uniform(-1, 1)) / size
    for _ in range(40):
        cx, cy, r = rng.uniform(0, size), rng.uniform(0, size), rng.uniform(5, 30)
        kind, v = rng.integers(3), rng.uniform(0, 1)
        if kind == 0:
            m = (xx - cx) ** 2 + (yy - cy) ** 2 < r ** 2
        elif kind == 1:
            m = (abs(xx - cx) < r) & (abs(yy - cy) < rng.uniform(0.3, 1) * r)
        else:
            t = rng.uniform(0, np.pi)
            m = ((xx - cx) ** 2 + (yy - cy) ** 2 < r ** 2) & \
                (np.sin((xx * np.cos(t) + yy * np.sin(t)) / rng.uniform(1.5, 4)) > 0)
        img[m] = v
    img = gaussian_filter(img, 1.0) + 0.01 * rng.standard_normal((size, size))
    return np.clip(img, 0, 1)

def descriptor(p):
    """128-D SIFT-like descriptor of a 16x16 patch: 4x4 cells x 8 orientation bins of
    gradient magnitude, l2-normalized, clipped at 0.2 and renormalized (as in Lowe 2004)."""
    gy, gx = np.gradient(p)
    mag, ang = np.hypot(gx, gy), np.arctan2(gy, gx) % (2 * np.pi)
    b = (ang / (2 * np.pi) * 8).astype(int) % 8
    h = np.zeros((4, 4, 8))
    for i in range(4):
        for j in range(4):
            s = np.s_[4 * i:4 * i + 4, 4 * j:4 * j + 4]
            h[i, j] = np.bincount(b[s].ravel(), mag[s].ravel(), minlength=8)
    h = h.ravel() / (np.linalg.norm(h) + 1e-12)
    h = np.minimum(h, 0.2)
    return h / (np.linalg.norm(h) + 1e-12)

def make_data(seed, n_images=20, per_image=500, n_queries=500, stride=6):
    """Database: descriptors of the 500 most textured 16x16 windows (by gradient energy, on a
    6-px grid) of each image. Queries: the same windows of randomly chosen database entries in a
    second rendering of the image with changed contrast and brightness and fresh noise."""
    rng = np.random.default_rng(seed)
    imgs = [synthetic_image(rng) for _ in range(n_images)]
    db, src = [], []
    for k, im in enumerate(imgs):
        e = np.hypot(*np.gradient(im))
        locs = [(y, x) for y in range(2, 238, stride) for x in range(2, 238, stride)]
        energy = [e[y:y + 16, x:x + 16].sum() for y, x in locs]
        for i in np.argsort(energy)[::-1][:per_image]:
            y, x = locs[i]
            db.append(descriptor(im[y:y + 16, x:x + 16])); src.append((k, y, x))
    q = []
    for i in rng.choice(len(src), n_queries, replace=False):
        k, y, x = src[i]
        view = rng.uniform(0.7, 1.3) * imgs[k] + rng.uniform(-0.1, 0.1)
        view = view + 0.02 * rng.standard_normal(view.shape)
        q.append(descriptor(view[y:y + 16, x:x + 16]))
    return np.array(db, np.float32), np.array(q, np.float32), rng

def sqdist(Q, X):
    return (Q ** 2).sum(1)[:, None] - 2 * Q @ X.T + (X ** 2).sum(1)[None, :]

# ---------- k-d tree with best-bin-first search ----------
def build_kdtree(X, idx, rng=None, leaf=8):
    """Split at the median of the highest-variance coordinate; with rng, of one drawn at random
    from the five highest-variance coordinates (randomized k-d tree as in FLANN, Muja & Lowe 2009)."""
    if len(idx) <= leaf:
        return ("leaf", idx)
    var = X[idx].var(0)
    d = int(np.argmax(var)) if rng is None else int(rng.choice(np.argsort(var)[-5:]))
    v = np.median(X[idx, d])
    left, right = idx[X[idx, d] <= v], idx[X[idx, d] > v]
    if len(right) == 0 or len(left) == 0:
        return ("leaf", idx)
    return ("node", d, v, build_kdtree(X, left, rng, leaf), build_kdtree(X, right, rng, leaf))

def kd_search(trees, X, q, max_leaves=None):
    """Best-bin-first search of one tree or a forest with a single priority queue: visit bins in
    order of their distance bound from q (Beis & Lowe 1997). max_leaves=None gives exact search.
    Returns (index of the nearest point found, number of points compared)."""
    best, best_d, n_leaves, seen = -1, np.inf, 0, set()
    heap = [(0.0, t, tree, np.zeros_like(q)) for t, tree in enumerate(trees)]
    tie = len(heap)                           # tie-breaker so the heap never compares nodes
    while heap:
        bound, _, node, off = heapq.heappop(heap)
        if bound >= best_d or (max_leaves is not None and n_leaves >= max_leaves):
            break
        while node[0] == "node":              # descend to the query's leaf, queueing far sides
            _, d, v, l, r = node
            near, far = (l, r) if q[d] <= v else (r, l)
            off2 = off.copy(); off2[d] = q[d] - v
            tie += 1
            heapq.heappush(heap, (bound - off[d] ** 2 + off2[d] ** 2, tie, far, off2))
            node = near
        pts = [p for p in node[1] if p not in seen]   # a forest shares points between trees
        seen.update(pts); n_leaves += 1
        if pts:
            dd = ((X[pts] - q) ** 2).sum(1)
            j = int(np.argmin(dd))
            if dd[j] < best_d:
                best, best_d = pts[j], dd[j]
    return best, len(seen)

# ---------- random-hyperplane LSH (Charikar 2002) ----------
def lsh_search(X, Q, n_bits, n_tables, rng):
    cand = [set() for _ in Q]
    for _ in range(n_tables):
        W = rng.standard_normal((X.shape[1], n_bits)).astype(np.float32)
        pw = 1 << np.arange(n_bits)
        kx, kq = ((X @ W) > 0) @ pw, ((Q @ W) > 0) @ pw
        buckets = {}
        for i, k in enumerate(kx):
            buckets.setdefault(k, []).append(i)
        for i, k in enumerate(kq):
            cand[i].update(buckets.get(k, []))
    return rerank(X, Q, cand)

def rerank(X, Q, cand):
    out, scanned = [], 0
    for q, c in zip(Q, cand):
        c = np.fromiter(c, int, len(c))
        scanned += len(c)
        out.append(c[np.argmin(((X[c] - q) ** 2).sum(1))] if len(c) else -1)
    return np.array(out), scanned / (len(Q) * len(X))

# ---------- k-means, IVF and product quantization (Jegou, Douze & Schmid 2011) ----------
def kmeans(X, k, rng, iters=15):
    C = X[rng.choice(len(X), k, replace=False)].copy()
    for _ in range(iters):
        a = np.argmin(sqdist(X, C), 1)
        for j in range(k):
            if np.any(a == j):
                C[j] = X[a == j].mean(0)
    return C, np.argmin(sqdist(X, C), 1)

def ivf_search(X, Q, C, assign, n_probe):
    lists = [np.flatnonzero(assign == j) for j in range(len(C))]
    probes = np.argsort(sqdist(Q, C), 1)[:, :n_probe]
    return rerank(X, Q, [np.concatenate([lists[j] for j in p]) for p in probes])

def pq_train_encode(X, m, ks, rng):
    D = X.shape[1] // m
    books, codes = [], []
    for j in range(m):
        Cj, aj = kmeans(X[:, j * D:(j + 1) * D], ks, rng, iters=10)
        books.append(Cj); codes.append(aj)
    return books, np.stack(codes, 1).astype(np.uint8)

def pq_adc(Q, books, codes):
    """Asymmetric distance: query uncompressed, database as codes; one m x ks table per query."""
    D = books[0].shape[1]
    tables = np.stack([sqdist(Q[:, j * D:(j + 1) * D], B) for j, B in enumerate(books)], 1)
    return sum(tables[:, j, codes[:, j]] for j in range(len(books)))   # (n_queries, n)

# ---------- experiment ----------
def run(seed):
    X, Q, rng = make_data(seed)
    n = len(X)
    D = np.maximum(sqdist(Q, X), 0)                      # exact search: every distance
    true_nn = np.argmin(D, 1)
    d12 = np.sqrt(np.sort(D, 1)[:, :2])
    distinct = d12[:, 0] < 0.8 * d12[:, 1]               # passes Lowe's ratio test
    rows = {}
    def add(name, found, frac, overhead=0):
        """overhead: other 128-D distances or dot products per query (centroids, hyperplanes)."""
        frac += overhead / n
        found = np.asarray(found)
        hit, ok = found == true_nn, found >= 0           # found = -1: no candidate at all
        ratio = np.sqrt(D[np.flatnonzero(ok), found[ok]]) / d12[ok, 0]
        rows[name] = (hit.mean(), hit[distinct].mean(), ratio.mean(), frac)
    tree = [build_kdtree(X, np.arange(n))]
    forest = [build_kdtree(X, np.arange(n), rng) for _ in range(4)]
    for name, trees, ml in (("k-d tree BBF, 125 leaves", tree, 125),
                            ("k-d tree BBF, 25 leaves", tree, 25),
                            ("4 random k-d trees, 125 leaves", forest, 125),
                            ("4 random k-d trees, 25 leaves", forest, 25)):
        res = [kd_search(trees, X, q, ml) for q in Q]
        add(name, [r[0] for r in res], np.mean([r[1] for r in res]) / n)
    for bits, tables in ((16, 4), (12, 8), (8, 8)):
        add(f"LSH {bits} bits x {tables} tables", *lsh_search(X, Q, bits, tables, rng),
            overhead=bits * tables)
    C, assign = kmeans(X, 128, rng)
    for w in (1, 4, 16):
        add(f"IVF 128 lists, probe {w}", *ivf_search(X, Q, C, assign, w), overhead=128)
    books, codes = pq_train_encode(X, 8, 256, rng)
    order = np.argsort(pq_adc(Q, books, codes), 1)
    for R in (1, 10, 100):                               # re-rank the top R exactly
        add(f"PQ 8-byte codes, re-rank {R}", *rerank(X, Q, [o[:R] for o in order]),
            overhead=256)                                # 8 x 256 table of 16-D distances
    return rows, distinct.mean(), n, codes.shape[1]

def kd_vs_dim(seed=0, dims=(2, 4, 8, 16, 32, 64, 128)):
    """Exact k-d tree search on the first d principal components of the descriptors (and on
    the original 128 coordinates): fraction of the database compared before the search stops."""
    X, Q, _ = make_data(seed)
    mu = X.mean(0)
    P = np.linalg.eigh(np.cov((X - mu).T))[1][:, ::-1]
    out = {}
    for d in dims + ("128, original axes",):
        if isinstance(d, str):
            Xd, Qd = X, Q[:200]
        else:
            Xd, Qd = (X - mu) @ P[:, :d], (Q[:200] - mu) @ P[:, :d]
        tree = [build_kdtree(Xd, np.arange(len(Xd)))]
        out[d] = np.mean([kd_search(tree, Xd, q)[1] for q in Qd]) / len(Xd)
    return out

if __name__ == "__main__":
    print("exact k-d tree, fraction of the database compared, seed 0, by dimension d:")
    print("  " + "  ".join(f"{d}: {f:.3f}" for d, f in kd_vs_dim().items()))
    runs = [run(s) for s in range(5)]
    _, _, n, nbytes = runs[0]
    print(f"\n{n} database descriptors (128-D float32, 512 bytes each), PQ codes {nbytes} bytes, "
          f"500 queries; fraction passing the ratio test per seed: "
          + ", ".join(f"{r[1]:.2f}" for r in runs))
    print(f"{'':32s} {'recall@1 of the true nearest neighbor':^41s}")
    print(f"{'':32s} {'seed 0 (min-max over seeds 0-4)':^41s} {'mean over seeds':^18s}")
    print(f"{'method':32s} {'all queries':^20s} {'ratio-test queries':^20s} {'dist ratio':>10s} {'cost':>7s}")
    for name in runs[0][0]:
        a = np.array([r[0][name] for r in runs])
        print(f"{name:32s} {a[0, 0]:.3f} ({a[:, 0].min():.3f}-{a[:, 0].max():.3f}) "
              f"{a[0, 1]:.3f} ({a[:, 1].min():.3f}-{a[:, 1].max():.3f}) {a[:, 2].mean():10.3f} {a[:, 3].mean():7.4f}")

exact k-d tree, fraction of the database compared, seed 0, by dimension d:
  2: 0.001  4: 0.003  8: 0.008  16: 0.017  32: 0.045  64: 0.123  128: 0.276  128, original axes: 0.741

10000 database descriptors (128-D float32, 512 bytes each), PQ codes 8 bytes, 500 queries; fraction passing the ratio test per seed: 0.70, 0.75, 0.68, 0.74, 0.71
                                   recall@1 of the true nearest neighbor  
                                      seed 0 (min-max over seeds 0-4)       mean over seeds  
method                               all queries       ratio-test queries  dist ratio    cost
k-d tree BBF, 125 leaves         0.738 (0.734-0.804) 0.763 (0.758-0.839)      1.169  0.0616
k-d tree BBF, 25 leaves          0.484 (0.474-0.578) 0.503 (0.503-0.608)      1.478  0.0123
4 random k-d trees, 125 leaves   0.848 (0.832-0.872) 0.863 (0.858-0.891)      1.094  0.0580
4 random k-d trees, 25 leaves    0.666 (0.590-0.674) 0.677 (0.639-0.691)      1.288  0.0119
LSH 16 bits x 4 tables           0.522 (0.522-0.592) 0.546 (0.546-0.626)      1.489  0.0080
LSH 12 bits x 8 tables           0.874 (0.858-0.912) 0.860 (0.860-0.923)      1.071  0.0250
LSH 8 bits x 8 tables            0.980 (0.976-0.992) 0.980 (0.980-0.992)      1.007  0.0989
IVF 128 lists, probe 1           0.924 (0.892-0.924) 0.906 (0.864-0.917)      1.091  0.0217
IVF 128 lists, probe 4           0.998 (0.996-1.000) 0.997 (0.994-1.000)      1.002  0.0453
IVF 128 lists, probe 16          1.000 (1.000-1.000) 1.000 (1.000-1.000)      1.000  0.1311
PQ 8-byte codes, re-rank 1       0.846 (0.846-0.878) 0.983 (0.981-0.989)      1.022  0.0257
PQ 8-byte codes, re-rank 10      0.986 (0.982-1.000) 1.000 (1.000-1.000)      1.000  0.0266
PQ 8-byte codes, re-rank 100     1.000 (1.000-1.000) 1.000 (1.000-1.000)      1.000  0.0356

The first line shows the curse of dimensionality on image data. An exact k-d tree over the first 2 to 16 principal components compares under 2% of the database, but the fraction grows steadily with $d$, and on the full 128 coordinates the search compares about three quarters of the database to certify its answer, little better than a scan. Rotating the same 128-dimensional data onto its principal axes cuts that to under 30%, the effect Silpa-Anan and Hartley exploited; in our reading, the descriptors’ variance is concentrated in a few directions, and splitting along them prunes far more.

In the table, best-bin-first search in a single tree, stopped after 25 or 125 leaves, misses the true neighbor often, and four randomized trees with the same leaf budget do clearly better in every seed. Random-hyperplane LSH spans the whole range: long keys in few tables are cheap and miss often, short keys in eight tables find the neighbor for 98% of queries but compare almost a tenth of the database. Among the methods that compare full vectors, the inverted file is the most efficient here, reaching recall 0.996 or more at under 5% cost with four probes, a sign, in our reading, that 20 source images with recurring structures give clusters that k-means captures well; uniformly spread data would favor it less. Product quantization ranks the database from 8-byte codes, 64 times smaller than the 512-byte vectors, and re-ranking its top 10 with exact distances recovers the true neighbor for 98% to 100% of queries; its cost column leaves out the scan of all 10,000 codes, 8 table lookups each.

The two recall columns separate easy from ambiguous queries. Most methods do about as well on the queries that pass the ratio test as on all queries, but PQ without re-ranking misses mainly ambiguous ones: its recall@1 rises from about 0.85 to 0.98 on the distinctive queries, because a quantization error larger than the gap between the first and second neighbor swaps them. The distance-ratio column shows that misses are mostly near misses: averaged over all queries, every setting with recall@1 above 0.8 returns a neighbor within 10% of the true distance, and only the single tree with 25 leaves and the 16-bit LSH, at about 1.5, return distant ones. Lowe’s argument for pairing BBF with the ratio test, that the hardest queries have many neighbors at similar distances and need no exact answer, fits PQ’s misses here; the trees and LSH, however, missed distinctive queries about as often as ambiguous ones. The ranking of methods in this example depends on the data and on counting distances rather than measuring time; it is not a benchmark result.

Common Misconceptions

• “Approximate means inaccurate.” At practical operating points indexes return the exact nearest neighbor for most queries; Muja and Lowe (2014) call more than 95% correct neighbors at two or more orders of magnitude over linear search common. Recall is a parameter, bought with time or memory.
• “ANN recall is retrieval accuracy.” It measures agreement with exact search, not relevance; see Formal Definition.
• “High dimension makes nearest-neighbor search hopeless.” The negative results concern worst-case or uniformly distributed data. Descriptors and embeddings usually have a much lower intrinsic dimension than their length, as Muja and Lowe’s image patches show.
• “LSH has guarantees, so it is the best choice.” The guarantees are worst-case bounds. On real descriptors, data-adapted methods have usually been faster at equal recall (Jégou et al., 2011; Muja & Lowe, 2014).
• “An index is built once.” Learned centroids, codebooks and graphs fit one data distribution; a new embedding model means new vectors and a rebuilt index.

Where It Is Used

• Local feature matching. SIFT matching for object recognition used BBF (Lowe, 2004), and descriptor matching in feature matching, structure from motion and SLAM pipelines needs the same search.
• Visual vocabularies. Building and applying a bag of visual words vocabulary assigns millions of descriptors to their nearest centroids, which that article describes doing with k-d forests or a vocabulary tree.
• Image retrieval. Instance and category image retrieval rank a gallery of global embeddings or aggregated descriptors by distance to the query, with inverted files and product quantization when the gallery holds millions to billions of images.
• Cross-modal search. Image–text retrieval with models such as CLIP searches normalized embeddings by inner product, for example with FAISS.
• Verification and recognition at scale. Face recognition and person re-identification search a gallery of metric-learning embeddings; nearest-neighbor classifiers and the evaluation of self-supervised features with $k$-nearest-neighbor votes do the same over training sets.
• Deduplication and graph construction. Building a $k$-nearest-neighbor graph over a whole collection, for near-duplicate detection or clustering, is $n$ searches against one index; Johnson et al. (2017 arXiv version) built such a graph over one billion vectors.