P02 hands-on — Approximate nearest neighbours, block by block

Why a random graph fails, why a navigable one works, and what recall costs.

Source: handson/h02_ann.py --- run it with python3 handson/h02_ann.py
Full project spec: P02 — Approximate Nearest-Neighbour Index

This is the project where the baseline teaches more than the algorithm. Brute force is trivially correct and its cost is knowable in advance; everything an index does is buy latency by giving up recall, so the only interesting question is the exchange rate.

The file builds that comparison honestly. It measures relative contrast on the data before touching an index, because a dataset with low contrast has no nearest neighbour worth finding and no graph will rescue it. It then builds a deliberately random graph to show that connectivity alone is worthless, and only then the navigable small-world construction that makes greedy descent work.

The two-factor model in the assembly is the point of the whole exercise: it predicts the speedup from first principles and lands on the measured number.

Contents

• Block 1 — Vectors and the metric decision
• Block 2 — Relative contrast
• Block 3 — Brute force = the oracle
• Block 4 — recall@k
• Block 5 — A random graph, and why it fails
• Block 6 — NSW: edges that mean something
• Block 7 — The efSearch knob
• The assembly
• The design space
• Latency, bandwidth and the memory hierarchy
• Hardware: CPU, GPU, SSD
• Advanced algorithms and data structures
• How this connects to the rest of the track
• Failure modes at scale
• Primary sources
• Running it
• What to do with this

How to read this page

Each block below is a self-contained lego piece: it builds one mechanism, proves it works on its own, and returns what the next block needs. The code is the real source, sliced out of the script. The output underneath it is the real output, captured by running that script --- not transcribed, not idealised. Where a measurement contradicted what I expected, the contradiction is in the output and the prose says so.

The assembly at the end wires every block into one working thing and measures it.

Block 1 — Vectors and the metric decision

Teaches: normalise once, and three metrics collapse into one

The problem. Before any index exists there is a decision that silently constrains everything after it: which metric, and in what representation. Get it wrong and every recall number you produce afterwards measures the wrong question.

@block(1, "Vectors and the metric decision", "normalise once, and three metrics collapse into one")
def b1(s, show):
    def norm(x): return x / np.linalg.norm(x, axis=-1, keepdims=True)
    n, d = 8000, 48
    data = norm(rng.standard_normal((n, d), dtype=np.float32))
    q = norm(rng.standard_normal((64, d), dtype=np.float32))
    if show:
        dot = data @ q[0]
        l2 = np.linalg.norm(data - q[0], axis=1)
        ident = np.abs(l2**2 - (2 - 2*dot)).max()
        print(f"  {n} vectors, d={d}, unit norm")
        print(f"  ||a-b||^2 == 2-2<a,b>: max deviation {ident:.2e}")
        print(f"  argsort by -dot == argsort by L2: "
              f"{np.array_equal(np.argsort(-dot)[:20], np.argsort(l2)[:20])}")
        print("  so ONE distance function serves cosine, dot and L2 -- but only")
        print("  because we normalised. Skip it and rankings silently diverge.")
    return {"data": data, "q": q, "norm": norm, "n": n, "d": d}

Reading the implementation

The normalisation is not a preprocessing convenience — it is what collapses three different search problems into one. For unit vectors \(|a|=|b|=1\):

\[ |a-b|^2 = |a|^2 + |b|^2 - 2a!\cdot!b = 2 - 2a!\cdot!b \]

so Euclidean distance is a strictly decreasing function of inner product, and cosine similarity is the inner product. Ranking by any of the three gives an identical order (proofs.md P17). That is why the index can store one thing and answer all three questions.

The trap is that this equivalence holds only under normalisation. Maximum inner product search (MIPS) on un-normalised vectors is a genuinely different problem: it has no triangle inequality, so tree and graph methods lose their correctness arguments, and the standard fix is an asymmetric transformation that lifts vectors into \(d+1\) dimensions to restore a metric. If your embeddings carry magnitude information you actually want (document length, confidence), normalising discards it and you must use MIPS properly rather than pretending.

What the numbers say

