How I write SIMD kernels_

Astral's vector index finds neighbors by comparing a query against stored vectors. Each candidate visit ends in a dot product. One path reads ordinary f32 values. Another reads E5M2, an 8-bit floating-point format used for the first approximate comparisons. The best candidates are then rescored against their original f32 vectors.

I tried the same change in both dot-product kernels: keep several partial sums so the CPU has independent multiply-add chains to run. It helped E5M2 and did nothing for f32. Four accumulators cut E5M2 insert p99 by 1.5%. Six and eight were slower. Later, scoring two E5M2 records together cut top-10 p50 by 2.86% and p99 by 4.34% over 100,000 vectors.

The generated instructions explained the split. E5M2 must unpack and convert each byte before multiplication, so its loop has conversion work to overlap with those independent sums. The f32 loop loads values already in the needed format. More accumulators gave it little useful work to hide.

The result has to keep the answers right

The raw E5M2 measurements in this article hold 92.03% recall@10. I used that number to catch any quality change while comparing kernels. Astral's 99% requirement applies to the high-recall storage configurations that rerank candidates in f32, not to this raw E5M2 lane.

I read recall before latency, then p99 before p50. A faster kernel does not stay if it changes the neighbors returned, and an average can hide the slowest one-in-a-hundred operations.

I measured every result in this post on an Intel Core Ultra 7 265 with 20 cores, GCC 13.3.0, and Ubuntu 24.04. The numbers are specific to that machine and toolchain. The same changes need fresh measurements on different hardware.

The dot kernel

The read kernel operates on E5M2.

One byte holds 1 sign bit, 5 exponent bits, and 2 mantissa bits. The format is compact and coarse, so the byte scores route the graph search and the survivors are reranked against their original f32 vectors. The conversion post covers the bit layout and recall trade.

The kernel problem starts after storage.

x86 has no instruction that loads an E5M2 row and scores it directly. The loop has to turn bytes into f32 lanes while it computes the dot product, and any kernel change still has to preserve enough approximate ordering to meet the recall requirement.

The dot kernel converts while it scores. It loads 8 bytes, zero-extends them into 32-bit lanes, and shifts each lane left by 8. That places the bits where a half-float expects them. It repacks the lanes as half floats, converts them to f32, then runs the FMA loop.

The converter handles eight values at a time. Finding this sequence took five attempts in one day. The conversion post covers them.

Each iteration converts bytes into f32 lanes and accumulates their score with FMA instructions. Converting the next bytes does not depend on the previous FMA result, so the CPU can overlap the two kinds of work. The accumulator count determines how much independent work is available.

Accumulators

Set the conversion aside for a moment, because the accumulator question already exists in a plain f32 dot product. The naive loop keeps one running sum:

// one running sum: every fma depends on the one before it
float dot(const float* a, const float* b, int n) {
    float s = 0.0f;
    for (int i = 0; i < n; ++i)
        s = std::fma(a[i], b[i], s);
    return s;
}

Every FMA reads the previous sum, so the loop has one dependency chain. An FMA may take several cycles to finish even though the core can start independent FMAs every cycle. With an illustrative four-cycle latency and two starts per cycle, one chain uses only one of the eight starts available during those four cycles.

The fix is independent partial sums:

// four independent chains, folded once at the end (n a multiple of 4)
float dot4(const float* a, const float* b, int n) {
    float s0 = 0, s1 = 0, s2 = 0, s3 = 0;
    for (int i = 0; i < n; i += 4) {
        s0 = std::fma(a[i+0], b[i+0], s0);
        s1 = std::fma(a[i+1], b[i+1], s1);
        s2 = std::fma(a[i+2], b[i+2], s2);
        s3 = std::fma(a[i+3], b[i+3], s3);
    }
    return (s0 + s1) + (s2 + s3);
}

Four partial sums create four independent chains. Their latency overlaps, and the loop combines them only at the end. The pairwise reduction (s0 + s1) + (s2 + s3) also cuts the dependency tree from three levels to two.

More chains are not always better. They consume registers and compete with loads and conversion instructions. The useful count has to be measured.

The compiler will not make this change under the current flags. Splitting one float sum into four changes addition order, and floating-point addition is not associative. I write the partial sums explicitly so the accumulator count remains a measured source-level decision.

The scalar E5M2 kernel already used four partial sums before the SIMD version existed. The first AVX2 kernel carried two vector accumulators and processed 16 values per iteration. That width had not been measured.

Here is the two-accumulator kernel and the x86-64 loop emitted by Clang 18. It zero-extends the bytes, shifts them into half-float positions, packs them to 16 bits, converts them, and accumulates with vfmadd231ps.

// widen 8 E5M2 bytes to 8 f32 through the F16C unit
static inline __m256 widen_e5m2_x8(const uint8_t* p) {
    __m256i w = _mm256_slli_epi32(_mm256_cvtepu8_epi32(       // byte -> own 32-bit lane
                    _mm_loadl_epi64((const __m128i*)p)), 8);  // << 8 into half position
    __m128i h = _mm_packus_epi32(_mm256_castsi256_si128(w),
                                 _mm256_extracti128_si256(w, 1)); // pack lanes -> 8 halves
    return _mm256_cvtph_ps(h);                                   // 8 half -> 8 f32 (F16C)
}
// the E5M2 dot kernel, two vector accumulators (n a multiple of 16)
float dot_e5m2_2acc(const uint8_t* a, const float* b, int n) {
    __m256 s0 = _mm256_setzero_ps(), s1 = _mm256_setzero_ps();
    for (int i = 0; i < n; i += 16) {
        s0 = _mm256_fmadd_ps(widen_e5m2_x8(a + i),     _mm256_loadu_ps(b + i),     s0);
        s1 = _mm256_fmadd_ps(widen_e5m2_x8(a + i + 8), _mm256_loadu_ps(b + i + 8), s1);
    }
    return hsum256(_mm256_add_ps(s0, s1));   // hsum256 = extract-add + two hadds
}

