You've embedded a hundred million passages with a 768-dimensional model. A query arrives, and the product team wants the top 10 in 20 milliseconds. How hard can that be? Let's count. Scoring one query against every vector is $10^8 \times 768 \approx 7.7\times 10^{10}$ multiply-adds. Worse, it means reading all the vectors, which in float32 is 307 GB. A server CPU socket streams memory at very roughly 100 to 200 GB per second, so the scan alone takes a couple of seconds, a hundred times over budget, before we've even sorted anything.

This section is about doing much less work and still finding (almost) the right answers. That's approximate nearest neighbor (ANN) search, and four ideas cover most of what's deployed. Inverted files skip most of the data. Product quantization shrinks every vector to a few bytes. HNSW graphs walk toward the answer in a few hundred steps. Locality-sensitive hashing buckets similar vectors together. We'll build each from scratch, derive what it costs, and finish with how to choose.

Exact search, and when it's enough

Set the problem up precisely. We have a database $\mathbf{X} \in \mathbb{R}^{n \times d}$ (row $\mathbf{x}_i$), a query $\mathbf{y} \in \mathbb{R}^d$, and we want the $k$ rows with the smallest squared distance $\lVert \mathbf{y} - \mathbf{x}_i\rVert^2$, or equivalently for unit vectors the largest inner product (§1.3 showed the rankings agree). Brute force computes all $n$ scores in $O(nd)$ time and selects the top $k$ in $O(n)$ with a partial sort (Exercise 1.3b).

The expansion $\lVert \mathbf{y} - \mathbf{x}\rVert^2 = \lVert\mathbf{y}\rVert^2 - 2\,\mathbf{y}^\top\mathbf{x} + \lVert\mathbf{x}\rVert^2$ turns a batch of queries into one matrix multiply, which is why exact search on a GPU is fast. With a batch of $b$ queries, each byte of the database read from memory is reused $b$ times, so the work goes from memory-bound to compute-bound, the roofline story of §11.2. For a single query it stays memory-bound, and the bytes are what you pay for.

So exact search is not a strawman. Below a million or so vectors it's often the right choice, and it's always the reference you measure ANN against. The standard measure is recall@$k$, the fraction of the true top-$k$ that the approximate method returns, $|\text{ANN}_k \cap \text{exact}_k|/k$. Notice what it compares against. This recall treats exact search as the truth and needs no relevance labels, so it measures the index alone. The retrieval recall of §7.5 compares against human relevance labels and measures the whole system, embedding model included. A perfect index over a poor embedding model scores 1.0 on the first and badly on the second.

One practical detail even exact search needs is blocking. Scoring 10 million vectors for a batch of 1,000 queries would need a $1000 \times 10^7$ score matrix, 40 GB. Instead you scan the database in blocks and keep a running top-$k$ per query, merging each block's best candidates into it.

import numpy as np

def blocked_topk(Qs, X, k, block):
    raise NotImplementedError

Inverted files, search only where the answer probably is

The first idea is the one a librarian would have. Don't scan every book, go to the right shelf. Cluster the database with k-means into $n_{\text{list}}$ cells, with centroids $\mathbf{c}_1, \dots, \mathbf{c}_{n_{\text{list}}}$. Store each vector in the list of its nearest centroid. That's an inverted file index (IVF), named after the inverted indexes of text search (§7.5), with cells playing the role of words.

At query time, compare the query with all centroids, pick the $n_{\text{probe}}$ nearest cells, and scan only the vectors in those lists. Why more than one cell? Because the true nearest neighbor of a query near a cell boundary often lives in the cell next door. $n_{\text{probe}}$ is the knob that trades speed for recall.

The cost is easy to write down. If the lists are balanced, each holds about $n/n_{\text{list}}$ vectors, and a query computes

multiply-adds. In plain words, more cells make the lists shorter but the centroid scan longer. Setting the derivative with respect to $n_{\text{list}}$ to zero, $d - n_{\text{probe}}\,n\,d/n_{\text{list}}^2 = 0$, gives

So IVF turns $O(n)$ into $O(\sqrt{n})$ per query, for a fixed $n_{\text{probe}}$. That fine print matters. As $n$ grows with $n_{\text{list}} \propto \sqrt{n}$, every cell holds more vectors and the space is cut into more cells, so in practice $n_{\text{probe}}$ has to creep up to keep the same recall, and the real scaling is somewhat worse than $\sqrt{n}$. For $n = 10^8$ and $n_{\text{probe}} = 16$, that's about 40,000 cells and roughly 80,000 distance computations (40,000 centroids plus 40,000 vectors in the probed lists) instead of 100 million. FAISS's guidelines suggest $n_{\text{list}}$ between $4\sqrt{n}$ and $16\sqrt{n}$, which is exactly our optimum $\sqrt{n_{\text{probe}}}\,\sqrt{n}$ for $n_{\text{probe}}$ between 16 and 256, and in practice the centroid scan itself is often accelerated with a small HNSW graph (more on HNSW below) once $n_{\text{list}}$ reaches the tens of thousands.