Output:

  8000 vectors, d=48, unit norm
  ||a-b||^2 == 2-2<a,b>: max deviation 7.15e-07
  argsort by -dot == argsort by L2: True
  so ONE distance function serves cosine, dot and L2 -- but only
  because we normalised. Skip it and rankings silently diverge.

Beyond the toy

Float32 is the default and is usually wasteful. Distance rankings are robust to substantial quantisation — that is the entire premise of PQ (P02 deep dive) — so production systems store float16 (2× saving, negligible recall change) or 8-bit scalar quantisation (4×) for the resident copy and keep float32 only for re-ranking. Since this workload is memory-bound at ~0.25 FLOP/byte, halving the bytes very nearly halves the latency, which makes precision the highest-leverage knob in the system before any algorithmic change.

Block 2 — Relative contrast

Teaches: measure how hard your dataset is BEFORE benchmarking an index

The problem. Every ANN paper reports recall on a benchmark dataset. Almost none report whether the dataset had a findable nearest neighbour in the first place. Relative contrast is the diagnostic that tells you, and it costs one pass over the data.

@block(2, "Relative contrast", "measure how hard your dataset is BEFORE benchmarking an index")
def b2(s, show):
    data, q = s["data"], s["q"]
    def rc(data, q, k=10):
        dist = np.sqrt(np.maximum(0, 2 - 2 * (data @ q.T)))
        srt = np.sort(dist, axis=0)
        return float(np.mean(dist.mean(axis=0) / srt[0]))
    r = rc(data, q)
    if show:
        print(f"  RC = mean distance / nearest distance = {r:.3f}")
        for dd in (8, 48, 256):
            sub = s["norm"](rng.standard_normal((3000, dd), dtype=np.float32))
            sq = s["norm"](rng.standard_normal((32, dd), dtype=np.float32))
            print(f"    d={dd:>4}: RC={rc(sub, sq):.3f}")
        print("  RC falls toward 1 with dimension: every point equidistant, no")
        print("  gradient for greedy search. RC predicts difficulty; d does not.")
    return {"rc": r}

Reading the implementation

Relative contrast is the ratio of the mean distance to the nearest distance:

\[ C = \frac{\bar{d}}{d_{\min}} \]

As \(C \to 1\), every point is roughly equidistant from the query and there is no nearest neighbour to find — not because the algorithm is weak but because the data carries no signal at that dimensionality. He, Kumar and Chang showed that the difficulty of nearest-neighbour search is essentially a function of this number, and that it degrades with dimension for any distribution whose components are independent.

The implementation detail that matters: compute this on the full data with the actual metric. An earlier version of this file computed contrast on a 2000-point subsample using similarity ratios and returned ≈1.30 for every configuration — a broken metric that looked plausible. The fix was to use proper Euclidean distances over the whole set, and it changed the conclusions.

What the numbers say

Output:

  RC = mean distance / nearest distance = 1.439
    d=   8: RC=3.693
    d=  48: RC=1.410
    d= 256: RC=1.130
  RC falls toward 1 with dimension: every point equidistant, no
  gradient for greedy search. RC predicts difficulty; d does not.

Beyond the toy

• Intrinsic dimension is the companion diagnostic. Real embeddings from a trained model occupy a manifold of far lower dimension than their ambient size — a 768-dim CLIP embedding might have intrinsic dimension 20--40 — which is precisely why ANN works at all on real data and fails on uniform random data of the same width. Estimate it with the two-NN or MLE estimators before concluding an index is broken.
• The synthetic-data trap. The clustered generator in this file originally produced uniform data because the cluster spread \(\sigma\sqrt{d} = 2.0\) swamped unit-norm centres, giving contrast 1.393 against uniform's 1.356. Any synthetic benchmark needs this check, and the file now warns when \(\sigma\sqrt{d} > 0.7\).
• Practical use. If contrast is near 1 on your production queries, the correct action is not to tune ef — it is to fix the embedding, add a re-ranking stage with a richer model, or accept that top-k is meaningless and return a diversified set instead.

Block 3 — Brute force = the oracle

Teaches: you cannot measure recall without exact ground truth

