Every vector database you have ever used — Qdrant, Weaviate, Milvus — relies on the same core algorithm for fast approximate nearest neighbor search: HNSW (Hierarchical Navigable Small World). Most developers treat it as a black box. This post tears it open, implements it in pure Python, and shows you exactly what levers you pull when you tune a production vector index.
Why Not Just Use Brute Force?
For small datasets (under 10k vectors), brute-force k-NN works fine: compute the distance from your query to every stored vector, sort, return top-k. At 1M+ vectors it breaks down. A single 1536-dimensional embedding query requires roughly 1.5 billion float multiplications against the full dataset. HNSW reduces that to a few thousand operations at 95–99% recall.
The accuracy cost is measured by recall@k: the fraction of the true top-k nearest neighbors that the approximate index returns. Production HNSW configurations typically achieve recall@10 > 0.97 at 50–200x the speed of brute force.
How HNSW Works
HNSW builds a layered proximity graph — think of it as a hierarchy of skip lists in high-dimensional space.
-
Layer 0 contains every inserted vector, each connected to its
Mnearest neighbors. - Layer 1 contains a random ~1/e subset of layer 0, similarly connected.
- Higher layers contain progressively sparser subsets, up to a single entry point.
Search starts at the entry point in the highest layer, descends greedily to the closest node, drops one layer, and repeats until layer 0. At layer 0, a beam search with width ef collects the final candidates.
Insert mirrors search: find the neighborhood of the new vector at each layer, then add bidirectional edges.
Two parameters control the quality-speed tradeoff:
-
M— connections per node. Higher = better recall, more memory and build time. -
ef_construction— beam width at index build. Higher = better graph quality, slower inserts.
The Implementation
import heapq
import math
import random
from typing import Optional
class HNSWIndex:
def __init__(self, dim: int, M: int = 16, ef_construction: int = 200):
self.dim = dim
self.M = M
self.M_max0 = 2 * M
self.ef_construction = ef_construction
self.ml = 1.0 / math.log(M)
self.vectors: list[list[float]] = []
self.graphs: list[dict[int, list[int]]] = []
self.entry_point: Optional[int] = None
self.max_layer = -1
def _dist(self, a: list[float], b: list[float]) -> float:
return sum((x - y) ** 2 for x, y in zip(a, b)) ** 0.5
def _random_level(self) -> int:
return int(-math.log(random.random()) * self.ml)
def _search_layer(
self, query: list[float], ep: int, ef: int, layer: int
) -> list[tuple[float, int]]:
visited = {ep}
d = self._dist(query, self.vectors[ep])
candidates = [(d, ep)] # min-heap
found = [(-d, ep)] # max-heap (negated distances)
while candidates:
dist, node = heapq.heappop(candidates)
worst = -found[0][0]
if dist > worst:
break
for nb in self.graphs[layer].get(node, []):
if nb not in visited:
visited.add(nb)
nd = self._dist(query, self.vectors[nb])
if nd < worst or len(found) < ef:
heapq.heappush(candidates, (nd, nb))
heapq.heappush(found, (-nd, nb))
if len(found) > ef:
heapq.heappop(found)
return sorted((-neg_d, node) for neg_d, node in found)
def add(self, vector: list[float]) -> int:
idx = len(self.vectors)
self.vectors.append(vector)
level = self._random_level()
while len(self.graphs) <= level:
self.graphs.append({})
for lyr in range(len(self.graphs)):
self.graphs[lyr].setdefault(idx, [])
if self.entry_point is None:
self.entry_point = idx
self.max_layer = level
return idx
ep = self.entry_point
for lyr in range(self.max_layer, level, -1):
results = self._search_layer(vector, ep, 1, lyr)
if results:
ep = results[0][1]
for lyr in range(min(level, self.max_layer), -1, -1):
neighbors = self._search_layer(vector, ep, self.ef_construction, lyr)
M_max = self.M_max0 if lyr == 0 else self.M
selected = [i for _, i in neighbors[:M_max]]
self.graphs[lyr][idx] = selected
for nb in selected:
self.graphs[lyr][nb].append(idx)
if len(self.graphs[lyr][nb]) > M_max:
nb_vec = self.vectors[nb]
pruned = sorted(
(self._dist(nb_vec, self.vectors[n]), n)
for n in self.graphs[lyr][nb]
)[:M_max]
self.graphs[lyr][nb] = [n for _, n in pruned]
if neighbors:
ep = neighbors[0][1]
if level > self.max_layer:
self.max_layer = level
self.entry_point = idx
return idx
def search(
self, query: list[float], k: int = 10, ef: int = 50
) -> list[tuple[float, int]]:
if self.entry_point is None:
return []
ep = self.entry_point
for lyr in range(self.max_layer, 0, -1):
results = self._search_layer(query, ep, 1, lyr)
if results:
ep = results[0][1]
return self._search_layer(query, ep, max(ef, k), 0)[:k]
The _search_layer function is the heart of the algorithm. It maintains two heaps in parallel: a min-heap of candidates to explore and a max-heap of the ef best results found so far. The early-exit condition (dist > worst) is what makes HNSW fast — once the closest unvisited candidate is farther than the worst result in your best-ef set, you stop.
Measuring Recall Against Brute Force
Never deploy an ANN index without first measuring its actual recall on your data distribution. Here is a simple benchmark harness:
import random
import time
def brute_force_knn(
vectors: list[list[float]], query: list[float], k: int
) -> set[int]:
dists = [
(sum((a - b) ** 2 for a, b in zip(query, v)) ** 0.5, i)
for i, v in enumerate(vectors)
]
return {i for _, i in sorted(dists)[:k]}
def benchmark(n: int = 5_000, dim: int = 64, k: int = 10):
rng = random.Random(42)
data = [[rng.gauss(0, 1) for _ in range(dim)] for _ in range(n)]
index = HNSWIndex(dim=dim, M=16, ef_construction=100)
t0 = time.perf_counter()
for v in data:
index.add(v)
print(f"Build: {time.perf_counter() - t0:.2f}s for {n} vectors")
queries = [[rng.gauss(0, 1) for _ in range(dim)] for _ in range(200)]
recall_total = 0.0
t0 = time.perf_counter()
for q in queries:
hnsw_hits = {i for _, i in index.search(q, k=k, ef=50)}
true_hits = brute_force_knn(data, q, k)
recall_total += len(hnsw_hits & true_hits) / k
avg_ms = (time.perf_counter() - t0) / len(queries) * 1000
print(f"Recall@{k}: {recall_total / len(queries):.3f}")
print(f"Avg query: {avg_ms:.2f}ms")
benchmark()
On a modern laptop with 5k vectors at dim=64, expect ~0.92 recall@10 at under 2ms per query with ef=50. Increase ef to 200 to push recall past 0.98 at 4–6x the latency cost. The tradeoff is tunable at query time without rebuilding the index.
Production Notes
The pure-Python implementation above is for understanding, not production throughput. For real workloads, use hnswlib (C++ backend, Python bindings, ~100x faster) or a dedicated vector database.
Before you commit to a hosted vector search service, review the security posture of your embedding pipeline — embedding user data and storing it in a third-party index has data residency and access control implications. We have a concise checklist covering this in our security hardening resources.
A few practical tuning notes:
- Memory: at M=16, expect ~800 bytes per vector for the graph structure alone. Budget accordingly.
- Dimensionality: embeddings above 512 dimensions benefit from PCA reduction to 256 before indexing — roughly 6x speedup at ~2% recall loss.
- Concurrent access: HNSW search is read-safe for parallelism; inserts require exclusive write locking.
-
Persistence: serialize
vectorsandgraphswith msgpack or pickle. hnswlib has built-insave_index/load_index.
The Takeaway
HNSW is not magic — it is a hierarchical proximity graph where you exchange a small fraction of accuracy for orders-of-magnitude speed gains. Implementing it from scratch makes the parameters M, ef_construction, and ef concrete rather than abstract knobs. When your vector similarity queries return unexpected results in production, you will know immediately whether to increase ef, rebuild with higher M, check your distance metric, or look at the data quality upstream.
The pure-Python code here is a learning artifact. The mental model it gives you is a production asset.
I run AYI NEDJIMI Consultants, a cybersecurity consulting firm. We publish free security hardening checklists — PDF and Excel.
Top comments (0)