# Building a Vector Search Engine from Scratch with HNSW in Python

> Source: <https://dev.to/ayinedjimi-consultants/building-a-vector-search-engine-from-scratch-with-hnsw-in-python-18hj>
> Published: 2026-09-09 10:03:05+00:00

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:

At query time:

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:

``` python
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]]
python
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](https://ayinedjimi-consultants.fr/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:

`np.save` for vectors, JSON or msgpack for the neighbor maps).
For any production workload, use [hnswlib](https://github.com/nmslib/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](https://ayinedjimi-consultants.fr), a cybersecurity consulting firm. We publish [free security hardening checklists](https://ayinedjimi-consultants.fr/checklists) — PDF and Excel.*