9.2 Vector Databases & Approximate Nearest Neighbor Search¶

In the previous chapter, Embeddings & Representation Learning, we learned how to turn a chunk of text into a dense vector \(\mathbf{x} \in \mathbb{R}^d\) such that semantically similar passages land close together in vector space. That is half of a retrieval system. The other half is the problem this chapter solves: given a query vector \(\mathbf{q}\) and a corpus of \(N\) vectors — where \(N\) might be ten million, or ten billion — find the \(k\) vectors closest to \(\mathbf{q}\) fast enough to put in the critical path of an LLM request. A user will not wait two seconds for retrieval before the model even starts generating. You have a few milliseconds.

This is the nearest neighbor search problem, and it is deceptively hard. The naive solution — compare the query to every vector — is correct but linear in \(N\). At a billion vectors, that is a billion dot products per query. The entire field of Approximate Nearest Neighbor (ANN) search exists to trade a sliver of correctness for orders of magnitude in speed. The systems that productionize these algorithms — FAISS, Milvus, Qdrant, Weaviate, pgvector, ScaNN — are vector databases. By the end of this chapter you will understand the three algorithm families that power essentially all of them (graphs, inverted files, and quantization), you will have built a brute-force searcher, a product quantizer, and a working HNSW index from scratch, and you will be able to reason quantitatively about the recall-latency-memory triangle that governs every deployment decision.

This is foundational plumbing for Retrieval-Augmented Generation Architectures, Chunking, Reranking & Hybrid Search, and Memory Systems for Agents. It is also one of the most reliably asked ML-system-design topics in interviews, which we address directly.

The Problem: Exact Nearest Neighbor Search¶

Formalizing the query¶

We have a database of vectors \(\mathcal{D} = \{\mathbf{x}_1, \dots, \mathbf{x}_N\}\), each \(\mathbf{x}_i \in \mathbb{R}^d\). Given a query \(\mathbf{q} \in \mathbb{R}^d\) and a distance function \(\operatorname{dist}(\cdot, \cdot)\), the \(k\)-nearest-neighbor (\(k\)-NN) problem is to return the set

i.e. the \(k\) database vectors with the smallest distance to \(\mathbf{q}\). The two distance functions that matter in practice are squared Euclidean (L2) distance and, far more commonly for normalized text embeddings, cosine similarity (equivalently, inner product on unit vectors):

A crucial identity ties these together. For unit-normalized vectors (\(\|\mathbf{q}\| = \|\mathbf{x}\| = 1\)):

So minimizing L2 distance is equivalent to maximizing inner product when everything is normalized. This is why most retrieval pipelines L2-normalize embeddings at index time: it lets you use a single index that is simultaneously a cosine index and a Euclidean index. The flip side — maximum inner product search (MIPS) on un-normalized vectors — is genuinely harder, because inner product is not a metric (it violates the triangle inequality), which breaks many tree- and graph-based pruning arguments.

Brute force: the baseline you must always have¶

The exact answer is computed by scanning every vector. Done in a vectorized way, it is just a matrix-vector (or matrix-matrix, for a batch of queries) product:

import numpy as np

def brute_force_knn(data: np.ndarray, q: np.ndarray, k: int = 10):
    """
    Exact k-NN under squared-L2 distance.
    data: (N, d) float32 corpus.   q: (d,) float32 query.
    Returns (indices, distances) of the k closest vectors, sorted ascending.
    """
    # ||x - q||^2 = ||x||^2 - 2 x.q + ||q||^2.  The ||q||^2 term is constant
    # across all x, so it does not affect the ranking and we can drop it.
    # We keep the full form here for an honest distance value.
    diffs = data - q                  # (N, d) broadcast subtraction
    dists = np.einsum("nd,nd->n", diffs, diffs)  # (N,) squared L2, fused

    # argpartition is O(N): it finds the k smallest without fully sorting.
    idx = np.argpartition(dists, k)[:k]
    # Then sort just those k for a clean ranked list.
    idx = idx[np.argsort(dists[idx])]
    return idx, dists[idx]

# For a *batch* of queries, prefer the GEMM trick:
def brute_force_batch(data: np.ndarray, queries: np.ndarray, k: int = 10):
    """data: (N,d), queries: (Q,d) -> top-k indices per query, (Q,k)."""
    # ||x-q||^2 = ||x||^2 - 2 q.x + ||q||^2. Compute the cross term with one GEMM.
    x_norm = np.einsum("nd,nd->n", data, data)        # (N,)
    cross  = queries @ data.T                          # (Q, N)  the expensive part
    dists  = x_norm[None, :] - 2.0 * cross             # (Q, N), q_norm dropped
    return np.argpartition(dists, k, axis=1)[:, :k]

Two lessons hide in this code. First, the dominant cost is the queries @ data.T GEMM (general matrix multiply), which is \(O(Q \cdot N \cdot d)\) FLOPs but maps perfectly onto BLAS / a GPU — see GPU Architecture & The Memory Hierarchy. On a modern GPU, brute force over a few million vectors is genuinely fast and is often the right answer; FAISS ships a IndexFlatIP/IndexFlatL2 precisely for this. Second, brute force gives recall = 1.0 by definition, so it is your ground-truth oracle for measuring how approximate the approximate methods are. Never ship an ANN index without measuring its recall against brute force on a held-out query set.

The Curse of Dimensionality¶

Why can’t we just use a \(k\)-d tree, the classic exact NN data structure from any algorithms course? Because in high dimensions, space itself stops behaving the way our 2-D and 3-D intuitions expect. This phenomenon is the curse of dimensionality, and understanding it is the key to understanding why ANN methods look the way they do.

Distances concentrate¶

Consider \(N\) points drawn uniformly at random in the unit hypercube \([0,1]^d\). As \(d\) grows, the distance from a query to its nearest neighbor and to its farthest neighbor become nearly identical. Concretely, a classic result (Beyer et al., When Is “Nearest Neighbor” Meaningful?) shows that under broad conditions,

The contrast between near and far vanishes. We can feel this with a one-line experiment:

import numpy as np
rng = np.random.default_rng(0)