The problem. Recall is defined against exact ground truth. Without an oracle there is no number, and "it returns plausible results" is not a measurement. This block builds the oracle and, incidentally, the baseline that the index must beat.

@block(3, "Brute force = the oracle", "you cannot measure recall without exact ground truth")
def b3(s, show):
    data, q = s["data"], s["q"]
    def brute(qv, k=10):
        sims = data @ qv
        idx = np.argpartition(-sims, k)[:k]
        return idx[np.argsort(-sims[idx])].tolist()
    t0 = time.perf_counter()
    truth = [brute(v) for v in q]
    ms = (time.perf_counter() - t0) / len(q) * 1e3
    if show:
        naive = sorted(range(s["n"]), key=lambda i: -float(data[i] @ q[0]))[:10]
        print(f"  agrees with a naive sort: {naive == truth[0]}")
        print(f"  {ms:.3f} ms/query, {1000/ms:>6.0f} qps, recall 1.0 BY CONSTRUCTION")
        print(f"  cost model: {s['n']}x{s['d']} = {s['n']*s['d']:,} MACs, "
              f"{s['n']*s['d']*4/1e6:.1f} MB streamed -- one BLAS call, prefetchable")
    return {"brute": brute, "truth": truth, "brute_ms": ms}

Reading the implementation