import numpy as np

def ivf_build(X, C):
    raise NotImplementedError

def ivf_search(y, X, C, lists, nprobe, k):
    raise NotImplementedError

Here's the speed and recall trade-off on 10,000 vectors in 32 dimensions with only mild cluster structure, which is closer to real embeddings than neat, well-separated blobs.

Probing one cell finds about half the true neighbors. Recall climbs steadily with $n_{\text{probe}}$, and the last few percent are expensive, which is the shape you'll see on every ANN benchmark. The widget below lets you see why in two dimensions. The neighbors you miss are the ones across a cell boundary.

Product quantization, a few bytes per vector

IVF reduces how many vectors we touch. It doesn't reduce how big they are, and at a billion vectors the index must also fit in memory. Product quantization (PQ, Jégou, Douze & Schmid, 2011) compresses each vector to a handful of bytes while still allowing fast approximate distances. It's the most important compression idea in vector search, so let's derive it carefully.

Start with ordinary vector quantization. Run k-means with $K$ centroids on the database and store each vector as the id of its nearest centroid, $\log_2 K$ bits. Distances to the query are approximated by distances to the centroid. To make the approximation good you need a huge $K$, and a codebook of $2^{64}$ centroids can be neither trained nor stored.

The product trick. Split each vector into $m$ contiguous subvectors of $d_s = d/m$ dimensions, $\mathbf{x} = [\mathbf{u}_1(\mathbf{x}), \dots, \mathbf{u}_m(\mathbf{x})]$ with $\mathbf{u}_j(\mathbf{x}) \in \mathbb{R}^{d_s}$. Learn a separate small codebook for each subspace, $\mathbf{C}^{(j)} \in \mathbb{R}^{K \times d_s}$ (row $\mathbf{C}^{(j)}_c$), by k-means on that subspace alone. A vector is encoded as $m$ ids,

and reconstructed by concatenation, $q(\mathbf{x}) = [\mathbf{C}^{(1)}_{c_1}, \dots, \mathbf{C}^{(m)}_{c_m}]$.

In plain words, PQ quantizes each slice of the vector separately with its own small k-means, and the code is the list of slice ids.

The implied codebook is the Cartesian product of the small ones, $K^m$ centroids, but we only store $m \cdot K \cdot d_s = K d$ numbers. With $K = 256$ (one byte per id) and $m = 8$, that's $2^{64}$ effective centroids in 8 bytes per vector.

Why is training each subspace independently the right thing to do? Because squared distance splits over the subvectors,

Each term involves only one codebook, so the total distortion is minimized by minimizing each term separately, which is exactly what k-means on each subspace does. The only approximation is the structural one, assuming the product of small codebooks is a good codebook. It works best when the subspaces carry similar amounts of variance, which is why optimized PQ (OPQ, Ge et al., 2013) first applies a learned rotation to balance them.

Asymmetric distances, computing with codes

Now the clever part. To rank the database for a query $\mathbf{y}$, we never decompress anything. The same split gives

For each query, precompute a lookup table $\boldsymbol{\Delta} \in \mathbb{R}^{m \times K}$ with $\Delta_{j,c} = \lVert \mathbf{u}_j(\mathbf{y}) - \mathbf{C}^{(j)}_c\rVert^2$. That costs $K d$ operations, once per query. Then the distance to any database vector is $m$ table lookups and $m - 1$ additions. In plain words, the query stays exact, the database is compressed, and distances become sums of precomputed numbers. This is called asymmetric distance computation (ADC).

The alternative, symmetric distance computation (SDC), quantizes the query too and looks up centroid-to-centroid distances. It costs the same, but it's less accurate, and the triangle inequality tells you exactly how much.

For inner-product search the same trick works, since $\mathbf{y}^\top q(\mathbf{x}) = \sum_j \mathbf{u}_j(\mathbf{y})^\top \mathbf{C}^{(j)}_{c_j}$ also splits into per-subspace lookups.

import numpy as np

def pq_train(X, m, K):
    raise NotImplementedError

def pq_encode(X, books):
    raise NotImplementedError

Let's put numbers on it with the 32-dimensional data from before, cut into $m = 8$ subspaces of 4 dimensions with $K = 64$ centroids each, and compare ADC with SDC.