for d in [2, 10, 100, 1000]:
    X = rng.random((10_000, d))          # 10k uniform points in [0,1]^d
    q = rng.random(d)
    dists = np.linalg.norm(X - q, axis=1)
    dmin, dmax = dists.min(), dists.max()
    # "contrast": how much farther the farthest point is vs the nearest
    print(f"d={d:5d}  dmin={dmin:7.3f}  dmax={dmax:7.3f}  "
          f"ratio dmax/dmin={dmax/dmin:5.2f}")
# d=    2  ... ratio ~ 50x      (nearest is dramatically nearer)
# d= 1000  ... ratio ~ 1.1x     (everything is roughly equidistant)

At \(d=2\) the farthest point is dozens of times farther than the nearest, so “nearest neighbor” is meaningful and a tree can prune aggressively. At \(d=1000\) the ratio collapses toward 1: the nearest and farthest neighbors are almost the same distance away.

Why this kills exact structures¶

A \(k\)-d tree prunes a branch when the splitting hyperplane is farther from the query than the best candidate found so far. In high dimensions, because distances concentrate and because a query’s “ball” of radius (best-so-far) intersects almost every cell, almost no branches can be pruned. The tree degrades to a linear scan — with worse constants than brute force, because of the pointer chasing. Empirically, \(k\)-d trees lose to brute force somewhere around \(d \approx 20\)–\(30\).

There is a second, geometric way to see the curse. The volume of a \(d\)-dimensional unit ball relative to the enclosing cube goes to zero super-exponentially:

Almost all the volume of a high-dimensional cube is in its corners, far from the center. Random data is spread thin across an enormous space, and “locality” — the property tree structures exploit — barely exists.

Inverted Files (IVF): Search by Coarse Partition¶

The first practical idea is brutally simple: don’t search the whole database, search a neighborhood of it. Partition the vector space into cells, figure out which cell the query falls in, and scan only that cell (plus a few neighbors). This is the Inverted File index, IVF, named by analogy to the inverted index of classical text search.

Building an IVF index¶

1. Train by running \(k\)-means on a sample of the corpus to find \(n_{\text{list}}\) centroids \(\{\mathbf{c}_1, \dots, \mathbf{c}_{n_{\text{list}}}\}\). These centroids are the coarse quantizer. A common heuristic is \(n_{\text{list}} \approx \sqrt{N}\).
2. Assign every database vector to its nearest centroid. The index is then a set of \(n_{\text{list}}\) lists (the “inverted lists” or posting lists), where list \(j\) holds the ids of all vectors assigned to centroid \(\mathbf{c}_j\).
3. Search by finding the \(n_{\text{probe}}\) centroids nearest to \(\mathbf{q}\), then doing brute force only over the vectors in those \(n_{\text{probe}}\) lists.

import numpy as np

class IVFFlat:
    """Inverted-file index with exact (flat) distances inside each cell."""
    def __init__(self, n_list: int = 256, n_iter: int = 25, seed: int = 0):
        self.n_list = n_list
        self.n_iter = n_iter
        self.rng = np.random.default_rng(seed)

    def _kmeans(self, X):
        n = X.shape[0]
        # init centroids from random data points (k-means++ is better in prod)
        C = X[self.rng.choice(n, self.n_list, replace=False)].copy()
        for _ in range(self.n_iter):
            # assign: nearest centroid for each point via the (a-b)^2 GEMM trick
            d2 = (np.einsum("nd,nd->n", X, X)[:, None]
                  - 2 * X @ C.T
                  + np.einsum("kd,kd->k", C, C)[None, :])   # (n, n_list)
            assign = d2.argmin(1)
            for j in range(self.n_list):                    # update step
                pts = X[assign == j]
                if len(pts):
                    C[j] = pts.mean(0)
        return C

    def train(self, X):
        self.centroids = self._kmeans(X)
        # assign every vector to its nearest centroid -> build posting lists
        d2 = (np.einsum("nd,nd->n", X, X)[:, None]
              - 2 * X @ self.centroids.T
              + np.einsum("kd,kd->k", self.centroids, self.centroids)[None, :])
        assign = d2.argmin(1)
        self.lists = [np.where(assign == j)[0] for j in range(self.n_list)]
        self.data = X
        return self

    def search(self, q, k=10, n_probe=8):
        # 1) find the n_probe nearest centroids to the query
        cd = np.einsum("kd,kd->k", self.centroids - q, self.centroids - q)
        probes = np.argpartition(cd, n_probe)[:n_probe]
        # 2) gather candidate ids from those lists and do exact distances
        cand = np.concatenate([self.lists[j] for j in probes])
        if len(cand) == 0:
            return np.array([], dtype=int)
        dists = np.einsum("nd,nd->n", self.data[cand] - q, self.data[cand] - q)
        order = np.argpartition(dists, min(k, len(cand) - 1))[:k]
        order = order[np.argsort(dists[order])]
        return cand[order]

The IVF tradeoff knob: n_probe¶

n_probe is the single most important search-time parameter. With n_probe = 1 you scan only the query’s own cell — blazing fast, but you miss any true neighbor that happens to sit just across a cell boundary (the boundary problem). With n_probe = n_list you scan everything and recover exact brute force. In between, recall rises monotonically with n_probe while latency rises roughly linearly. The expected number of vectors scanned is approximately

so with \(n_{\text{list}} = \sqrt{N}\) and \(n_{\text{probe}} = 8\), you scan on the order of \(8\sqrt{N}\) vectors instead of \(N\). For \(N = 10^9\) that is \(8 \cdot 31{,}623 \approx 2.5 \times 10^5\) vectors per query — a 4000× reduction.

IVF alone solves the speed problem but not the memory problem: every vector is still stored in full float32. For a billion 768-d vectors that is 3 TB of RAM. To break the memory wall we need quantization.

Product Quantization (PQ): Compressing the Vectors¶

Product Quantization, introduced by Jégou, Douze & Schmid (Product Quantization for Nearest Neighbor Search, 2011), is the idea that lets a vector database hold billions of vectors in gigabytes of RAM. It is one of the most elegant ideas in the field.

The core trick: split, then quantize each piece¶