Brute force is \(O(Nd)\) per query with perfect recall and, crucially, a known cost — you can compute it exactly in advance from the corpus size and dimension, which no approximate method allows. It is also embarrassingly parallel, cache- friendly (a sequential scan, unlike the graph's pointer chasing), and trivially correct.

That combination is why brute force is not a strawman. On a GPU, scanning 1M × 768 fp16 vectors is 1.5 GB at ~3 TB/s ≈ 0.5 ms — faster than most CPU ANN implementations, and exact. The regime where an index is genuinely necessary starts higher than most people assume, and the honest version of "we need a vector database" is usually "we need one above N vectors, and here is N".

What the numbers say

Output:

  agrees with a naive sort: True
  0.129 ms/query,   7776 qps, recall 1.0 BY CONSTRUCTION
  cost model: 8000x48 = 384,000 MACs, 1.5 MB streamed -- one BLAS call, prefetchable

Beyond the toy

The ground truth itself becomes a cost at scale: computing exact top-k for 10k queries against 100M vectors is \(10^{12}\) distance computations. Standard practice is to compute it once on a GPU and ship it with the benchmark, which is what ANN-Benchmarks and the BigANN challenge datasets do. A benchmark that recomputes ground truth with an approximate method — and it happens — measures agreement between two approximations, not recall.

Block 4 — recall@k

Teaches: the contract: what exactly did approximation cost?

The problem. "Recall" is used loosely enough to be meaningless. This block pins it down, because every later number on this page is denominated in it.

@block(4, "recall@k", "the contract: what exactly did approximation cost?")
def b4(s, show):
    def recall_at_k(got, want, k):
        return len(set(got[:k]) & set(want[:k])) / k
    if show:
        print(f"  perfect: {recall_at_k([1,2,3],[1,2,3],3):.2f}   "
              f"two of three: {recall_at_k([1,2,9],[1,2,3],3):.2f}   "
              f"none: {recall_at_k([7,8,9],[1,2,3],3):.2f}")
        print("  always against EXACT ground truth, same metric, same data, stated k")
    return {"recall_at_k": recall_at_k}

Reading the implementation

\[ \text{recall@}k = \frac{|,\text{returned}_k \cap \text{true}_k,|}{k} \]

Three details that separate a usable metric from a misleading one:

• Set intersection, not order. Recall@10 does not care whether the true nearest neighbour came back at rank 1 or rank 10. If order matters for your product, measure NDCG or MRR instead — they are different questions and a system can be excellent at one and poor at the other.
• Ties. With duplicate or near-duplicate vectors, "the true top-k" is not unique and naive recall under-reports. Benchmarks handle this by comparing against a distance threshold rather than an ID set.
• The denominator is \(k\), not the number returned. A method that returns 3 results of which all 3 are correct has recall@10 of 0.3, not 1.0. This is exactly the post-filter failure mode in P03, and defining recall this way is what makes it visible.

What the numbers say

Output:

  perfect: 1.00   two of three: 0.67   none: 0.00
  always against EXACT ground truth, same metric, same data, stated k

Beyond the toy

Recall@k against exact neighbours is an intrinsic metric — it measures the index against itself, not against the task. A retrieval system with recall@10 = 0.85 may be indistinguishable from one at 0.99 in end-to-end product terms, because the ranker downstream reorders anyway (P08) and because the embedding's own error dwarfs the index's. Always pair the intrinsic number with an extrinsic one before spending weeks moving it.

Block 5 — A random graph, and why it fails

Teaches: the naive design, measured, so the fix is motivated

The problem. It is tempting to think a graph search works because the graph is connected. This block builds a connected graph with the right degree and the right number of edges, and shows it produces almost nothing — which is what makes the next block's construction motivated rather than arbitrary.

@block(5, "A random graph, and why it fails", "the naive design, measured, so the fix is motivated")
def b5(s, show):
    data, n = s["data"], s["n"]
    deg = 16
    g = [rng.choice(n, deg, replace=False).tolist() for _ in range(n)]
    def greedy(qv, graph, entry=0, ef=32):
        seen = {entry}; d0 = 1 - float(data[entry] @ qv)
        cand = [(d0, entry)]; res = [(-d0, entry)]; nd = 1
        while cand:
            d, node = heapq.heappop(cand)
            if -res[0][0] < d and len(res) >= ef: break
            for nb in graph[node]:
                if nb in seen: continue
                seen.add(nb); dn = 1 - float(data[nb] @ qv); nd += 1
                if len(res) < ef or dn < -res[0][0]:
                    heapq.heappush(cand, (dn, nb)); heapq.heappush(res, (-dn, nb))
                    if len(res) > ef: heapq.heappop(res)
        return [i for _, i in heapq.nsmallest(10, [(-a, b) for a, b in res])], nd
    rec = np.mean([s["recall_at_k"](greedy(v, g)[0], t, 10)
                   for v, t in zip(s["q"], s["truth"])])
    if show:
        print(f"  {deg} RANDOM edges per node, greedy beam ef=32")
        print(f"  recall@10 = {rec:.4f}   <- near useless, and that is the point")
        print("  a random graph has no locality, so greedy descent has nothing to")
        print("  descend. The fix is not a bigger beam; it is better EDGES.")
    return {"greedy": greedy, "random_recall": float(rec)}

Reading the implementation

greedy is a best-first search with a bounded candidate set: pop the closest unexpanded node, evaluate its neighbours, keep the best ef seen. Two lines carry the algorithm:

• if -res[0][0] < d and len(res) >= ef: break — the termination rule. Stop when the closest unexplored candidate is farther than the worst result already held. This is the standard best-first bound, and it is what makes the search adaptive: easy queries terminate early, hard ones explore more.
• nd += 1 counts distance computations. That counter is the real cost metric — wall-clock time varies with implementation and machine, but distance computations are the invariant that lets you compare an index against brute force's \(N\) honestly. Reporting speedup in terms of nd is what makes P15's capacity model possible.

What the numbers say

Output:

  16 RANDOM edges per node, greedy beam ef=32
  recall@10 = 0.0609   <- near useless, and that is the point
  a random graph has no locality, so greedy descent has nothing to
  descend. The fix is not a bigger beam; it is better EDGES.

Beyond the toy

Greedy descent on a graph is a local algorithm: it only ever moves downhill in distance, so it terminates at a local minimum of the distance function restricted to the graph. A random graph has no correlation between edge structure and geometry, so almost every node is a local minimum and the search stops immediately. The fix is therefore not more search (a bigger ef barely helps — worth measuring) but better edges.

The theoretical frame is Kleinberg's small-world result: a lattice augmented with long-range links whose length distribution follows \(P(u\to v) \propto d(u,v)^{-r}\) is navigable by a greedy decentralised algorithm in \(O(\log^2 n)\) steps only when \(r\) equals the lattice dimension. Too few long links and you cannot cross the space; too many and you cannot home in. NSW's insertion procedure produces that distribution approximately and for free, which is the elegance of the next block.

Block 6 — NSW: edges that mean something

Teaches: insert by searching what you have built so far

The problem. Build a graph whose edges encode geometry, using nothing but the search procedure you already have. The construction is three lines and the reason it works is a genuinely deep result.

@block(6, "NSW: edges that mean something", "insert by searching what you have built so far")
def b6(s, show):
    data, n = s["data"], s["n"]
    M, efC = 12, 60
    graph = [[] for _ in range(n)]
    order = rng.permutation(n)
    entry = int(order[0])
    for count, node in enumerate(order[1:], 1):
        node = int(node)
        found, _ = s["greedy"](data[node], graph, entry, efC)
        cands = sorted(found, key=lambda i: 1 - float(data[i] @ data[node]))[:M]
        graph[node] = cands
        for c in cands:
            graph[c].append(node)
            if len(graph[c]) > 2 * M:               # degree cap
                dd = data[graph[c]] @ data[c]
                graph[c] = [graph[c][i] for i in np.argsort(-dd)[:2 * M]]
    if show:
        edges = sum(len(x) for x in graph)
        print(f"  M={M} efConstruction={efC}, {edges/n:.1f} edges/node")
        print("  early insertions land in a near-empty graph, so their edges are")
        print("  LONG. Those accidental long links are the small-world property.")
    return {"graph": graph, "entry": entry, "M": M}

Reading the implementation

The insertion rule is: to insert a point, search the graph you have built so far, and link to what you find. That single recursive idea produces navigability without any global structure:

• Points inserted early, when the graph is sparse and the search is inaccurate, get long-range links — they connect regions that are far apart.
• Points inserted late, when the graph is dense and the search is accurate, get short-range links — they refine the local neighbourhood.

The result is a length distribution across scales, which is exactly Kleinberg's navigability condition arrived at by accident of construction. Nobody computes it; it falls out of using the search to build the index.

• M is the degree bound and the memory knob: M × 4 bytes per node for the edge list. It is also the recall knob that most benchmarks under-report — sweeping M moves the memory/recall curve far more than sweeping ef moves the latency/recall one.
• Edges are added bidirectionally, and this matters: a directed graph built this way has poor in-degree for early nodes and the search cannot get back out of a region it descends into.
• HNSW's addition on top of this is the layer hierarchy — long links in upper layers, taken first — which reduces the number of distance computations to reach the right region, not the number of hops.

What the numbers say

Output:

  M=12 efConstruction=60, 17.1 edges/node
  early insertions land in a near-empty graph, so their edges are
  LONG. Those accidental long links are the small-world property.

Beyond the toy

• Pruning is the missing piece. Real HNSW does not keep the M nearest neighbours; it applies a heuristic that keeps an edge \(u\to v\) only if no already-selected \(w\) is closer to \(v\) than \(u\) is. That is a relative-neighbourhood-graph relaxation, and it produces diverse edges — one per direction rather than M edges into the same dense cluster. Without it, recall on clustered data is substantially worse.
• Construction cost is \(O(N \log N \cdot M)\) distance computations, which for 100M vectors is hours on many cores. This makes index build a batch job (P06) and makes incremental insertion valuable, which is why the NSW-style construction — which is naturally incremental — beat the tree methods that require a global build.
• The entry point matters more than it looks. A fixed entry point means every search starts in the same place and the nodes near it are traversed on every query — hot in cache, but also a concentration of load. HNSW's top layer solves this; some implementations use multiple random entry points.

Block 7 — The efSearch knob

Teaches: recall is concave in ef; pick the recall FIRST

The problem. ef is the dial between latency and recall, and its shape decides how you configure the system. This block measures the curve rather than assuming it, and the shape has a direct operational consequence.

@block(7, "The efSearch knob", "recall is concave in ef; pick the recall FIRST")
def b7(s, show):
    if show:
        print(f"  {'ef':>5}{'recall@10':>11}{'ms':>9}{'dists/q':>10}{'vs brute':>10}")
        for ef in (10, 32, 64, 128):
            t0 = time.perf_counter(); recs = []; nds = 0
            for v, t in zip(s["q"], s["truth"]):
                got, nd = s["greedy"](v, s["graph"], s["entry"], ef)
                recs.append(s["recall_at_k"](got, t, 10)); nds += nd
            ms = (time.perf_counter() - t0) / len(s["q"]) * 1e3
            print(f"  {ef:>5}{np.mean(recs):>11.4f}{ms:>9.3f}"
                  f"{nds/len(s['q']):>10.0f}{s['brute_ms']/ms:>9.2f}x")
    return {}

Reading the implementation

ef (the search-time beam width) bounds the candidate set: the search keeps the ef best nodes seen and terminates when no unexplored candidate can improve them. Larger ef means more distance computations and a lower chance of getting stuck in a local minimum.

The critical property is that recall is concave in ef — the first doubling buys a lot, the fourth buys almost nothing, while cost grows roughly linearly. The operational consequence is the reverse of how people usually configure it:

Pick the recall you need first, then find the smallest ef that reaches it. Configuring ef to a latency budget and reporting whatever recall results is how systems end up paying 4× the cost for 1% of recall.

What the numbers say

Output:

     ef  recall@10       ms   dists/q  vs brute
     10     0.3344    0.223       263     0.58x
     32     0.6219    0.540       595     0.24x
     64     0.7891    0.875       978     0.15x
    128     0.9328    1.540      1672     0.08x

Beyond the toy

• ef is a per-query dial, not a global one. Since the search terminates adaptively, you can raise ef for queries the system judges hard (low top-1 similarity, high tail-index) and lower it for easy ones. That is real, deployed practice and it buys a better recall-per-millisecond than any global setting.
• The two-factor speedup model in the assembly is what makes this predictable rather than empirical: speedup is (fraction of vectors examined) × (per-vector cost ratio), and predicting 0.15× against a measured 0.15× is what turns a benchmark into a model you can extrapolate with.
• Recall is not the tail. The mean recall at ef=64 hides that some queries return nothing useful. Report the distribution — recall@10 at p10 across queries — because a system with mean recall 0.9 and 5% of queries at recall 0.0 is a very different product from one with uniform 0.9.

The assembly

Every block above, wired together into one working system:

def assembly(s):
    print("\nThe seven blocks are an ANN index. Now the measurement that matters.\n")
    ef = 64
    t0 = time.perf_counter(); recs = []; nds = 0
    for v, t in zip(s["q"], s["truth"]):
        got, nd = s["greedy"](v, s["graph"], s["entry"], ef)
        recs.append(s["recall_at_k"](got, t, 10)); nds += nd
    ms = (time.perf_counter() - t0) / len(s["q"]) * 1e3
    dpq = nds / len(s["q"])
    algo = s["n"] / dpq
    const = (ms * 1e6 / dpq) / (s["brute_ms"] * 1e6 / s["n"])
    print(f"  at efSearch={ef}: recall {np.mean(recs):.4f}, {ms:.3f} ms/query")
    print(f"  random-graph recall was {s['random_recall']:.4f} -> NSW is "
          f"{np.mean(recs)/max(s['random_recall'],1e-9):.0f}x better on the SAME search code.\n")
    print("  THE TWO-FACTOR MODEL -- decompose before you tune anything:")
    print(f"    algorithmic win : {s['n']:,} / {dpq:.0f} distances = {algo:>7.1f}x fewer")
    print(f"    constant factor : {ms*1e6/dpq:>7.1f} ns/dist (python) vs "
          f"{s['brute_ms']*1e6/s['n']:.1f} ns (BLAS) = {const:>6.1f}x slower each")
    print(f"    predicted speedup: {algo:.1f} / {const:.1f} = {algo/const:.2f}x")
    print(f"    measured speedup : {s['brute_ms']/ms:.2f}x")
    print("\n  The algorithm is RIGHT and the implementation is WRONG. Those are")
    print("  different bugs. Only the distance counter can tell you which you have.")
    print(f"\n  Dataset difficulty: RC = {s['rc']:.3f}. Report it with every recall")
    print("  number, or the result does not transfer to anyone else's corpus.")
    print("\n  Built: metric choice -> RC -> oracle -> recall -> random graph ->")
    print("  NSW -> ef sweep -> the decomposition.")
    print("  Missing, and on the project page: HNSW layers (m6), Algorithm 4 (m7),")
    print("  a compiled inner loop (m8), persistence (m9), hnswlib comparison (E12).")

Output:

The seven blocks are an ANN index. Now the measurement that matters.

  at efSearch=64: recall 0.7891, 0.895 ms/query
  random-graph recall was 0.0609 -> NSW is 13x better on the SAME search code.

  THE TWO-FACTOR MODEL -- decompose before you tune anything:
    algorithmic win : 8,000 / 978 distances =     8.2x fewer
    constant factor :   915.3 ns/dist (python) vs 16.1 ns (BLAS) =   56.9x slower each
    predicted speedup: 8.2 / 56.9 = 0.14x
    measured speedup : 0.14x

  The algorithm is RIGHT and the implementation is WRONG. Those are
  different bugs. Only the distance counter can tell you which you have.

  Dataset difficulty: RC = 1.439. Report it with every recall
  number, or the result does not transfer to anyone else's corpus.

  Built: metric choice -> RC -> oracle -> recall -> random graph ->
  NSW -> ef sweep -> the decomposition.
  Missing, and on the project page: HNSW layers (m6), Algorithm 4 (m7),
  a compiled inner loop (m8), persistence (m9), hnswlib comparison (E12).

The design space

Every ANN index is a bet about where you can afford to lose information. The four families make different bets, and they fail differently.

In practice the production answer is a composite: IVF to narrow, PQ to compress, a graph over centroids to route, exact re-ranking on the raw vectors of the top few hundred. FAISS's IVF65536_HNSW32,PQ64 is that sentence as a factory string.

Why the graph works at all

Greedy descent on a proximity graph is a decentralised search: at each step, move to the neighbour closest to the query. It terminates at a local minimum, so the entire design problem is making local minima rare. NSW achieves that by inserting nodes in random order and linking each to its M nearest already-inserted nodes — long-range links form early while the graph is sparse, short-range links late, producing a small-world graph of \(O(\log N)\) diameter. HNSW adds explicit layers so long hops are taken first, cutting the number of distance computations rather than the number of hops.

Block 3's result — a random graph gives 0.06 recall — is the control that proves connectivity is not the active ingredient. Navigability is.

Latency, bandwidth and the memory hierarchy

Graph search is pointer chasing, the worst access pattern the memory hierarchy has. Each hop is a dependent load: you cannot prefetch hop \(k+1\) until hop \(k\) resolves.

Using this machine's measured constants (numbers.md):

A search touching 600 distinct vectors of 128 dims × 4 B reads ~300 KB scattered across memory: ~600 line-missing loads at 121 ns ≈ 73 µs of pure latency against perhaps 5 µs of arithmetic. The search is latency-bound — not bandwidth-bound and nowhere near compute-bound — which is why the roofline efficiency in P15 comes out below 1%.

Three consequences follow directly:

1. Layout dominates. Storing vectors in visit order, or inlining a compressed code beside each edge list, converts random reads into sequential ones. This makes quantisation a latency optimisation, not merely a memory one.
2. TLB reach matters. A 100 GB index with 4 KiB pages needs 25M page-table entries; a 1536-entry TLB covers 6 MB. Huge pages (2 MiB) extend reach 512× and are routinely worth 10--30% — P12 block 2's mechanism applied to a data structure.
3. Batching queries helps for a non-obvious reason. It barely raises arithmetic intensity; what it creates is memory-level parallelism, many independent dependent-load chains in flight at once, which is the only way to hide 121 ns.

Hardware: CPU, GPU, SSD

• CPU is the natural home for graph search: branchy, pointer-heavy, latency-sensitive, helped by a large out-of-order window and hardware prefetchers — which the Sattolo-cycle trap in numbers.md §14 shows are easily fooled into flattering a benchmark.
• GPU is the natural home for brute force. A 1M × 768 fp16 matrix is 1.5 GB; scanning it at ~3 TB/s takes 0.5 ms. Below a few million vectors, GPU brute force beats CPU ANN on latency and is exact — which is most of why FAISS-GPU exists. Graph traversal maps badly to SIMT because each query diverges.
• SSD changes the algorithm, not the parameters. DiskANN/Vamana prune with an α-RNG rule to get a small diameter, keep PQ codes of everything in RAM to guide the walk, and read full vectors for only a few hundred candidates. That design exists because NVMe random reads are ~50 µs and parallel (deep queues), whereas HDD random reads are 10 ms and effectively serial.

Advanced algorithms and data structures

• Product quantisation splits a \(d\)-dim vector into \(m\) subvectors, k-means each subspace to 256 centroids, stores \(m\) bytes. Distance becomes a sum of \(m\) lookups in an L1-resident table. OPQ learns a rotation first, because PQ assumes subspace independence and raw embedding dimensions are not independent.
• ScaNN's anisotropic loss notes that for maximum-inner-product search, quantisation error parallel to the query matters far more than orthogonal error, and weights training accordingly — a rare case where changing the loss, not the structure, buys large recall.
• Relative contrast (He, Kumar & Chang) is what block 2 computes: \(C = \bar{d}/d_{\min}\). As \(C \to 1\) no algorithm can separate the nearest neighbour from the crowd, because the data has no signal. Measure it before blaming an index.
• α-RNG pruning — keep edge \(u \to v\) unless some \(w\) is closer to both — is a relative-neighbourhood-graph relaxation that bounds out-degree while preserving navigability. It is the theoretical spine under both Vamana and HNSW's pruning heuristic.

How this connects to the rest of the track

• P03 wraps this index in a database: filters, persistence, planning.
• P08 uses it as a candidate generator; its block 7 shows the index's recall@C is a hard ceiling on the recommender's recall@k.
• P01's attention is a soft version of the same lookup and meets the same high-dimensional contrast problem.
• P14 explains the ≈0.25 FLOP/byte intensity that makes this memory-bound.
• P04's Bloom filter is the same idea as a PQ code: a cheap approximate test in fast memory that avoids an expensive exact one in slow memory.

Failure modes at scale

• Deletions break connectivity. Tombstoning without repair slowly disconnects regions; most systems rebuild segments instead (P03, P04 compaction).
• Distribution shift after an embedding upgrade invalidates centroids and graph simultaneously. There is no incremental fix, only a rebuild.
• Recall measured on indexed vectors overstates production recall, because real queries are out-of-distribution relative to the corpus.
• The p99 is not p50 × constant. Hop counts have a long tail; a query landing in a sparse region takes 5--10× the median.

Primary sources

• Malkov & Yashunin, Efficient and Robust ANN Search Using HNSW (2016).
• Jégou, Douze & Schmid, Product Quantization for NN Search (2011).
• Subramanya et al., DiskANN (NeurIPS 2019) — the SSD design.
• Guo et al., Anisotropic Vector Quantization (ScaNN, ICML 2020).
• He, Kumar & Chang, On the Difficulty of Nearest Neighbor Search (ICML 2012).
• Johnson, Douze & Jégou, Billion-scale Similarity Search with GPUs (2017).

Running it

python3 handson/h02_ann.py            # every block, then the assembly
python3 handson/h02_ann.py --block 3  # just block 3 and its prerequisites
python3 handson/h02_ann.py --quiet    # the assembly only

What to do with this

The parameter that matters most is the one this file leaves at a default. Sweep M (edges per node) rather than ef and you will find the memory/recall curve that most benchmark tables omit entirely. Then re-run the contrast measurement on a real embedding set --- CLIP or a sentence encoder --- and compare it with the synthetic number here. If the real contrast is lower, every recall figure you have read about that dataset is easier than it looks.

Milestones, experiments, readings and exit criteria for this project: P02 — Approximate Nearest-Neighbour Index.