Brute force, IVF and HNSW, written from scratch — and then actually measured.
Three vector search algorithms in readable Python, with a benchmark that answers the question the marketing material never does: at the recall you need, on the data you have, which one is actually faster?
The results were not what I expected, and the interesting part of this repository is the two places the obvious answer turned out to be wrong.
Everyone reaches for HNSW. On 20,000 vectors of dimension 128, it was beaten at every single operating point — not by a rival graph index, but by IVF, and at high recall by plain brute force:
| index | recall@10 | queries/sec | vs. exact |
|---|---|---|---|
| IVF nlist=256, nprobe=4 | 0.956 | 39,555 | 6.3× |
| IVF nlist=256, nprobe=8 | 0.988 | 24,607 | 3.9× |
| IVF nlist=64, nprobe=2 | 0.977 | 19,456 | 3.1× |
| flat (exact) | 1.000 | 6,273 | 1.0× |
| HNSW M=16, ef=10 | 0.964 | 6,590 | 1.05× |
| HNSW M=16, ef=50 | 0.994 | 3,076 | 0.49× |
HNSW at 99.4% recall runs at half the speed of scanning everything.
This is not a bug in the implementation, and it is not a claim that HNSW is overrated. It is the constant factor, and it is worth being precise about why:
- IVF's inner loop is one
numpymatmul over a contiguous slice of the corpus. It runs in BLAS, multi-threaded, at near peak FLOPs, with sequential memory access. - HNSW's inner loop is Python-level graph traversal — dictionary lookups, heap operations, pointer chasing. It performs far fewer distance computations and takes far longer to do them.
HNSW does asymptotically less work. It just does that work in the slowest possible way here. Ported to SIMD C++ over a flat adjacency array — which is what every production implementation actually is — the ranking inverts.
Measured rather than assumed. make crossover grows the corpus with HNSW held at
a fixed configuration, so the trend is the corpus growing rather than the
parameters moving:
| corpus | flat qps | HNSW qps | HNSW / flat | HNSW recall | HNSW build |
|---|---|---|---|---|---|
| 1,000 | 109,029 | 3,292 | 0.03× | 1.000 | 5s |
| 2,000 | 62,446 | 3,673 | 0.06× | 1.000 | 11s |
| 5,000 | 23,021 | 3,841 | 0.17× | 1.000 | 29s |
| 10,000 | 9,882 | 4,000 | 0.40× | 1.000 | 57s |
| 20,000 | 6,979 | 3,894 | 0.56× | 0.994 | 113s |
| 50,000 | 2,216 | 2,155 | 0.97× | 0.970 | 294s |
| 100,000 | 1,173 | 1,935 | 1.65× | 0.979 | 597s |
The crossover is at roughly 50,000 vectors. Below it, brute force wins; above it, HNSW pulls away.
The shapes are exactly what the theory predicts, which is the point. Flat throughput falls as 1/n — 109,029 → 1,173 across a 100× corpus is a 93× drop, near-perfectly linear. HNSW throughput decays far more slowly (3,292 → 1,935 over the same 100×), because a graph descent is logarithmic in the corpus. Two curves with different exponents cross exactly once, and this is where.
Note what that means in practice: for a corpus of 20,000 embeddings, adding an HNSW index would have made the system almost twice as slow while costing two minutes of build time and a large pile of new code. For this implementation, on this machine — a SIMD C++ HNSW moves the crossover down by more than an order of magnitude, which is precisely why the constant factor deserves a benchmark rather than an assumption.
The lesson generalises past vector search: an algorithm's complexity class tells you where the curves eventually cross, and nothing whatsoever about whether you are past that point. Below the crossover, the constant factor is the performance. Benchmark before you adopt.
Building an HNSW graph means choosing which M neighbours each node keeps.
Keeping the M closest is the obvious choice, and it is catastrophically
wrong. Run python bench/ablation.py:
selection ef_search recall@10 qps build s
------------------------------------------------------
heuristic 10 0.9946 9953 28.3
heuristic 100 1.0000 2098 28.3
closest 10 0.6290 11645 3.3
closest 100 0.6352 3538 3.3
0.995 versus 0.629. And notice that raising ef_search tenfold moves the
broken version from 0.629 to 0.635 — the missing recall is not insufficient
search effort. You cannot fix it by looking harder.
The reason: inside a cluster, every node's nearest neighbours are other members of that same cluster. Keep only the closest and the graph becomes dense cliques with nothing bridging them — greedy search enters a cluster and can never leave. The real heuristic keeps a candidate only if it is closer to the node than to any neighbour already kept, which deliberately preserves long-range edges.
The part that would have cost me a week in production:
heuristic mean_degree=32.0 max=32 isolated=0 layers=4
closest mean_degree=32.0 max=32 isolated=0 layers=4
Identical. Fully connected, no isolated nodes, correct layer structure. Every health metric you would think to check says the broken index is fine. Connectivity is not navigability, and only a recall measurement finds the difference. This is why the ablation runs in CI.
The heuristic also costs 8.5× the build time (28.3s vs 3.3s). Worth it, but it should be a decision rather than a surprise.
flat.py — exhaustive search, blocked. A much stronger baseline than its
reputation: one large matrix multiply, in BLAS, perfectly sequential. The
blocking exists to bound peak memory, because a naive queries @ corpus.T
allocates a matrix that exists only to have its top-10 taken.
ivf.py — partition into nlist cells with k-means, scan only the nprobe
most promising. nprobe is a genuine dial from approximate to exact: at
nprobe == nlist it scans everything, so recall is 1.0 by construction — which
is asserted in a test. The failure mode is queries near a cell boundary, whose
true neighbours sit just across the line in a cell that ranked too low to probe.
kmeans.py — Lloyd's algorithm with k-means++ seeding, written out because
the seeding is what matters for IVF and what a library call hides. On power-law
clustered data, random initialisation reliably strands several centroids in the
same dense region, and an overloaded partition is exactly what destroys IVF's
latency guarantee.
hnsw.py — the multi-layer proximity graph. Skip lists, with distance in
place of ordering. Measured layer occupancy on 5,000 vectors: 5000 → 321 → 22 → 2,
the ~M-fold geometric decay that keeps the descent logarithmic.
Three rules the harness enforces, because each is easy to violate by accident and each one flatters the approximate methods:
- Every index answers the same queries against the same corpus, with ground truth computed exhaustively once, up front.
- Build time is reported separately from query time. An index that takes two minutes to build and answers instantly is a different product from one that builds instantly, and a single "speed" number hides that.
- A warm-up pass runs before timing. The first query pays for page faults and BLAS thread-pool spin-up.
The data is generated, not downloaded, and that is a deliberate choice.
Uniform random vectors are the wrong benchmark: in high dimensions they are all
roughly equidistant, so every index looks equally bad and nearest-neighbour
search stops meaning anything. Real embeddings sit on a low-dimensional manifold
in clusters of wildly uneven size. So the generator builds clustered vectors with
a controllable intrinsic dimensionality — generated in intrinsic_d
dimensions, then rotated into d by an orthonormal basis — because intrinsic
dimensionality is what ANN difficulty actually depends on. Cluster sizes follow a
power law, since uniform cluster sizes are the friendliest possible case for IVF
and no real corpus looks like that.
That the knob does something is itself a test:
def test_lower_intrinsic_dimension_makes_an_easier_index():
assert scores[4] > scores[64] # same ambient d, different difficulty61 tests, 97% coverage. The parametrised tests are a shared contract every index must satisfy — k valid ids per query, ordered nearest-first, monotone improvement as the accuracy dial turns — and the rest pin the properties specific to each algorithm:
test_flat_is_exact_by_construction
test_ivf_probing_every_cell_is_exact
test_ivf_partitions_cover_the_corpus_exactly_once
test_hnsw_graph_is_layered_with_geometric_decay
test_hnsw_respects_its_degree_budget
test_hnsw_leaves_no_node_unreachable
test_the_diversity_heuristic_beats_keeping_the_closest
test_ground_truth_matches_a_naive_computation
test_the_warmup_pass_runs_before_timing
That last-but-one is Result 2 as an assertion. The last one checks the ground truth — which everything else is scored against — against the most obvious possible implementation, rather than against itself.
CI lints, tests on 3.11 and 3.12, then re-runs the ablation and asserts the
exact-search invariants (flat recall is 1.0; IVF at nprobe == nlist agrees
with brute force).
make install && make test
make bench # recall/throughput Pareto -> docs/pareto.png
make ablation # the neighbour-selection experiment
make crossover # brute force vs HNSW as the corpus grows- Not a FAISS competitor. FAISS is SIMD C++ with quantisation and GPU kernels. This is legible Python that measures itself honestly. If you need a vector index in production, use FAISS or a vector database — and the numbers above are part of the argument for why.
- No quantisation. PQ and scalar quantisation are where the big memory wins live, and they are a separate project.
- In-memory, single node, no deletes. Deletion in a graph index is genuinely hard — tombstones degrade the graph and repair means rebuilding neighbourhoods.
- Pure-Python HNSW build is slow: ~2 minutes for 20,000 vectors. The build is the honest cost of the implementation language, and it is reported rather than hidden.
MIT. Built by Denina Vincent.