A scalar quantizer maps a single number to one of \(K\) codebook values. To quantize a whole \(d\)-dimensional vector to \(K\) codewords, you would need \(K\) exponentially large to cover the space — infeasible. PQ sidesteps this: split the vector into \(m\) contiguous subvectors and quantize each subvector independently with its own small codebook of \(k^*\) entries (typically \(k^* = 256\), so each code fits in one byte).

Split \(\mathbf{x} \in \mathbb{R}^d\) into \(m\) subvectors \(\mathbf{x} = (\mathbf{u}^1, \dots, \mathbf{u}^m)\), each of dimension \(d^* = d/m\). Learn a separate codebook \(\mathcal{C}^j = \{\mathbf{c}^j_1, \dots, \mathbf{c}^j_{k^*}\}\) per subspace via \(k\)-means. The PQ code of \(\mathbf{x}\) is the tuple of nearest-centroid indices:

The genius is in the arithmetic. The reconstruction codebook implicitly represents \((k^*)^m\) distinct vectors — for \(k^* = 256\), \(m = 8\) that is \(256^8 = 2^{64}\) possible reconstructions — but you only store \(m \cdot k^* \cdot d^*\) floats of codebooks and \(m\) bytes per database vector. A 768-d float32 vector (3072 bytes) compresses to \(m = 8\) bytes: a 384× memory reduction.

Asymmetric Distance Computation (ADC)¶

How do you compute the distance from a full-precision query to a compressed database vector without decompressing? You exploit the separability of squared L2 distance across the subspaces. Approximate

where \(\mathbf{q}^j\) is the \(j\)-th sub-block of the query. The trick: at query time, precompute a small lookup table (LUT) of size \(m \times k^*\) holding the squared distance from each query sub-block to every centroid in that subspace. Then the distance to any database vector is just \(m\) table lookups and \(m{-}1\) additions — no multiplications, no decompression. This is ADC (the query stays “asymmetric”: full precision; only the database is quantized, which keeps accuracy higher than quantizing both sides).

import numpy as np

class ProductQuantizer:
    def __init__(self, m=8, k_star=256, n_iter=20, seed=0):
        self.m, self.k_star, self.n_iter = m, k_star, n_iter
        self.rng = np.random.default_rng(seed)

    def train(self, X):
        n, d = X.shape
        assert d % self.m == 0, "d must be divisible by m"
        self.dsub = d // self.m
        # one codebook per subspace: (m, k_star, dsub)
        self.codebooks = np.zeros((self.m, self.k_star, self.dsub), np.float32)
        for j in range(self.m):
            sub = X[:, j*self.dsub:(j+1)*self.dsub]          # (n, dsub)
            C = sub[self.rng.choice(n, self.k_star, replace=False)].copy()
            for _ in range(self.n_iter):                     # k-means per subspace
                d2 = (np.einsum("nd,nd->n", sub, sub)[:, None]
                      - 2 * sub @ C.T
                      + np.einsum("kd,kd->k", C, C)[None, :])
                a = d2.argmin(1)
                for c in range(self.k_star):
                    pts = sub[a == c]
                    if len(pts): C[c] = pts.mean(0)
            self.codebooks[j] = C
        return self

    def encode(self, X):
        """Map each vector to m uint8 codes."""
        n = X.shape[0]
        codes = np.zeros((n, self.m), np.uint8)
        for j in range(self.m):
            sub = X[:, j*self.dsub:(j+1)*self.dsub]
            C = self.codebooks[j]
            d2 = (np.einsum("nd,nd->n", sub, sub)[:, None]
                  - 2 * sub @ C.T
                  + np.einsum("kd,kd->k", C, C)[None, :])
            codes[:, j] = d2.argmin(1)                       # nearest centroid id
        return codes

    def search(self, q, codes, k=10):
        """ADC: build the per-subspace LUT, then sum lookups for every code."""
        # 1) precompute LUT: distance from each query sub-block to all centroids
        lut = np.empty((self.m, self.k_star), np.float32)
        for j in range(self.m):
            qsub = q[j*self.dsub:(j+1)*self.dsub]            # (dsub,)
            diff = self.codebooks[j] - qsub                  # (k_star, dsub)
            lut[j] = np.einsum("kd,kd->k", diff, diff)       # (k_star,)
        # 2) approximate distance to every db vector = sum of m table lookups
        approx = np.zeros(codes.shape[0], np.float32)
        for j in range(self.m):
            approx += lut[j, codes[:, j]]                    # vectorized gather
        idx = np.argpartition(approx, k)[:k]
        return idx[np.argsort(approx[idx])]

OPQ: rotate before you split¶

PQ carves \(\mathbf{x}\) into contiguous blocks, which quietly assumes that adjacent coordinates belong together and that variance is spread evenly across the \(m\) blocks. Real embeddings satisfy neither: some coordinate groups carry far more variance than others, so a uniform byte budget wastes bits on the flat blocks and starves the informative ones. Optimized Product Quantization (Ge et al., 2013) fixes this by learning an orthogonal rotation \(R\) and quantizing \(R\mathbf{x}\) instead of \(\mathbf{x}\), alternating between updating \(R\) and re-fitting the codebooks. Because \(R\) is orthogonal it preserves L2 distances and inner products exactly, so it costs one \(d \times d\) matrix-vector product on the query and nothing at ranking time. In FAISS it is the OPQ64_768 prefix of a factory string, and for a few extra training minutes it is almost always worth it.

Scalar and binary quantization: the cheaper cousins¶

PQ is not the only compressor, and it is often not the first one you should reach for.

Scalar quantization (SQ) quantizes each coordinate independently. At train time you record a per-dimension range (min/max, or better, a clipped percentile range to survive outliers); at index time each float32 becomes one uint8, giving a flat 4× reduction with essentially no training cost and a decode that is one fused multiply-add. Recall loss is typically under a point for normalized text embeddings. This is FAISS’s SQ8, Qdrant’s default scalar quantization, and — one step milder — pgvector’s halfvec (float16, 2×). Start here. If SQ8 plus a re-rank clears your eval bar, you never need PQ’s \(k\)-means machinery.