The vectors shrank 16 times. On its own, PQ ranking gets about half of the exact top 10. But more than 90 percent of the true top 10 sit within PQ's top 100, so the standard recipe is to shortlist with PQ and rerank the shortlist with exact distances, from full vectors kept on SSD or from a finer code. ADC beats SDC on both counts, as the bound predicted.

In production, PQ is almost always combined with IVF, a design Jégou et al. called IVFADC. The coarse quantizer picks cells, and PQ encodes each vector's residual $\mathbf{x} - \mathbf{c}_{\text{cell}}$, which is smaller and more uniform than $\mathbf{x}$ itself, so the same bytes buy more accuracy. FAISS writes this as IVF65536,PQ64. A billion 768-dimensional vectors (3 TB in float32) then fit in about 64 GB of codes plus 8 GB of ids.

HNSW, navigating a graph of neighbors

The third idea comes from small-world networks. Connect each vector to some of its near neighbors, forming a proximity graph. To search, start anywhere, look at your current node's neighbors, move to whichever is closest to the query, and repeat until no neighbor is closer. That's greedy routing, and on a good graph it reaches the neighborhood of the answer in surprisingly few hops.

Plain proximity graphs have a problem. If every edge is short, getting across the space takes many hops. Hierarchical Navigable Small World graphs (HNSW, Malkov & Yashunin, 2018) fix it with layers, in the spirit of a skip list.

• Every vector lives in layer 0, a graph where each node links to up to $M_0 = 2M$ neighbors.
• Each vector is also placed in layers $1, 2, \dots, \ell$, where its top layer $\ell$ is random, $\ell = \lfloor -\ln(U)\, m_L \rfloor$ with $U$ uniform on $(0, 1]$. Upper layers hold exponentially fewer nodes with links of up to $M$ neighbors, so their edges are long.
• Search starts at a fixed entry point in the top layer, greedily descends to the closest node there, uses it as the entry to the layer below, and repeats. In layer 0 it switches from greedy to a beam search that keeps the $ef$ best candidates found so far, and returns the best $k$ of them.

The level distribution is where the logarithm comes from. For any level $\lambda \ge 0$, $\mathbb{P}(\ell \ge \lambda) = \mathbb{P}(-\ln U \ge \lambda/m_L) = e^{-\lambda/m_L}$, so choosing $m_L = 1/\ln M$, the paper's recommended setting, gives $\mathbb{P}(\ell \ge \lambda) = M^{-\lambda}$. Each layer has $M$ times fewer nodes than the one below, so there are about $\log_M n$ layers, and with a bounded number of hops per layer the search visits $O(\log n)$ nodes. The paper reports logarithmic complexity scaling empirically. It's not a worst-case guarantee.

Building the graph uses search. To insert a vector, draw its level, search down to that level, and in each layer from there to 0 run a beam search with width $ef_{\text{construction}}$ to find candidate neighbors, then connect it to $M$ of them. Which $M$? Not the closest. HNSW's neighbor-selection heuristic keeps a candidate only if it's closer to the new node than to every neighbor already selected. That prefers neighbors in different directions, which keeps long bridges between clusters, and the paper found that it helps most at high recall and on highly clustered data, where plain nearest neighbors would leave each cluster connected only to itself.

The layer-0 beam search is the inner loop of every graph index, so it's worth implementing once.

import numpy as np, heapq

def beam_search(X, graph, y, entry, ef):
    raise NotImplementedError

On a simple nearest-neighbor graph (no layers), here is how the beam width trades recall for work. The work is counted in distance computations against a database of 2,000 vectors.

Even without layers, a beam of 10 finds 97 percent of the true top 10, and a beam of 20 finds all of them while computing distances to about a fifth of the database. Wider beams buy nothing more here and keep costing distance computations. The layers of HNSW mostly speed up the approach from a far-away entry point, which matters more as $n$ grows.

HNSW is the default in most vector databases because its recall-versus-latency curve is among the best on public benchmarks (ann-benchmarks, Aumüller et al., 2020) and because it supports inserting vectors one at a time. Its costs are memory, since each vector carries about $2M$ links in layer 0 (with $M = 32$ and 4-byte ids, 256 bytes per vector on top of the vector itself), slow builds, and awkward deletions (usually handled by marking nodes deleted and rebuilding periodically). The graph must also live in RAM. DiskANN (Subramanya et al., 2019) builds a single-layer graph designed for SSDs, keeping compressed vectors in memory and full vectors on disk, to serve a billion vectors from one machine.

Locality-sensitive hashing

