Building a Custom Vector Database - Part 2: Contiguous Storage and AVX2 Vectorization

In our previous article, we established a baseline vector search engine using SQLite’s Virtual Table API. Our profiling results revealed that 77.64% of the search execution time was spent inside a single scalar function: calculate_l2_distance.

If we want to speed up our search, we have to look closely at how modern CPUs handle math. In standard code, a CPU processes numbers one at a time (scalar execution). But modern CPUs have a hidden superpower: they can process multiple numbers simultaneously using a technique called SIMD (Single Instruction, Multiple Data).

This article is designed as a practical, beginner-friendly introduction to hardware-level parallelization. We will:

1. Flatten our database storage to ensure our data sits side-by-side in memory.
2. Deconstruct the CPU’s vector registers to understand how they work.
3. Rewrite our distance calculation loop using Intel AVX2 instructions to calculate 8 dimensions in a single step.
4. Benchmark the performance and profile new hardware-level bottlenecks.

Part 1: Why Memory Layout Rules Hardware Performance

Before we can ask the CPU to perform parallel calculations, we must arrange our data correctly.

In our original database, we stored coordinates using a nested structure:

// Previous State: An array of dynamically allocated arrays
std::vector<std::vector<float>> vectors_;

In standard programming, this is highly intuitive. However, this structure forces the actual float coordinates to be scattered across different locations on the system heap. Accessing each vector requires pointer chasing—reading the memory address of the sub-array first, then jumping to the actual coordinates.

To utilize parallel hardware, we need our data to be completely sequential. We achieve this by flattening our layout into a single, one-dimensional array:

// Current State: All coordinates stored in one contiguous block of memory
std::vector<float> vectors_;

When we insert a new vector, we simply append its coordinates directly to the end of this single array:

void VectorEngine::insert(int64_t id, const std::vector<float>& vec) {
  if (vec.size() != dimensions_) {
    throw std::invalid_argument("Vector dimension mismatch.");
  }
  ids_.push_back(id);
  vectors_.insert(vectors_.end(), vec.begin(), vec.end());
}

Nested Storage (Pointer Chasing):
vectors_ -> [ Ptr 0 ] > Heap Block 0: [ x0, y0, z0 ] (Scattered)
            [ Ptr 1 ] > Heap Block 1: [ x1, y1, z1 ] (Scattered)

Flat Storage (Contiguous & Sequential):
vectors_ -> [ x0, y0, z0, x1, y1, z1, x2, y2, z2 ... ] (All side-by-side)

Now, instead of passing high-overhead C++ std::vector references to our distance calculation, we can pass raw, sequential pointers and step through them using simple memory offsets:

const float* query_ptr = query.data();
const float* dataset_ptr = vectors_.data();

for (size_t i = 0; i < ids_.size(); i++) {
  // Step directly to the memory address of vector 'i'
  float dist = calculate_l2_distance(query_ptr, dataset_ptr + i * dimensions_);
  results.push_back({ids_[i], dist});
}

Part 2: What is SIMD and AVX2?

To understand how we parallelize calculations, we need to look at CPU registers.

Typically, when you perform math in C++, the compiler uses scalar registers. For example, adding two floats together happens in a 32-bit register.

AVX2 (Advanced Vector Extensions 2) is an instruction set that unlocks larger, 256-bit registers called YMM registers.

      A SINGLE 256-BIT YMM REGISTER
+-----+-----+-----+-----+-----+-----+-----+-----+
| 32b | 32b | 32b | 32b | 32b | 32b | 32b | 32b |  <- 8 Slots (Lanes)
|Float|Float|Float|Float|Float|Float|Float|Float|
+-----+-----+-----+-----+-----+-----+-----+-----+

Because a single-precision floating-point number (float) is 32 bits wide, a single 256-bit register can hold exactly eight floats at the same time.

Using SIMD, we can load eight floats from Vector A into one register, eight floats from Vector B into another register, and ask the CPU to subtract or multiply them all in a single clock cycle.

Register A:  [ a0 ][ a1 ][ a2 ][ a3 ][ a4 ][ a5 ][ a6 ][ a7 ]
                 -     -     -     -     -     -     -     -    (Single SUB instruction)
Register B:  [ b0 ][ b1 ][ b2 ][ b3 ][ b4 ][ b5 ][ b6 ][ b7 ]
                 =     =     =     =     =     =     =     =
Result:      [ diff0 ][ diff1 ][ diff2 ][ diff3 ][ diff4 ][ diff5 ][ diff6 ][ diff7 ]