Binary quantization (BQ) takes the idea to its limit: store \(\operatorname{sign}(x_j)\) as a single bit, so a 768-d vector becomes 96 bytes — 32× — and distance becomes Hamming distance, computed with xor + popcount at many gigabytes per second. On its own it is crude, but because it is so cheap it makes an excellent first stage: shortlist a few hundred candidates by Hamming, then re-score them with full-precision (or SQ8) vectors. RaBitQ (Gao & Long, SIGMOD 2024) is the principled version of this: randomly rotate the residual, normalize it onto the unit sphere, and binarize — yielding an unbiased inner-product estimator with a provable error bound, and reported to match or beat PQ at the same memory budget. Its extension generalizes the construction to any bit budget, and it is being adopted in production engines as a PQ replacement.

# Real-library usage: the two cheap quantizers, both with a re-rank stage.
# pip install faiss-cpu
import faiss, numpy as np

d, N = 768, 100_000
rng = np.random.default_rng(0)
xb = rng.standard_normal((N, d)).astype("float32"); faiss.normalize_L2(xb)
xq = rng.standard_normal((4, d)).astype("float32"); faiss.normalize_L2(xq)

# --- 4x: int8 scalar quantization inside IVF cells (train is fast) ---
sq = faiss.index_factory(d, "IVF1024,SQ8", faiss.METRIC_INNER_PRODUCT)
sq.train(xb); sq.add(xb)
faiss.ParameterSpace().set_index_parameter(sq, "nprobe", 16)
D, I = sq.search(xq, 10)

# --- 32x: 1 bit per dimension, Hamming shortlist, exact re-rank ---
bb = np.packbits(xb > 0, axis=1)            # (N, d/8) uint8 -- sign bits packed
bq = np.packbits(xq > 0, axis=1)
bidx = faiss.IndexBinaryFlat(d)             # d = number of BITS here
bidx.add(bb)
_, cand = bidx.search(bq, 200)              # (4, 200) Hamming shortlist
# Stage 2: re-score only the shortlist against full-precision vectors.
top = []
for qi in range(len(xq)):
    c = cand[qi]
    sims = xb[c] @ xq[qi]                   # exact inner products, 200 of them
    top.append(c[np.argsort(-sims)[:10]])
print(I.shape, np.array(top).shape)         # (4, 10) (4, 10)

IVF + PQ: the workhorse of billion-scale search¶

The legendary FAISS index IVF{n_list},PQ{m} combines all of the above. IVF restricts the search to a handful of cells (solving speed); PQ compresses each vector to \(m\) bytes (solving memory). A common refinement, IVFADC with residuals, quantizes not the raw vector but its residual after subtracting its coarse centroid, which dramatically improves accuracy because residuals are small and well-clustered. A final, optional re-ranking step (IndexRefineFlat) fetches the top candidates’ full vectors from disk and re-scores them exactly to recover the recall that quantization lost.

HNSW: Navigable Small-World Graphs¶

Quantization compresses; IVF partitions. The third family — graph-based ANN — takes a different route and is, for most in-memory workloads under ~100M vectors, the recall-latency champion. The dominant algorithm is Hierarchical Navigable Small World (HNSW) graphs, from Malkov & Yashunin (Efficient and robust approximate nearest neighbor search using HNSW graphs, 2016/2018).

Greedy search on a proximity graph¶

Build a graph where each vector is a node connected to some of its near neighbors. To search, start at any node and greedily walk downhill: repeatedly move to whichever neighbor of the current node is closer to the query, until no neighbor improves. On a well-constructed navigable graph, this greedy walk converges to (approximately) the true nearest neighbor in a number of hops that grows only logarithmically with \(N\).

The naive version gets stuck in local minima and has long-range connectivity problems. HNSW fixes both with two ideas:

1. Small-world structure. Inspired by Kleinberg’s navigable small-world networks (the “six degrees of separation” insight), the graph mixes short links (to immediate neighbors, for accuracy) with long links (for fast traversal across the space). This gives the logarithmic hop count.
2. A hierarchy of layers. Like a skip list, HNSW stacks multiple graph layers. The top layer is sparse with long-range links; each lower layer is denser. A node’s maximum layer is drawn from an exponentially decaying distribution, \(\ell = \lfloor -\ln(\text{Uniform}(0,1)) \cdot m_L \rfloor\) with \(m_L = 1/\ln(M)\). Search starts at the single top-layer entry point, greedily descends to the nearest node in each layer, then drops down — zooming in like a coarse-to-fine map.

  Layer 2:   E ------------------- A                (sparse, long hops)
                 \               /
  Layer 1:   E --- B ----- C --- A --- D            (medium)
              \    |       |     |    /
  Layer 0:   E-B-F-G-B-C-H-I-C-A-J-D-K...           (dense; every node lives here)
                          ^
              query enters at top entry E, greedily descends,
              refines in layer 0 to collect the final top-k.

From-scratch HNSW¶

Here is a complete, runnable HNSW. It is faithful to the paper’s two routines — SEARCH-LAYER (a best-first beam search with a candidate min-heap and a result max-heap) and the layered INSERT — and it achieves real recall. The two build parameters are \(M\) (neighbors per node) and ef_construction (beam width during insertion); the one search parameter is ef (beam width during query, which trades recall for latency).

import numpy as np, heapq, math

def l2(a, b):
    d = a - b
    return float(d @ d)               # squared L2