The oldest idea gives the cleanest theory. You built a random-hyperplane index in §1.3; here we derive its collision probability and the amplification that makes it useful. Design a random hash function so that similar vectors are likely to collide and dissimilar ones aren't, then look only at vectors in the query's buckets. For cosine similarity the classic choice is the random hyperplane (Charikar, 2002). Draw $\mathbf{r} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ and set $h_{\mathbf{r}}(\mathbf{x}) = \operatorname{sign}(\mathbf{r}^\top\mathbf{x})$. Then

where $\theta$ is the angle between $\mathbf{x}$ and $\mathbf{y}$.

One bit is a weak filter. LSH amplifies it in two ways. Concatenate $r$ bits into one bucket key, so two vectors share a bucket with probability $p^r$ where $p = 1 - \theta/\pi$ (an AND, which sharpens). Then build $L$ independent tables and take the union of the query's buckets, so a pair collides in at least one table with probability

In plain words, $r$ makes buckets strict and $L$ gives each near neighbor many chances to be found, and together they turn the gentle line $1 - \theta/\pi$ into a steep S-curve that's near 1 for close pairs and near 0 for far ones. You'll meet exactly this banding construction, with MinHash in place of hyperplanes, for near-duplicate detection in §8.3.

With 16 bits and 40 tables, pairs 20 degrees apart are found 99.9 percent of the time and pairs 60 degrees apart only 6 percent. LSH has provable sublinear query time for the approximate problem, returning a point within a factor $c \gt 1$ of the nearest distance with high probability (Indyk & Motwani, 1998), and it is simple to shard and update. For dense learned embeddings, though, graph and quantization methods usually reach the same recall with far less memory and time, because LSH ignores the data's actual distribution and needs many tables. You'll still meet it in streaming systems, deduplication, and behind the sign-bit codes of §7.3, which are hyperplane hashes whose hyperplanes are the coordinate hyperplanes (the sign of coordinate $j$ is the sign of $\mathbf{e}_j^\top\mathbf{x}$). Because the axes of a trained model are not random directions, the $1 - \theta/\pi$ formula holds only approximately for them, and a random rotation before binarizing restores it bit by bit.

Choosing an index

There's no best index, only a best index for a budget. Here's a map, with FAISS's factory names in parentheses.

A few rules of thumb that hold up in practice.

• Measure recall against exact search on a sample of real queries before tuning anything. A latency number without a recall number is meaningless.
• The three-way trade-off is recall, latency, and memory. You can usually buy any two. Compression buys memory with recall, which reranking buys back with latency.
• Filters are hard. "Top 10 documents from this customer, newer than last week" breaks naive ANN. Filtering after search can return almost nothing, and filtering before search defeats the index. Vector databases handle this with filter-aware graph traversal or partitioning by the filter key.
• Normalize consistently. If the model is trained for cosine, normalize at index and query time and use inner product. If lengths carry meaning (raw dot products, as in many two-tower recommenders), you need maximum inner-product search, so use an index built for inner product or the norm-augmentation reduction of §1.3, which turns it into nearest-neighbor search. ScaNN (Guo et al., 2020) goes further and trains its quantizer to preserve inner products rather than distances.

In PyTorch

On a GPU, exact search is a matrix multiply plus topk, done in chunks to bound memory, and it's often fast enough to skip ANN entirely below tens of millions of vectors.

With FAISS, the same ideas are a few lines each. This needs the faiss package, so it's shown and not run.

Read more

• Jégou, Douze & Schmid (2011). Product Quantization for Nearest Neighbor Search. PQ, ADC versus SDC, and IVFADC. The derivations in this section follow it.
• Malkov & Yashunin (2018). Efficient and robust approximate nearest neighbor search using Hierarchical Navigable Small World graphs. HNSW, the level distribution, and the neighbor heuristic.
• Douze et al. (2024), The Faiss library, and Johnson, Douze & Jégou (2019), Billion-scale similarity search with GPUs. A modern overview of the design space with the trade-offs measured, and how exact and IVF-PQ search map onto GPUs.
• Charikar (2002). Similarity Estimation Techniques from Rounding Algorithms. Random hyperplane hashing.
• Subramanya et al. (2019). DiskANN: Fast Accurate Billion-point Nearest Neighbor Search on a Single Node. Billion-scale graph search from SSD.
• Guo et al. (2020). Accelerating Large-Scale Inference with Anisotropic Vector Quantization (ScaNN). Quantization that preserves inner products, behind one of the fastest MIPS libraries.
• Aumüller, Bernhardsson & Faithfull (2020). ANN-Benchmarks: A Benchmarking Tool for Approximate Nearest Neighbor Algorithms. How to compare ANN methods fairly, and how they compare.
• The FAISS wiki. Guidelines to choose an index. Practical rules of thumb, including the cell-count guideline.

Watch next

▶

▶

▶

▶