For further information on AVX2, check out the wikipedia page: Advanced Vector Extensions

Part 3: Implementing AVX2 in C++

To write SIMD instructions directly in C++, compilers provide intrinsics. Intrinsics are special compiler-specific functions that look like standard library calls but translate directly into single CPU instructions. Intrinsics are accessible by including <immintrin.h>.

Because AVX2 registers process exactly eight floats at a time, we must ensure our vector dimensions are a multiple of eight:

VectorEngine::VectorEngine(size_t dimensions) : dimensions_(dimensions) {
  if (dimensions % 8 != 0) {
    throw std::invalid_argument("Vector dimension must be a multiple of 8 for AVX2 processing.");
  }
}

Here is the updated, AVX2-accelerated $L_2$ distance calculation function:

#include <immintrin.h>

float VectorEngine::calculate_l2_distance(const float* a, const float* b) const {
  // 1. Create a 256-bit register initialized with 8 zeros
  __m256 sum_vec = _mm256_setzero_ps();

  // 2. Loop through dimensions in steps of 8
  for (size_t i = 0; i < dimensions_; i += 8) {
    // Load 8 floats from pointers a and b into vector registers
    __m256 va = _mm256_loadu_ps(a + i);
    __m256 vb = _mm256_loadu_ps(b + i);
    
    // Subtract va - vb across all 8 slots simultaneously
    __m256 diff = _mm256_sub_ps(va, vb);
    
    // Fused Multiply-Add: (diff * diff) + sum_vec
    sum_vec = _mm256_fmadd_ps(diff, diff, sum_vec);
  }

  // 3. Move vector data back to a standard C++ array
  alignas(32) float buffer[8];
  _mm256_store_ps(buffer, sum_vec);

  // 4. Sum the 8 final accumulated elements together
  float sum = 0.0f;
  for (int i = 0; i < 8; i += 1) {
    sum += buffer[i]; 
  }
  return sum;
}

Deconstructing the Code: Step-by-Step

Let’s translate Intel’s naming convention and explore exactly what each intrinsic does.

Naming Conventions Decoded

All AVX2 intrinsics follow a structured naming pattern:

• _mm256: Operates on 256-bit registers (YMM).
• _ps: Operates on “Packed Single-precision” values (i.e., multiple float values packed into one register).

__m256 sum_vec = _mm256_setzero_ps();

The type __m256 represents a physical 256-bit vector register. The function _mm256_setzero_ps clears this register, setting all 8 float slots (or “lanes”) to 0.0f. This register acts as our running accumulator.

_mm256_loadu_ps(a + i)

This instruction loads eight sequential float values from a memory address into a YMM register.

• The u in loadu stands for unaligned. In computer systems, loading data from memory addresses aligned to 32-byte physical boundaries is fastest. Since we cannot guarantee that every vector in our flat array starts on a clean 32-byte boundary, loadu safely loads the memory regardless of alignment.

_mm256_sub_ps(va, vb)

This instruction performs element-wise subtraction across all eight slots simultaneously: va[n] - vb[n].

_mm256_fmadd_ps(diff, diff, sum_vec)

This instruction executes a Fused Multiply-Add (FMA) operation. Normally, squaring the differences and adding them to the running total requires two distinct operations: a multiplication followed by an addition. FMA does both in a single execution step:

$\text{sum\_vec} = (\text{diff} \times \text{diff}) + \text{sum\_vec}$

This combines two arithmetic operations into one instruction pipeline, minimizing rounding errors and doubling the mathematical throughput of our loop.

Gathering the Results: Horizontal Addition

After processing all dimensions, our squared differences are calculated and stored inside the eight separate slots of the sum_vec register:

sum_vec: [ Slot 0 ][ Slot 1 ][ Slot 2 ][ Slot 3 ][ Slot 4 ][ Slot 5 ][ Slot 6 ][ Slot 7 ]

We need to add these eight independent numbers together to obtain our single scalar distance value. To do this, we copy the data out of the CPU’s register and back into standard system memory:

alignas(32) float buffer[8];
_mm256_store_ps(buffer, sum_vec);

We declare our target float buffer using alignas(32). This guarantees that our stack buffer is aligned to a physical 32-byte boundary, allowing us to use the highly efficient, aligned store instruction _mm256_store_ps (which lacks the u of loadu) to write the data back instantly.

Finally, we use a simple loop to add those eight numbers together into our final scalar float:

float sum = 0.0f;
for (int i = 0; i < 8; i += 1) {
  sum += buffer[i]; 
}
return sum;