class HNSW:
    def __init__(self, data, M=16, ef_construction=100, seed=2):
        self.data = data
        self.M = M                        # max neighbors per node (layer > 0)
        self.M0 = 2 * M                   # layer 0 gets more links (paper heuristic)
        self.efc = ef_construction
        self.mL = 1.0 / math.log(M)       # level-generation normalizer
        self.rng = np.random.default_rng(seed)
        self.graph = []                   # graph[node][layer] -> list of neighbor ids
        self.levels = []                  # top layer index of each node
        self.entry = None                 # current global entry point
        self.top = -1                     # current top layer

    def _rand_level(self):
        # exponentially-decaying layer assignment (skip-list style)
        return int(-math.log(self.rng.random()) * self.mL)

    def _search_layer(self, q, entry_points, ef, layer):
        """Best-first search within one layer; returns ef nearest as (dist,id)."""
        visited = set(entry_points)
        cand = []   # min-heap of (dist, id): the frontier to expand
        res  = []   # max-heap of (-dist, id): the ef best found so far
        for e in entry_points:
            de = l2(q, self.data[e])
            heapq.heappush(cand, (de, e))
            heapq.heappush(res, (-de, e))
        while cand:
            dc, c = heapq.heappop(cand)
            # if the closest frontier node is worse than our current worst
            # kept result, we can stop expanding (the beam has converged).
            if -res[0][0] < dc:
                break
            for nb in self.graph[c][layer]:
                if nb in visited:
                    continue
                visited.add(nb)
                dn = l2(q, self.data[nb])
                if dn < -res[0][0] or len(res) < ef:
                    heapq.heappush(cand, (dn, nb))
                    heapq.heappush(res, (-dn, nb))
                    if len(res) > ef:
                        heapq.heappop(res)     # keep only the ef best
        return [(-nd, n) for nd, n in res]

    def _select_neighbors(self, node, candidates, M):
        """Simple heuristic: keep the M closest candidates to `node`."""
        scored = sorted((l2(self.data[node], self.data[c]), c) for c in candidates)
        return [c for _, c in scored[:M]]

    def insert(self, i):
        level = self._rand_level()
        self.graph.append([[] for _ in range(level + 1)])
        self.levels.append(level)

        if self.entry is None:                  # first node ever
            self.entry, self.top = i, level
            return

        ep = [self.entry]
        # Phase 1: greedily descend the layers ABOVE the new node's top layer,
        # using a beam of width 1 (just find the single best entry per layer).
        for lc in range(self.top, level, -1):
            ep = [min(self._search_layer(self.data[i], ep, 1, lc))[1]]

        # Phase 2: for each layer the node lives in, find efc candidates and link.
        for lc in range(min(level, self.top), -1, -1):
            W = self._search_layer(self.data[i], ep, self.efc, lc)
            Mmax = self.M0 if lc == 0 else self.M
            neighbors = self._select_neighbors(i, [n for _, n in W], Mmax)
            self.graph[i][lc] = neighbors
            for n in neighbors:                 # add the reverse edges
                self.graph[n][lc].append(i)
                if len(self.graph[n][lc]) > Mmax:   # prune n if it got too full
                    self.graph[n][lc] = self._select_neighbors(
                        n, self.graph[n][lc], Mmax)
            ep = [n for _, n in W]

        if level > self.top:                    # this node becomes the new entry
            self.entry, self.top = i, level

    def search(self, q, k=10, ef=50):
        ep = [self.entry]
        for lc in range(self.top, 0, -1):       # descend with beam width 1
            ep = [min(self._search_layer(q, ep, 1, lc))[1]]
        W = sorted(self._search_layer(q, ep, ef, 0))   # widen the beam at layer 0
        return [n for _, n in W[:k]]

When we ran this on 8,000 random 32-d vectors and measured recall@10 against brute force, sweeping the search beam width ef produced exactly the behavior the theory predicts — a smooth recall-latency dial:

ef= 10   recall@10 = 0.500     (fast, low recall)
ef= 25   recall@10 = 0.673
ef= 50   recall@10 = 0.823
ef=100   recall@10 = 0.895
ef=200   recall@10 = 0.955     (slow, high recall)

That table is the recall-latency tradeoff, made concrete. You turn one knob, ef, and slide along the curve.

def _select_neighbors_heuristic(self, node, candidates, M):
    """HNSW Algorithm 4: keep *diverse* neighbors, not merely the closest."""
    scored = sorted((l2(self.data[node], self.data[c]), c)
                    for c in set(candidates) if c != node)
    kept = []
    for d_nc, c in scored:              # consider the nearest candidate first
        # accept c only if it is closer to `node` than to every kept neighbor;
        # otherwise some kept neighbor already "covers" that direction.
        if all(d_nc < l2(self.data[c], self.data[r]) for r in kept):
            kept.append(c)
        if len(kept) >= M:
            break
    return kept

Why HNSW dominates — and its costs¶

HNSW gives the best recall-per-latency of any method for in-memory data, with no training step (unlike IVF/PQ’s \(k\)-means). Its costs are real, though:

• Memory. It stores full vectors plus the graph. Each node holds up to \(M\) neighbor ids per layer; layer 0 alone is \(\sim 2M\) ids \(\times\) 4–8 bytes. For \(M=16\), that is ~64–128 bytes of graph overhead per vector on top of the \(4d\) bytes of the vector itself. HNSW is memory-hungry, which is why billion-scale systems combine it with PQ (IndexHNSWPQ).
• Build time. Insertion is \(O(\log N)\) per node but with a large constant; building a 100M-node graph takes hours and is the dominant offline cost.
• Updates. Deletes are awkward (you typically tombstone and periodically rebuild). High churn favors IVF, which re-clusters more gracefully.

ScaNN and the Anisotropic Insight¶

Google’s ScaNN (Scalable Nearest Neighbors; Guo et al., Accelerating Large-Scale Inference with Anisotropic Vector Quantization, ICML 2020) is worth a section because it changed how the field thinks about what quantization should preserve. Classic PQ minimizes reconstruction error \(\|\mathbf{x} - q(\mathbf{x})\|^2\) uniformly. ScaNN’s insight: for maximum inner product search, not all reconstruction errors are equally harmful. An error parallel to the database vector changes its inner product with queries; an error orthogonal to it mostly does not. So you should penalize parallel (anisotropic) error more heavily.

ScaNN replaces the isotropic loss with an anisotropic one that weights the parallel residual component more than the orthogonal one:

By learning codebooks under this loss, ScaNN preserves the inner products that determine ranking, achieving higher recall at the same compression. Combined with a partitioning tree and a heavily SIMD-optimized in-register distance computation, ScaNN is one of the fastest libraries on standard benchmarks. The broader lesson generalizes well beyond ScaNN: optimize the quantizer for the task metric (ranking by inner product), not for a generic surrogate (reconstruction MSE).

The Recall–Latency–Memory Triangle¶

