Building a Vector Search Engine from Scratch: HNSW Graphs, Cosine Distance, and Quantization

In modern artificial intelligence and machine learning engineering, almost all unstructured data—text passages, customer support transcripts, code repositories, audio spectrograms, and product catalog images—is represented as high-dimensional numerical vectors known as embeddings. Modern language models (such as OpenAI's text-embedding-3, Cohere Embed, or open-source BAAI/bge models) transform complex semantic concepts into coordinate arrays spanning 768, 1536, or 3072 floating-point dimensions.

Yet finding the nearest neighbors to a query vector in sub-millisecond latencies across millions of records presents a profound algorithmic challenge. A naive linear scan (calculating cosine similarity against every stored record) has time complexity $O(N \cdot D)$. At ten million records, a single nearest-neighbor query requires billions of floating-point operations, saturating memory bandwidth and driving latencies past several seconds.

To overcome this, industrial vector databases (such as Qdrant, Milvus, and pgvector) rely on Approximate Nearest Neighbor (ANN) algorithms. The undisputed state of the art is HNSW (Hierarchical Navigable Small World) graphs paired with Scalar Quantization (SQ8). In this systems guide, we will build a production-grade vector search engine from scratch in C# and .NET 8, implementing SIMD distance kernels, multi-layer skip-graph traversal, dynamic heuristic pruning, and 75% memory compression.

1. The Metric Space: Distance Functions and SIMD Hardware Acceleration

Before constructing graph topologies, we must establish how distance is computed in high-dimensional vector spaces. The three standard metrics are:

1. Euclidean Distance ($L_2$): Measures straight-line geometric distance:

   $$d_{L_2}(u, v) = \sqrt{\sum_{i=1}^{D} (u_i - v_i)^2}$$

2. Dot Product (Inner Product): Measures coordinate alignment and magnitude:

   $$\langle u, v \rangle = \sum_{i=1}^{D} u_i \cdot v_i$$

3. Cosine Distance: Measures the angular difference between vectors, independent of magnitude:

   $$d_{\text{cosine}}(u, v) = 1 - \frac{\sum_{i=1}^{D} u_i v_i}{\sqrt{\sum_{i=1}^{D} u_i^2} \sqrt{\sum_{i=1}^{D} v_i^2}}$$

The Normalization Optimization: If vectors are normalized to unit length ($\|u\|_2 = 1.0$) upon ingestion, cosine distance simplifies to: $d_{\text{cosine}}(u, v) = 1.0 - \langle u, v \rangle$. This eliminates expensive square root operations at query time and reduces the calculation to a pure dot product.

SIMD Hardware Vectorization

Because an inner product of two 1536-dimensional vectors involves 1,536 multiplications and 1,535 additions, scalar CPU loops bottleneck memory bandwidth. Using AVX2 / AVX-512 or ARM NEON, we process 8 to 16 floating-point numbers in a single clock cycle:

using System.Numerics;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;

public static unsafe float CosineDistanceSimd(float* a, float* b, int length)
{
    if (Avx2.IsSupported)
    {
        Vector256 sumVector = Vector256.Zero;
        int i = 0;

        // Process 8 float32 values per hardware instruction
        for (; i <= length - 8; i += 8)
        {
            Vector256 va = Avx.LoadVector256(a + i);
            Vector256 vb = Avx.LoadVector256(b + i);
            sumVector = Avx.Add(sumVector, Avx.Multiply(va, vb));
        }

        // Horizontal sum across vector lanes
        float dotProduct = Vector256.Sum(sumVector);

        // Process remaining tail elements
        for (; i < length; i++)
        {
            dotProduct += a[i] * b[i];
        }

        return 1.0f - dotProduct;
    }

    // Fallback scalar implementation
    float dot = 0f;
    for (int i = 0; i < length; i++) dot += a[i] * b[i];
    return 1.0f - dot;
}

2. Why Traditional Trees Fail: The Curse of Dimensionality

In low-dimensional spaces ($D \le 3$), spatial partitioning data structures such as k-d trees, R-trees, or quadtrees provide clean $O(\log N)$ nearest-neighbor search by recursively bisecting space with hyperplanes.

However, when dimensionality explodes into hundreds or thousands of dimensions ($D \ge 128$), the geometry of metric space changes radically:

• Volume Concentration: In high dimensions, almost all the volume of a hypersphere is concentrated in a razor-thin shell near its surface.
• Distance Equidistance: The distance between the closest pair of points and the farthest pair of points converges to nearly the same value.
• Hyperplane Degradation: Any hyperplane query in a k-d tree must visit nearly every leaf partition, degrading the search back to a linear scan ($O(N)$) while incurring immense pointer indirection overhead.

To scale search beyond $D = 128$, we must abandon geometric tree partitioning and transition to Proximity Graphs.

3. The Architecture of HNSW Graphs

Pioneered by Yury Malkov and Dmitry Yashunin, the Hierarchical Navigable Small World (HNSW) algorithm adapts the concept of 1D Skip-Lists into multi-dimensional metric graphs.

The Skip-Graph Concept

An HNSW index consists of multiple layers of directed proximity graphs:

• Top Layers (Layer $L_{\max} \dots 1$): Contain an exponentially decreasing subset of vectors. The edges in these upper layers bridge massive distances across the embedding space, acting as an express highway network.
• Base Layer (Layer 0): Contains 100% of the ingested vectors. The graph is tightly clustered, where every node is connected to its $M$ nearest topological neighbors.

During insertion, a node is assigned a maximum layer $l$ drawn from an exponential decay distribution governed by parameter $m_L = 1 / \ln(M)$:

$$l = \lfloor -\ln(\text{uniform}(0, 1)) \cdot m_L \rfloor$$

The Greedy Search Walk

Querying an HNSW graph is a process of greedy routing descent:

1. The search begins at the fixed Entry Point (EP) node on the highest layer ($L_{\max}$).
2. On each layer, the algorithm inspects all neighbors of the current node and greedily hops to whichever neighbor is closest to the query vector. It repeats this until no neighbor is closer than the current node (reaching a local minimum for that layer).
3. Once a local minimum is reached on Layer $k$, the algorithm drops down to Layer $k-1$ at the same coordinate and resumes the greedy walk.
4. Upon reaching Layer 0, the algorithm switches from a simple greedy walk to a bounded priority-queue search of size efSearch, returning the top-$k$ nearest neighbors with over 98% recall accuracy.

4. Memory Compression: Scalar Quantization (SQ8)

A single 1536-dimensional vector stored as 32-bit IEEE 754 floating-point numbers consumes:

$$1536 \text{ dimensions} \times 4 \text{ bytes} = 6,144 \text{ bytes / vector}$$

A dataset of 10,000,000 vectors requires 61.4 Gigabytes of RAM purely for raw embedding vectors, plus another 15 GB for HNSW graph edge adjacency lists. This requires expensive cloud instances with high memory provisions.

How Scalar Quantization (SQ8) Compresses Vectors

Scalar Quantization (SQ8) maps continuous 32-bit floating-point coordinates into discrete 8-bit unsigned integers ($0 \dots 255$):

$$q_i = \text{round}\left( \frac{x_i - \min_i}{\max_i - \min_i} \times 255 \right)$$

Each dimension now consumes exactly 1 byte instead of 4 bytes. A 1536-dimensional vector shrinks from 6,144 bytes to 1,536 bytes—an immediate 75% reduction in memory footprint!

Asymmetric Distance Computation (ADC)

To preserve precision during queries, high-performance engines use Asymmetric Distance Computation: the stored database vectors remain compressed as 8-bit integers, but the incoming search query vector is kept in full 32-bit floating-point precision. Distance calculations decode the quantized bytes on-the-fly inside CPU registers using precalculated scaling factors, retaining over 98.4% top-10 recall.

5. Building the Complete Vector Search Engine in C#

Below is a production-grade, multi-threaded C# implementation of an HNSW vector index supporting dynamic insertion, greedy multi-layer search, and SIMD distance calculations:

using System;
using System.Collections.Concurrent;
using System.Collections.Generic;
using System.Runtime.CompilerServices;

public sealed class HnswVectorIndex
{
    private readonly int _dimensions;
    private readonly int _m;             // Max neighbors per node on layer > 0
    private readonly int _m0;            // Max neighbors on layer 0 (typically 2 * M)
    private readonly int _efConstruction;// Beam search width during indexing
    private readonly double _mL;         // Normalization factor for level generation

    private readonly ConcurrentDictionary _vectors = new();
    private readonly ConcurrentDictionary[]> _graph = new();
    private int _entryPointId = -1;
    private int _maxCurrentLayer = -1;
    private readonly object _entryPointLock = new();

    public HnswVectorIndex(int dimensions, int m = 16, int efConstruction = 64)
    {
        _dimensions = dimensions;
        _m = m;
        _m0 = 2 * m;
        _efConstruction = efConstruction;
        _mL = 1.0 / Math.Log(m);
    }

    public void Insert(int id, float[] vector)
    {
        if (vector.Length != _dimensions)
            throw new ArgumentException($"Vector dimension must be {_dimensions}");

        _vectors[id] = vector;
        int nodeLevel = GenerateRandomLevel();
        var nodeLayers = new List[nodeLevel + 1];
        for (int l = 0; l (l == 0 ? _m0 : _m);
        }
        _graph[id] = nodeLayers;

        lock (_entryPointLock)
        {
            if (_entryPointId == -1)
            {
                _entryPointId = id;
                _maxCurrentLayer = nodeLevel;
                return;
            }
        }

        int currObj = _entryPointId;
        int maxLayer = _maxCurrentLayer;

        // Phase 1: Fast greedy descent from top layer down to nodeLevel + 1
        for (int l = maxLayer; l > nodeLevel; l--)
        {
            currObj = SearchLayerGreedy(vector, currObj, l);
        }

        // Phase 2: Insert into layers from min(maxLayer, nodeLevel) down to 0
        int startLayer = Math.Min(maxLayer, nodeLevel);
        for (int l = startLayer; l >= 0; l--)
        {
            int maxM = (l == 0) ? _m0 : _m;
            var candidates = SearchLayerEf(vector, currObj, _efConstruction, l);
            var selectedNeighbors = SelectNeighborsHeuristic(vector, candidates, maxM);

            // Connect bidirectionally
            foreach (var neighbor in selectedNeighbors)
            {
                _graph[id][l].Add(neighbor);
                lock (_graph[neighbor][l])
                {
                    _graph[neighbor][l].Add(id);
                    if (_graph[neighbor][l].Count > maxM)
                    {
                        PruneConnections(neighbor, l, maxM);
                    }
                }
            }

            if (candidates.Count > 0) currObj = candidates[0].Id;
        }

        lock (_entryPointLock)
        {
            if (nodeLevel > _maxCurrentLayer)
            {
                _maxCurrentLayer = nodeLevel;
                _entryPointId = id;
            }
        }
    }

    public List SearchKnn(float[] query, int k, int efSearch = 64)
    {
        if (_entryPointId == -1) return new List();

        int currObj = _entryPointId;
        for (int l = _maxCurrentLayer; l > 0; l--)
        {
            currObj = SearchLayerGreedy(query, currObj, l);
        }

        var results = SearchLayerEf(query, currObj, Math.Max(efSearch, k), 0);
        var finalNeighbors = new List();
        for (int i = 0; i  SearchLayerEf(float[] query, int startNode, int ef, int layer)
    {
        var visited = new HashSet { startNode };
        var candidates = new PriorityQueue();
        var results = new PriorityQueue();

        float startDist = ComputeDistance(query, _vectors[startNode]);
        candidates.Enqueue(startNode, startDist);
        results.Enqueue(startNode, -startDist); // Max-heap behavior for furthest result

        while (candidates.Count > 0)
        {
            candidates.TryDequeue(out int current, out float currentDist);
            results.TryPeek(out int furthestId, out float furthestDistNeg);
            float furthestDist = -furthestDistNeg;

            if (currentDist > furthestDist) break;

            var neighbors = _graph[current][layer];
            for (int i = 0; i  ef) results.Dequeue();
                    }
                }
            }
        }

        var sorted = new List();
        while (results.Count > 0)
        {
            results.TryDequeue(out int id, out float distNeg);
            sorted.Add((id, -distNeg));
        }
        sorted.Reverse();
        return sorted;
    }

    private List SelectNeighborsHeuristic(float[] baseVector, List candidates, int m)
    {
        var result = new List(m);
        for (int i = 0; i 

6. Production Benchmarks Across 1,000,000 Vectors

When evaluated against 1,000,000 OpenAI 1536-dimensional embeddings on an AMD EPYC 9654 server, our custom HNSW engine achieved the following operational metrics: