Building a Vector Search Engine from Scratch with HNSW in Python

Every vector database tutorial opens with pip install pinecone and skips the part where you understand why vector search works. That's fine for shipping fast, but it costs you later — when recall drops and you don't know where to look, or when a tuning parameter turns your latency from 10ms to 800ms. HNSW (Hierarchical Navigable Small World) is the algorithm behind almost every ANN (Approximate Nearest Neighbor) library in production today: FAISS, Qdrant, Weaviate, and hnswlib all use it internally. This post builds a working HNSW index in pure Python — not production-grade, but real enough to give you a concrete mental model. The naive approach to nearest neighbor search is brute force: compare your query vector to every stored vector, return the k closest. That's O(n·d) per query — for 1M vectors at 1536 dimensions, expect around 6 seconds per lookup. HNSW solves this with a layered graph. Vectors live across multiple layers: Higher layers: sparse "highway" connections between distant nodes Lower layers: dense "local" connections between nearby nodes Layer 0: the complete graph with the richest connectivity At query time: Enter at the top layer via a fixed entry point Greedily walk toward nodes closer to the query Drop to the next layer, repeat At layer 0, collect the k nearest neighbors Expected complexity drops from O(n) to O(log n). The trade-off: more memory, and approximate results — you may occasionally miss the single true nearest neighbor, but you'll be well within a small margin. Here's a minimal HNSW implementation that captures the core structure: import numpy as np import heapq import math import random from dataclasses import dataclass, field from typing import Optional @dataclass class HNSWNode: vector: np.ndarray neighbors: dict = field(default_factory=dict) # layer -> list of node ids class HNSWIndex: def __init__(self, dim: int, M: int = 16, ef_construction: int = 200, max_layers: int = 6): self.dim = dim self.M = M # max neighbors per layer self.M0 = M * 2 # max neighbors at layer 0 self.ef = ef_construction self.max_layers = max_layers self.nodes: list[HNSWNode] = [] self.entry_point: Optional[int] = None self.top_layer = 0 def _distance(self, a: np.ndarray, b: np.ndarray) -> float: diff = a - b return float(np.dot(diff, diff)) # squared L2 def _random_level(self) -> int: level = 0 while random.random() < (1.0 / self.M) and level < self.max_layers - 1: level += 1 return level def _search_layer(self, query: np.ndarray, entry_id: int, ef: int, layer: int) -> list[int]: visited = {entry_id} entry_dist = self._distance(query, self.nodes[entry_id].vector) candidates = [(entry_dist, entry_id)] dynamic_list = [(entry_dist, entry_id)] while candidates: c_dist, c_id = heapq.heappop(candidates) worst_dist = max(d for d, _ in dynamic_list) if c_dist > worst_dist: break for nb_id in self.nodes[c_id].neighbors.get(layer, []): if nb_id not in visited: visited.add(nb_id) nb_dist = self._distance(query, self.nodes[nb_id].vector) if nb_dist < worst_dist or len(dynamic_list) < ef: heapq.heappush(candidates, (nb_dist, nb_id)) dynamic_list.append((nb_dist, nb_id)) if len(dynamic_list) > ef: dynamic_list.remove(max(dynamic_list)) return [nid for _, nid in sorted(dynamic_list)[:ef]] def add(self, vector: np.ndarray) -> int: node_id = len(self.nodes) node = HNSWNode(vector=vector) self.nodes.append(node) level = self._random_level() if self.entry_point is None: self.entry_point = node_id self.top_layer = level for l in range(level + 1): node.neighbors[l] = [] return node_id ep = self.entry_point for l in range(self.top_layer, level, -1): candidates = self._search_layer(vector, ep, ef=1, layer=l) ep = candidates[0] for l in range(min(level, self.top_layer) + 1): neighbors = self._search_layer(vector, ep, ef=self.ef, layer=l) M = self.M0 if l == 0 else self.M neighbors = sorted( neighbors, key=lambda nid: self._distance(vector, self.nodes[nid].vector) )[:M] node.neighbors[l] = neighbors for nb_id in neighbors: nb_node = self.nodes[nb_id] if l not in nb_node.neighbors: nb_node.neighbors[l] = [] nb_node.neighbors[l].append(node_id) if len(nb_node.neighbors[l]) > M: nb_node.neighbors[l] = sorted( nb_node.neighbors[l], key=lambda nid: self._distance(nb_node.vector, self.nodes[nid].vector) )[:M] ep = neighbors[0] if neighbors else ep if level > self.top_layer: self.top_layer = level self.entry_point = node_id return node_id def search(self, query: np.ndarray, k: int = 10, ef: int = 50) -> list[tuple[float, int]]: if self.entry_point is None: return [] ep = self.entry_point for l in range(self.top_layer, 0, -1): candidates = self._search_layer(query, ep, ef=1, layer=l) ep = candidates[0] candidates = self._search_layer(query, ep, ef=max(ef, k), layer=0) results = sorted(candidates, key=lambda nid: self._distance(query, self.nodes[nid].vector)) return [(self._distance(query, self.nodes[nid].vector), nid) for nid in results[:k]] import numpy as np from hnsw import HNSWIndex # Build an index with 128-dimensional vectors index = HNSWIndex(dim=128, M=16, ef_construction=200) # Insert 10k random vectors (documents, embeddings, etc.) np.random.seed(42) n = 10_000 vectors = np.random.randn(n, 128).astype(np.float32) print("Building index...") for i, vec in enumerate(vectors): index.add(vec) if (i + 1) % 1000 == 0: print(f" {i + 1}/{n} vectors inserted") # Query with a random vector query = np.random.randn(128).astype(np.float32) results = index.search(query, k=5, ef=50) print("\nTop 5 nearest neighbors (HNSW):") for dist, node_id in results: print(f" node_id={node_id}, sq_distance={dist:.4f}") # Verify against brute force all_dists = np.sum((vectors - query) ** 2, axis=1) brute_top5 = sorted(enumerate(all_dists), key=lambda x: x[1])[:5] print("\nBrute-force top 5:") for nid, dist in brute_top5: print(f" node_id={nid}, sq_distance={dist:.4f}") On a 10k-vector dataset, you'll typically see 4 or 5 out of 5 results overlap between HNSW and brute force. That's the "approximate" in ANN: for most real workloads, occasionally missing one true neighbor is an acceptable trade-off for a 50x–100x query speedup. M — max connections per node per layer. Higher M means better recall and a denser graph, at the cost of more memory and slower inserts. Start with 16. Go up to 32 if your recall consistently falls below 90%. ef_construction — candidate pool size during index build. Larger values produce a better-connected graph but slow down insertion. 200 is a solid default; lower it to 100 if insertion throughput matters more than recall. ef (search-time) — candidate pool size at query time. This is the most important lever for production tuning because you can adjust it without rebuilding the index. Setting ef = k gives the fastest queries; ef = 500 approaches brute-force recall. Plotting ef vs. recall typically reveals a knee around ef = 50–100 for most datasets — that's your sweet spot. Understanding these parameters matters the moment your AI-powered search starts returning wrong results. Before blaming your embedding model, check ef and whether your index was built with sufficient ef_construction. It's the same diagnostic discipline that applies across security tooling — knowing the internals saves time when something breaks. The security hardening checklists we publish follow the same philosophy: understand the mechanism, then apply the control. This implementation deliberately omits several things you'd need in production: Deletion support: HNSW graphs don't support clean node removal without a full rebuild. Libraries like Qdrant implement soft-deletion with tombstones. Persistence: serialize node vectors and neighbor lists to disk (numpy's np.save for vectors, JSON or msgpack for the neighbor maps). Thread safety: concurrent inserts will corrupt the neighbor lists. Use a lock or batch inserts single-threaded. SIMD distance computation: the Python distance function is the bottleneck. hnswlib's C++ core runs distance in vectorized instructions — easily 100x faster. For any production workload, use hnswlib directly (pip install hnswlib) or a vector database that wraps it. The API is almost identical to what's shown above, and you get persistence, filtering, and proper concurrency for free. HNSW is not magic — it's a graph with controlled connectivity. The layered structure lets queries skip large sections of the search space, and the ef parameter gives you a runtime dial between speed and accuracy. Building a stripped-down version reveals what every production vector database is doing under the hood. That understanding pays off when you're debugging recall issues at 3am — and it's faster than re-reading library documentation when you already know what the parameters mean. I run AYI NEDJIMI Consultants, a cybersecurity consulting firm. We publish free security hardening checklists — PDF and Excel.