Every ANN deployment is a point inside a triangle whose three corners are recall (how often you find the true neighbors), latency (how fast), and memory (how much RAM/disk). You can have any two cheaply; the third is the price.

Measuring recall¶

Recall@\(k\) is the fraction of the true top-\(k\) that your approximate search returns:

Always compute exact_k with brute force on a held-out query set. A subtle but important point: recall is not accuracy of the downstream task. In RAG, a missed neighbor only hurts if it was the document that contained the answer. Tune recall against your end-to-end eval (see The Evaluation Problem & Benchmark Landscape), not in a vacuum — recall@10 = 0.85 may be plenty.

The knobs, summarized¶

The mental model: HNSW spends memory to buy recall and latency. PQ spends recall to buy memory. IVF spends recall to buy latency. Real systems stack them — IVF65536_HNSW32,PQ64 is a real FAISS factory string that uses an HNSW graph as the coarse quantizer over IVF cells whose contents are PQ-compressed, hitting all three corners at once.

Vector Databases: From Algorithm to System¶

An ANN algorithm (FAISS, ScaNN, hnswlib) is a library. A vector database wraps that algorithm with the boring-but-essential machinery that makes it a production service: persistence, replication, CRUD with consistency, metadata filtering, multi-tenancy, sharding across machines, and horizontal scaling. Knowing the landscape matters for system-design interviews.

• FAISS (Facebook AI Similarity Search) — the canonical library, not a database. It implements every index above (Flat, IVF, PQ, HNSW, and combinations) with CPU and GPU backends. If you build retrieval infra, you will use or reimplement FAISS. It has no server, no persistence layer, no metadata filtering — you bring those.
• Milvus — a distributed vector database built on top of FAISS/HNSW/etc., with sharding, a write-ahead log, object-storage persistence, and a query/data/index node separation for independent scaling.
• Qdrant — Rust, HNSW-based, with strong payload filtering (combine vector search with structured predicates) and good single-node ergonomics.
• Weaviate — HNSW-based, with hybrid (vector + keyword) search and a built-in module ecosystem.
• pgvector — a PostgreSQL extension adding vector, halfvec (float16) and bit column types with IVF-Flat and HNSW indexes. The pragmatic winner when you already run Postgres and your \(N\) is modest: no new system to operate, transactional consistency for free, trivial joins between vectors and your relational data.
• Chroma and LanceDB — embedded, single-process stores you pip install and point at a directory, with no server to run. Chroma is the default in many LangChain/LlamaIndex tutorials; LanceDB stores vectors in the columnar Lance format and serves IVF-PQ indexes straight off disk. For a laptop-scale corpus — including the capstone’s agent — an embedded store plus a FAISS index is the entire retrieval tier.

The filtering problem (where naive ANN breaks)¶

Real queries are rarely “find similar vectors.” They are “find similar vectors where tenant_id = 42 and published_at > 2025-01-01.” This filtered ANN problem is genuinely hard, because the graph/IVF structure was built over all vectors, and the nearest neighbors may all be filtered out, forcing the search to wander. The two strategies — and their failure modes — are interview-classic:

• Pre-filter: compute the matching id set first, then search only within it. Correct, but if the filter is selective, the matching set may be scattered across the whole HNSW graph; you can lose graph connectivity and recall craters, or you fall back to brute force over the filtered set.
• Post-filter: run unfiltered ANN, then drop results that fail the predicate. Fast, but if the filter is selective you may retrieve \(k=10\) candidates and have zero survive — you have to over-fetch by an unknown factor.

Modern systems (Qdrant, Milvus) use filtered HNSW variants that consult the predicate during graph traversal, skipping non-matching nodes while still using them as bridges for connectivity. This is the production frontier and a great thing to mention in a design round.

# Minimal real-library usage: FAISS HNSW with a re-ranking refine stage.
# pip install faiss-cpu
import faiss, numpy as np

d, N = 768, 200_000
rng = np.random.default_rng(0)
xb = rng.standard_normal((N, d)).astype("float32")
faiss.normalize_L2(xb)                  # cosine == inner product on unit vectors

# Build: HNSW with M=32 links, then wrap with exact refinement.
index = faiss.IndexHNSWFlat(d, 32, faiss.METRIC_INNER_PRODUCT)
index.hnsw.efConstruction = 200         # higher => better graph, slower build
index.add(xb)
index.hnsw.efSearch = 64                # the recall-latency knob at query time

xq = rng.standard_normal((5, d)).astype("float32"); faiss.normalize_L2(xq)
D, I = index.search(xq, k=10)           # I: (5,10) ids,  D: (5,10) similarities
print(I.shape, D[0][:3])                # nearest-neighbor ids + scores for query 0

The FAISS factory string: composing an index in one line¶

Everything in this chapter composes, and FAISS exposes that composition as a small string language. index_factory(d, "OPQ16_128,IVF256,PQ16,RFlat") reads left to right as a pipeline: a learned rotation, then a coarse partition, then 16-byte PQ codes inside each cell, then an exact re-ranking stage. This is how real indexes are specified — in production you tune the string and the two runtime parameters, not your code.

# pip install faiss-cpu
import faiss, numpy as np

d, N = 128, 50_000
rng = np.random.default_rng(0)
xb = rng.standard_normal((N, d)).astype("float32"); faiss.normalize_L2(xb)
xq = rng.standard_normal((4, d)).astype("float32"); faiss.normalize_L2(xq)

# Read the string left to right as a pipeline:
#   OPQ16_128 -> learned rotation, output dim 128, tuned for m=16 subquantizers
#   IVF256    -> 256 coarse k-means cells        (keep >= ~39 train points/cell,
#   PQ16      -> 16-byte PQ codes inside a cell   or FAISS prints a warning)
#   RFlat     -> keep the full vectors and re-rank the shortlist exactly
index = faiss.index_factory(d, "OPQ16_128,IVF256,PQ16,RFlat",
                            faiss.METRIC_INNER_PRODUCT)
index.train(xb)      # learns rotation + coarse centroids + PQ codebooks (~1 min)
index.add(xb)        # rotates, assigns to a cell, encodes to 16 bytes

# set_index_parameter reaches THROUGH the RFlat/OPQ wrappers to the IVF inside;
# plain `index.nprobe = 32` fails here because `index` is not the IVF itself.
faiss.ParameterSpace().set_index_parameter(index, "nprobe", 32)
faiss.ParameterSpace().set_index_parameter(index, "k_factor_rf", 10)  # 10*k shortlist
D, I = index.search(xq, 10)

# Persistence: a FAISS index is just a file. This is what a vector DB wraps.
faiss.write_index(index, "/tmp/chunks.faiss")
reloaded = faiss.read_index("/tmp/chunks.faiss")
print(reloaded.ntotal, I.shape)         # 50000 (4, 10)

Note what the last two lines imply. FAISS gives you an index and a file; a vector database is what you get after someone adds a write-ahead log so a crash mid-insert is recoverable, a shard router so the index can exceed one machine’s RAM, a replica set so QPS can exceed one machine’s cores, and a metadata store so you can filter. Those are database problems, not ANN problems — which is why the honest build-versus-buy question is “do I need a database, or do I need an index?”

pgvector: the same algorithms, inside your existing Postgres¶

If your corpus is in the millions rather than the billions and you already run Postgres, the operationally cheapest vector store is a column. pgvector implements the same HNSW and IVF-Flat described above, and — the real win — the metadata filter is just a WHERE clause the planner handles:

CREATE EXTENSION IF NOT EXISTS vector;
CREATE TABLE chunks (
  id        bigserial PRIMARY KEY,
  doc_id    int,
  body      text,
  embedding vector(768)          -- halfvec(768) instead halves the storage
);

-- Bulk-load rows FIRST, then build the index: it is far faster than
-- inserting into a live HNSW graph one row at a time.
CREATE INDEX ON chunks USING hnsw (embedding vector_cosine_ops)
  WITH (m = 16, ef_construction = 64);

SET hnsw.ef_search = 100;         -- the per-session recall/latency knob

-- <=> is cosine distance, <-> is L2, <#> is negative inner product.
-- Use the operator matching your opclass or the index will not be used.
SELECT id, body, 1 - (embedding <=> $1) AS cosine_sim
FROM chunks
WHERE doc_id = 42                 -- filter and vector search in one query
ORDER BY embedding <=> $1
LIMIT 10;

The filtering caveat from above applies here too: with a selective WHERE, an HNSW scan that returns ef_search candidates may yield fewer than LIMIT survivors. Recent pgvector releases address this with iterative index scans (hnsw.iterative_scan), which re-enter the index for more candidates instead of silently under-returning; check your installed version’s documentation, because the behaviour changed across releases.

Putting It Together: A Tiny End-to-End Benchmark¶

To cement the tradeoffs, here is a self-contained harness that builds Flat, IVF, PQ, and HNSW over the same data and reports each method’s recall@10 against the exact baseline. Running it is the fastest way to feel the triangle.

import numpy as np, time

def evaluate(index_search_fn, data, queries, k=10, ground_truth=None):
    """Generic recall@k evaluator against a brute-force oracle."""
    if ground_truth is None:
        ground_truth = []
        for q in queries:
            dd = np.einsum("nd,nd->n", data - q, data - q)
            ground_truth.append(set(np.argpartition(dd, k)[:k].tolist()))
    recalls, t0 = [], time.perf_counter()
    for q, gt in zip(queries, ground_truth):
        pred = set(int(i) for i in index_search_fn(q, k))
        recalls.append(len(gt & pred) / k)
    dt = (time.perf_counter() - t0) / len(queries) * 1e3   # ms/query
    return float(np.mean(recalls)), dt, ground_truth

if __name__ == "__main__":
    rng = np.random.default_rng(7)
    data    = rng.standard_normal((20_000, 64)).astype("float32")
    queries = rng.standard_normal((200, 64)).astype("float32")

    # Flat oracle (recall = 1.0 by construction):
    def flat(q, k):
        dd = np.einsum("nd,nd->n", data - q, data - q)
        return np.argpartition(dd, k)[:k]
    r, ms, gt = evaluate(flat, data, queries)
    print(f"Flat     recall={r:.3f}  {ms:6.2f} ms/q")

    # IVF and HNSW (classes defined earlier in the chapter):
    ivf = IVFFlat(n_list=256).train(data)
    for npb in (4, 16, 64):
        r, ms, _ = evaluate(lambda q, k: ivf.search(q, k, n_probe=npb),
                            data, queries, ground_truth=gt)
        print(f"IVF np={npb:<3d} recall={r:.3f}  {ms:6.2f} ms/q")

    hnsw = HNSW(data, M=16, ef_construction=100)
    for i in range(len(data)):
        hnsw.insert(i)
    for ef in (16, 64, 200):
        r, ms, _ = evaluate(lambda q, k: hnsw.search(q, k, ef=ef),
                            data, queries, ground_truth=gt)
        print(f"HNSW ef={ef:<3d} recall={r:.3f}  {ms:6.2f} ms/q")

The shape of the output is always the same story: Flat is exact but its latency grows with \(N\); IVF and HNSW each expose one knob (n_probe, ef) that buys recall with latency; and if you swapped the data store to PQ codes you would watch memory drop while recall takes a measured hit. That is the entire chapter, in one runnable file.

Where this lands in the capstone. Stack-100M’s narrow auto-research agent retrieves over a corpus of roughly \(10^5\)–\(10^6\) chunks — three to four orders of magnitude below where ANN becomes necessary. The right index there is IndexFlatIP on unit-normalized vectors: exact, zero-parameter, no training, no recall to measure, a few milliseconds per query, and a few hundred megabytes of RAM. Resist the urge to reach for HNSW; the whole point of this chapter’s first section is knowing when brute force is the answer. See A Narrow Auto-Research Agent: ReAct, Tool-Use & Retrieval by Distillation for how that index is wired into the agent’s tool loop, and Evaluation & Serving: Honest Benchmarks, int4 Quantization, and Running on a Laptop for the laptop-scale latency budget it has to fit inside.