; clang-18 -O2 -mavx2 -mfma -mf16c (x86-64), dot_e5m2_2acc inner loop
.LBB0_5:
    vpmovzxbd    ymm2, qword ptr [rdi + rcx]             ; 8 E5M2 bytes -> 8x u32
    vpslld       ymm2, ymm2, 8                           ; byte << 8, into half position
    vextracti128 xmm3, ymm2, 1                           ; high 4 lanes
    vpackusdw    xmm2, xmm2, xmm3                        ; pack 8x u32 -> 8x u16 (halves)
    vcvtph2ps    ymm2, xmm2                              ; 8 half -> 8 f32   (F16C)
    vfmadd231ps  ymm1, ymm2, ymmword ptr [rsi + 4*rcx]   ; acc0 += widened * b
    vpmovzxbd    ymm2, qword ptr [rdi + rcx + 8]
    vpslld       ymm2, ymm2, 8
    vextracti128 xmm3, ymm2, 1
    vpackusdw    xmm2, xmm2, xmm3
    vcvtph2ps    ymm2, xmm2
    vfmadd231ps  ymm0, ymm2, ymmword ptr [rsi + 4*rcx + 32] ; acc1 += widened * b
    add          rcx, 16
    cmp          rcx, rax
    jb           .LBB0_5

Two accumulators means two vfmadd231ps instructions into ymm0 and ymm1, each with its own widen stream feeding it. Going to four duplicates the complete widen-plus-FMA block twice more and advances by 32 components. The four-accumulator excerpt shows only the accumulator destinations. The full loop has four copies of the earlier conversion sequence:

; accumulator destinations in the 4-accumulator loop:
    vfmadd231ps ymm0, ymm4, ymmword ptr [rsi + 4*rcx]      ; acc0
    vfmadd231ps ymm1, ymm4, ymmword ptr [rsi + 4*rcx + 32] ; acc1
    vfmadd231ps ymm3, ymm4, ymmword ptr [rsi + 4*rcx + 64] ; acc2
    vfmadd231ps ymm2, ymm4, ymmword ptr [rsi + 4*rcx + 96] ; acc3

A standalone Clang 18 reproduction kept all four sums and two widening temporaries in registers. The loop had no vector spill or stack reference.

Two accumulators to four

E5M2 has more independent work than f32. Its zero-extend, shift, pack, and conversion instructions do not wait on the FMA chains. Four accumulators give the scheduler four FMA chains alongside that conversion work. The arithmetic count stays the same, but more instructions can run while an earlier FMA is pending.

I widened both E5M2 dot kernels.

They moved from two vector accumulators to four, process 32 values per iteration, and reduce four sums at the end. One kernel scores bytes against an f32 query, while the other scores bytes against an E5M2 query.

I measured the change while inserting 100,000 vectors. Each insert searches the graph for neighbors and calls these kernels many times. P95 fell from 1.926 to 1.899 milliseconds, and p99 fell from 2.079 to 2.049 milliseconds. P50 and recall did not move.

The widening removed about 1.5% from insert p99 and nothing from the median or recall, so it stayed. That run did not say where additional accumulators would stop helping.

Scoring two records together

I tried six accumulators, then eight. Both were slower. The next kernel scored two records together.

The accepted candidate scores two adjacent records in one fixed kernel. Each record keeps two accumulators, so the kernel still has four live sums. Conversion and FMA work remain proportional to the number of records, but the outer record loop and its result-control cadence run once per pair. An odd final record falls back to the existing one-record kernel.

These are medians across three interleaved process runs, each containing 2,048 warm queries, on the Intel Core Ultra 7 265 with three runtime workers plus the caller restricted to CPUs 0-3. Both top-10 recall runs measured 92.03% over the same 64 deterministic queries.

Median cycles fell 2.92%, and branch instructions fell 16.96%. Total instructions fell 0.13%, branch misses fell 0.88%, and cache misses rose 0.46%. The pair kernel still performs twice the conversion and FMA work. It saves record-loop and result-control work.

Where the widening stops

Each accumulator occupies a vector register for the whole loop, and conversion needs scratch registers too. AVX2 has sixteen vector registers. Wider strides can also create more tail handling for dimensions that do not divide evenly.

Four accumulators helped the E5M2 dot kernels. A wider f32 loop did not, so I reverted it. The neighboring L2 kernels still use two accumulators. Each loop has a different mix of conversion, arithmetic, and loads.

I do not have a recorded wider test for the L2 kernels, so I cannot say whether two is best. The next change there needs a new measurement.

The read version now overlaps conversion with four independent sums. The write-side quantizers are a separate problem and need their own aligned measurements.

Validating the instructions

I generally check the generated instructions for every kernel in GCC, Clang, and MSVC. This is part of the normal work.

On Windows, that check showed the problem immediately. MSVC had compiled the slower fallback and emitted no vcvtph2ps. The F16C intrinsics were available, but the feature macro used by the build guard was missing. I fixed the guard with an internal build macro. Runtime CPUID and XCR0 checks still decide whether F16C is safe to use on the machine.

I also check that the compiler kept the conversion and fused multiply-add operations I expected. The conversion post covers the compiler guard and runtime checks in more detail.

The paired-record kernel remained in the E5M2 search path at 92.03% recall, with top-10 p50 down 2.86% and p99 down 4.34%.