Part 4: Benchmarking the Optimized Engine

By pinning the execution to a single physical core (taskset -c 2) and locking the frequency governor, we generated highly stable metrics across $10,000$ vectors:

=== VectorEngine (Direct Call - Search Space: 10,000 vectors) (Cycles) ===
  Min:    469690
  P50:    657160
  P95:    1681494
  Max:    2608174
  Average:786518

=== SQLite Virtual Table (Via Virtual Interface - Search Space: 10,000 vectors) (Cycles) ===
  Min:    495946
  P50:    668952
  P95:    1642746
  Max:    2579008
  Average:809567

Performance Analysis

Let’s compare the $P_{50}$ cycles of our Phase 2 implementation (AVX2 + contiguous memory) with our Phase 1 implementation (scalar + nested vector):

• Direct Call ($P_{50}$):
  • Phase 1 (Baseline): $1,364,802$ cycles
  • Phase 2 (Optimized): $657,160$ cycles ($2.07\times$ speedup / $51.8\%$ reduction in latency)
• SQLite Virtual Table ($P_{50}$):
  • Phase 1 (Baseline): $1,486,334$ cycles
  • Phase 2 (Optimized): $668,952$ cycles ($2.22\times$ speedup / $55.0\%$ reduction in latency)

Notably, the absolute SQLite virtualization overhead has shrunk from $121,532$ cycles to just $11,792$ cycles ($\approx 1.7\%$ relative overhead at $P_{50}$). The execution engine is now fast enough that memory subsystem interactions and system overhead have begun to dominate our execution profiles.

Part 5: Profiling Analysis and the New Bottleneck

While our computational optimizations succeeded, a perf report analysis highlights a new set of system bottlenecks that were previously masked by the slow scalar math.

+   11.51%     0.00%  benchmark_drive  [kernel.kallsyms]     [k] asm_exc_page_fault
+    8.23%     0.00%  benchmark_drive  [kernel.kallsyms]     [k] exc_page_fault
+    6.62%     0.34%  benchmark_drive  [kernel.kallsyms]     [k] do_anonymous_page
+   11.75%     4.37%  benchmark_drive  libc.so.6             [.] __memmove_avx_unaligned_erms
+    2.23%     0.23%  benchmark_drive  libstdc++.so.6.0.33   [.] operator new(unsigned long)

Analyzing these results reveals two major bottlenecks:

1. Page Fault and Memory Mapping Penalty: Over 11% of total execution cycles are spent inside the Linux kernel handling page faults and allocating anonymous memory pages (asm_exc_page_fault).
2. Memory Moving and Heap Allocation Overhead: Another 11.7% of cycles are spent in memmove (__memmove_avx_unaligned_erms), and $\approx 2.2\%$ are spent on heap allocation (operator new).

Diagnosing the Source: Heap Churn in VectorEngine::search

The primary cause of this overhead is the allocation and handling of results inside the search hot path:

std::vector<SearchResult> results;
results.reserve(ids_.size()); // Allocates 10,000 SearchResult structs on every single query

for (size_t i = 0; i < ids_.size(); ++i) {
    ...
    results.push_back({ids_[i], dist});
}

// Partial sort on the entire 10,000 element heap block
std::partial_sort(results.begin(), results.begin() + actual_k, results.end(), ...);

On every single query, our database:

1. Allocates space for $10,000$ elements on the heap.
2. Triggers page faults as we write to this unmapped heap memory.
3. Writes $10,000$ values.
4. Executes sorting operations over a large $10,000$-element footprint.
5. Deallocates the memory.

Even with reserved memory, allocating and clearing $10,000$ results on every single query causes significant heap churn and operating system overhead.

Phase 3: Memory Reuse and Algorithmic Sorting Optimization

To bypass these memory and page fault penalties, we must eliminate heap allocations inside our query paths. In the next part of this series, we will transition our sorting strategy:

1. Bounded Max-Heap (Priority Queue): Instead of storing all $10,000$ elements and then sorting them, we can maintain a Max-Heap of size $K$ (where $K = 10$).
2. Reduced Allocations: A heap of size $K$ only requires $K+1$ elements of storage. We can stack-allocate a small local buffer or reuse a pre-allocated thread-safe block, eliminating heap allocations, memmove operations, and page faults.
3. Algorithmic Scaling: Instead of pushing $10,000$ items and sorting them, we will only perform heap insertions when a vector’s distance is smaller than the current largest distance in our size-10 heap.

This shifts our algorithmic complexity and removes the OS memory allocator from our query execution path entirely.