For how these retrievers slot into a generation pipeline — query construction, multi-vector retrieval, fusion with keyword search — continue to Retrieval-Augmented Generation Architectures and Chunking, Reranking & Hybrid Search. For the embeddings that feed all of this, revisit Embeddings & Representation Learning. And because retrieval lives in the latency budget of an LLM call, the serving concerns in Inference Economics: Latency, Throughput & Cost apply directly.

Further reading¶

• Hervé Jégou, Matthijs Douze, Cordelia Schmid, Product Quantization for Nearest Neighbor Search (IEEE TPAMI, 2011) — the original PQ and ADC.
• Tiezheng Ge, Kaiming He, Qifa Ke, Jian Sun, Optimized Product Quantization (CVPR 2013 / IEEE TPAMI 2014) — the learned rotation that PQ’s contiguous split needs.
• Moses Charikar, Similarity Estimation Techniques from Rounding Algorithms (STOC 2002) — signed random projections, the cosine LSH and the ancestor of binary quantization.
• Yury Malkov, Dmitry Yashunin, Efficient and Robust Approximate Nearest Neighbor Search Using Hierarchical Navigable Small World Graphs (2016/2018) — HNSW.
• Ruiqi Guo et al., Accelerating Large-Scale Inference with Anisotropic Vector Quantization (ICML 2020) — ScaNN.
• Kevin Beyer, Jonathan Goldstein, Raghu Ramakrishnan, Uri Shaft, When Is “Nearest Neighbor” Meaningful? (ICDT 1999) — the curse of dimensionality, formalized.
• Jeff Johnson, Matthijs Douze, Hervé Jégou, Billion-Scale Similarity Search with GPUs (2017) — the FAISS paper.
• Jon Kleinberg, Navigation in a Small World (Nature, 2000) — the navigable small-world result underpinning graph ANN.
• The FAISS library (facebookresearch/faiss) and its wiki — the canonical reference implementation and index “factory” documentation.

Exercises¶

1. (Conceptual) The chapter claims that once vectors are L2-normalized, a single index serves simultaneously as a cosine index and a Euclidean index. Starting from the identity \(\|\mathbf{q} - \mathbf{x}\|_2^2 = 2 - 2\,\mathbf{q}^\top\mathbf{x}\) for unit vectors, show that ranking the corpus by ascending squared-L2 distance produces the exact same order as ranking by descending cosine similarity. Then explain what the “metric mismatch footgun” from the chapter’s warning box would do to your results, and describe the one-line sanity check the chapter recommends.

2. (Quantitative) You deploy an IVF-Flat index over \(N = 10^{8}\) vectors, choosing \(n_{\text{list}} = \sqrt{N}\) centroids as the chapter’s heuristic suggests, and search with \(n_{\text{probe}} = 16\). (a) Roughly how many database vectors does one query scan? (b) What speedup is that over exact brute force? © The heuristic keeps \(n_{\text{list}} = \sqrt{N}\); as \(N\) grows, does the speedup at fixed \(n_{\text{probe}}\) get better or worse, and why?

3. (Quantitative) You must fit \(N = 10^{8}\) embeddings of dimension \(d = 1024\) into RAM using Product Quantization with \(m = 32\) subquantizers and \(k^{*} = 256\) centroids per subspace. (a) How many bytes does one vector’s PQ code take, and what is the compression ratio versus float32? (b) How many distinct vectors can the codebook implicitly represent? © How much RAM do the PQ codes for all \(N\) vectors take, and how much do the codebooks themselves take? (d) What size is the per-query ADC lookup table, and what operations does scoring one database vector cost?

4. (Implementation) The chapter mentions a re-ranking stage (IndexRefineFlat) that “fetches the top candidates’ full vectors and re-scores them exactly to recover the recall that quantization lost.” Extend the ProductQuantizer class with a method search_rerank(self, q, codes, data, k=10, rerank=100) that (1) uses fast ADC to shortlist the rerank best candidates, then (2) re-scores only those against their full-precision vectors in data with exact squared-L2 and returns the true top-\(k\). Explain the cost you added and why it helps.

def search_rerank(self, q, codes, data, k=10, rerank=100):
    """Two-stage: ADC shortlist -> exact re-rank on full vectors."""
    rerank = min(rerank, codes.shape[0])
    # --- Stage 1: ADC (same LUT trick as `search`), cheap over ALL codes ---
    lut = np.empty((self.m, self.k_star), np.float32)
    for j in range(self.m):
        qsub = q[j*self.dsub:(j+1)*self.dsub]
        diff = self.codebooks[j] - qsub
        lut[j] = np.einsum("kd,kd->k", diff, diff)
    approx = np.zeros(codes.shape[0], np.float32)
    for j in range(self.m):
        approx += lut[j, codes[:, j]]
    # keep the `rerank` smallest approximate distances (O(N), no full sort)
    cand = np.argpartition(approx, rerank - 1)[:rerank]
    # --- Stage 2: exact squared-L2 on the shortlist's FULL vectors ---
    diffs = data[cand] - q                       # (rerank, d)
    exact = np.einsum("nd,nd->n", diffs, diffs)  # (rerank,)
    order = np.argpartition(exact, min(k, rerank - 1))[:k]
    order = order[np.argsort(exact[order])]
    return cand[order]                           # exact top-k among shortlist

5. (Quantitative) HNSW assigns each node a maximum layer by \(\ell = \lfloor -\ln(U) \cdot m_L \rfloor\) with \(U \sim \text{Uniform}(0,1)\) and \(m_L = 1/\ln(M)\). Take \(M = 16\). (a) Derive \(P(\ell \geq L)\) as a function of \(M\) and \(L\). (b) For a graph of \(N = 1{,}000{,}000\) nodes, how many nodes do you expect to appear in layer \(\geq 1\), layer \(\geq 2\), and layer \(\geq 3\)? © Explain in one or two sentences why this geometric thinning is what gives HNSW its logarithmic search cost.

6. (Conceptual) A tenant runs the query “documents similar to \(\mathbf{q}\) where tenant_id = 42”, and tenant_id = 42 matches only \(0.1\%\) of the corpus. Contrast the pre-filter and post-filter strategies for this filtered ANN query: give the failure mode of each under such a selective predicate, and explain why a filter-aware HNSW traversal (as in Qdrant/Milvus) avoids both.