---
title: An Approximation That Never Decides
date: '2026-08-24'
description: 'I wrote it down as a choice: the fast distance formula or the exact one. That was wrong. A formula that only shortlists can be as wrong as its speed allows, and choosing between three of them turned out to hinge on a number that runs opposite to throughput.'
author: Pablo Fernandez
tags:
  - Vector search
  - Quantization
  - Numerical precision
  - Rust
  - Benchmarking
source: 'https://www.pablofr.com/blog/an-approximation-that-never-decides'
---
I wrote down a trade that does not exist. The build path in [OpenVector](/blog/vector-search-on-object-storage) spends most of its arithmetic on distances between rows and centroids, computed one pair at a time, and the usual way to make a loop like that fast is to stop doing it one pair at a time and hand the whole batch to a matrix multiply. The matrix multiply gets there through an algebraic identity that loses precision on exactly the pairs that matter, the close ones. So I called it fast or exact, pick one, and moved on.

The framing was the error. The fast formula does not have to decide anything. It only has to narrow the field, and once it is only narrowing, its error stops being a correctness question and becomes a recall number you can count.

That fix is the easy half. The interesting half is what happens next: once the cheap pass is only a filter, you get to choose which cheap function fills the slot. There turned out to be three candidates, two of them already written and measured in this repo, and the one I had been arguing for is the only one nobody here has ever run.

## Where the arithmetic actually goes

A bulk build in OpenVector clusters a 100,000-row sample, then assigns every row in the corpus to one of the resulting leaves. The assignment walks a two-hop tree: rank the group centroids, then scan the seeds inside the best groups. At 8.7 million rows that comes out at roughly 100 groups of roughly 100 seeds each, because the group count is the square root of the seed count.

```
      one row
      │
      ▼
   ┌──────────────────────────────────────┐
   │  rank every group centroid           │   100 distances
   └──────────────────────────────────────┘
      │
      │  keep the nearest two groups
      ▼
   ┌──────────────────────────────────────┐
   │  scan every seed inside those groups │   200 distances
   └──────────────────────────────────────┘
      │
      │  near-tie at the top? redo with four groups
      ▼
   ┌──────────────────────────────────────┐
   │  the leaf this row belongs to        │
   └──────────────────────────────────────┘
```

*One row's path through the two-hop tree, and the two places it pays for distances.*

That is about 300 full distance computations per row, times 8.7 million rows, in 128 dimensions. The loop is written the obvious way, one call per pair:

```rust:ann/spfresh.rs
for &(_, g) in groups {
    for &si in &self.groups[g] {
        let d = distance(vector, &self.seeds[si], self.distance_metric);
        if d < best {
            second = best;
            best = d;
            best_seed = si;
        } else if d < second {
            second = d;
        }
    }
}
```

Each call streams 128 floats of the row and 128 floats of the seed through the arithmetic units, produces one number, and throws both operands away. The next call reloads the row. Two vectors in, one scalar out, no reuse of anything. The seed scan is two thirds of that work, and about a tenth of rows pay for it twice when the top two seeds come back near-tied.

## The filter, not the answer

Here is the move I missed. Nothing requires the cheap distance to produce the answer. It only has to produce a shortlist that contains the answer.

```
      200 candidate seeds
      │
      │   one cheap pass over all of them
      ▼
   ┌──────────────────────────────────────┐
   │  a ranked list, small but real error │
   └──────────────────────────────────────┘
      │
      │   keep the best 8, where 1 is needed
      ▼
   ┌──────────────────────────────────────┐
   │  8 exact distances, today's kernel   │
   └──────────────────────────────────────┘
      │
      ▼
      the winner, chosen on exact numbers only
```

*The cheap pass narrows 200 candidates to 8; the exact kernel decides among those 8.*

The 99% of the work, comparing a row against hundreds of candidates, goes down the cheap path. The 1% that settles anything goes down the precise one. It is a metal detector for a garden you did not want to dig up, and then a shovel where it beeps.

This is the standard shape, not an invention. [FAISS](https://github.com/facebookresearch/faiss) documents it plainly in its [implementation notes](https://github.com/facebookresearch/faiss/wiki/Implementation-notes): if you need k equals 10, query 100 and rerank the top-100 with exact distance computations, which is what their `IndexRefine` object does. The two-stage encoding is how production indexes are built.

## The candidate I was arguing for, and its failure mode

The matrix multiply route works like this. Stack the rows into an N by d matrix and the seeds into a d by C matrix, and one [GEMM](https://www.netlib.org/blas/) gives you every row-to-seed dot product at once. A blocked kernel loads a tile of each operand into registers and computes a whole tile of the output from it, so every value it loads takes part in many multiply-adds instead of one. The speedup is not fewer operations, it is the same operations with the memory traffic amortised.

But a GEMM gives dot products, and dot products are not distances. To get from one to the other you expand the square, so that `||x - c||^2` becomes `||x||^2 - 2*dot(x, c) + ||c||^2`. The norms are computed once per vector and once per centroid, so next to the N times C dot products they are free. FAISS ships both implementations and switches between them at `nq * d < 128000`, a threshold their notes say was found experimentally.

The expansion is exact in real arithmetic and lossy in floating point, and the loss is worst where you need it most.

```
      a row and a centroid
      │
      ├──▶  sum over i of (x[i] - c[i])^2
      │     every intermediate is the size of the answer
      │
      └──▶  ||x||^2  -  2*dot(x, c)  +  ||c||^2
            250000   -  499997       +  250000  = 3
            three big numbers, one small answer: the digits cancel
```

*Two routes to the same number: one keeps every intermediate small, the other builds three large numbers whose leading digits cancel.*

Those magnitudes are illustrative, not measured, but the arithmetic under them is exact. An f32 carries a 24-bit mantissa, so near 500,000 the gap between representable numbers is about 0.03. A squared distance of 3 computed that way arrives with roughly one part in a hundred of slop, purely from the final rounding. A neighbour sitting at 3.02 is a coin flip.

The kernel in the tree today has no such problem, because it never builds a big intermediate:

```rust:domain/distance.rs
fn euclidean_squared_scalar(a: &[f32], b: &[f32]) -> f32 {
    a.iter()
        .zip(b)
        .map(|(x, y)| {
            let d = x - y;
            d * d
        })
        .sum()
}
```

Every term here is already the size of the answer, so nothing large is subtracted from anything large and there is nothing to [cancel](https://en.wikipedia.org/wiki/Loss_of_significance). FAISS says the same about their own pair in one line: solution 2 may be less stable numerically than solution 1 for vectors of very different magnitudes, with a link to [the issue where it bit someone](https://github.com/facebookresearch/faiss/issues/297).

There are two repairs for that, and both are real. Subtracting the global mean from every vector and centroid once leaves the distance untouched, because squared distance does not care about translation, while shrinking the norms that every error term scales with. Keeping the norms and the final three-term combination in f64 removes two more roundings for the price of three operations per pair instead of 128. What neither repair removes is the error the f32 dot product accumulates inside the GEMM itself, which is the largest of the three terms. So the fast path stays approximate after both, and needs the exact rerank behind it regardless.

None of that is a reason not to do it. It is a reason to check whether there is a cheap function with no new failure mode at all.

## Two candidates were already in the repo

There are. Both of them read the same one-bit-per-dimension codes, and they turn those bits into a number in two entirely different ways.

OpenVector stores, next to every posting, a code for each row. The bit is the sign of the row's residual against its centroid, after a rotation. The rotation is not decoration: sign bits only estimate an angle when they behave like random hyperplanes, and along the raw coordinate axes correlated dimensions all flip together. The module records what that costs, measured on a 128-dimension cluster with a rerank of the top 50:

| what the sign bits are taken of | recall@10 |
| --- | --- |
| the raw vector, no centroid | 0.008 |
| the residual, no rotation | 0.448 |
| the rotated residual | the path that shipped |

The rotation itself is a random sign flip followed by a [Walsh-Hadamard transform](https://en.wikipedia.org/wiki/Fast_Walsh%E2%80%93Hadamard_transform), the standard cheap stand-in for a random orthogonal matrix. It runs in `O(d log d)` and preserves the Euclidean norm, which matters because the pruning bound is built from residual norms and a transform that moved them would quietly unsound it.

At 128 dimensions a row's vector is 512 bytes and its sign bits are 16, and the stored section works out at 29 bytes a row once a residual norm and an id hash ride along with them. Either way a candidate costs an order of magnitude fewer bytes to look at. What differs between the two candidates is what happens after those bytes are loaded.

**The table lookup, which is what a query runs.** The estimator on the default path keeps the query in full precision and puts only the stored side in bits. It precomputes, once per cluster, the partial sum for every one of the 256 sign patterns a byte of code can hold, so scoring a row is one table load and one add per byte instead of one multiply-add per dimension. At GIST's 960 dimensions that is 120 lookups rather than 960 adds. The estimate arrives inside a two-sided Cauchy-Schwarz bound, which is what makes pruning on it sound rather than merely plausible.

**The popcount, which nothing calls.** Binarise both sides, XOR, count the set bits. The kernels are public in the same file as the exact one:

```rust:domain/distance.rs
pub fn popcount_words(words: &[u64]) -> u64 {
    kernels::popcount_words(words)
}

pub fn bitwise_and_popcount(lhs: &[u64], rhs: &[u64]) -> u64 {
    kernels::bitwise_and_popcount(lhs, rhs)
}

pub fn bitwise_xor_popcount(lhs: &[u64], rhs: &[u64]) -> u64 {
    kernels::bitwise_xor_popcount(lhs, rhs)
}
```

Underneath they dispatch to an AVX-512 `VPOPCNTDQ` path that does the bitwise op and the popcount over eight `u64` at a time. Each has a benchmark and a unit test. What none of them has is a caller: the module that used to run this estimator in production is gone, and the symmetric version survives only as a measurement baseline under `tests/common/`.

So the shelf holds three filters, not one, and only one of them would need arithmetic written from scratch.

## The fastest one is the wrong one

Throughput first, because that is the number that nearly decided it. Both scan kernels are benchmarked on the same rows:

| filter | how it scores a candidate | measured here |
| --- | --- | --- |
| GEMM expansion | one matrix multiply, then the three-term combine | not measured, estimated 5 to 10x |
| codes table lookup | one load and one add per byte of code | 14 to 49 Mrows/s |
| XOR and popcount | one instruction per 64 dimensions | 88 to 276 Mrows/s |

The production counter agrees with the middle row: `ann:scan:rows` measured 114,838 rows in 3.03 ms, which is 37.9 Mrows/s, sitting inside the codes kernel's benchmarked range.

Seven times faster, on 21% fewer bytes per row. On throughput alone the popcount wins and it is not close, and that reading nearly concluded that three rounds of NEON work had gone into the wrong kernel. It is also why the repo grew a test that does not measure throughput at all.

```
   codes table lookup      29.0 B/row     14-49 Mrows/s     recall@10  0.9970
   XOR and popcount        24.0 B/row    88-276 Mrows/s     recall@10  0.7805
```

*Same posting, same 200 SIFT queries, same rerank width: throughput and recall point in opposite directions.*

Same cluster, same 200 SIFT queries, only the estimator differs, and the whole curve says the same thing:

| shortlist width | codes recall@10 | popcount recall@10 |
| --- | --- | --- |
| 10 | 0.4960 | 0.2485 |
| 20 | 0.6955 | 0.3540 |
| 50 | 0.8910 | 0.5165 |
| 100 | 0.9695 | 0.6535 |
| 200 | 0.9970 | 0.7805 |

For 21% more bytes the table lookup returns 0.997 where the popcount returns 0.781, and it is ahead at every width. The reason is the asymmetry: keeping the query in full precision and coding only the stored side preserves signal that binarising both sides throws away, for no saving in what has to be read.

That inverts how a filter should be graded. Throughput is what makes a filter attractive. What makes it safe is recall at a given shortlist width, because the filter's only job is to keep the true winner inside the shortlist. A kernel seven times faster that needs a much wider shortlist has spent its advantage on exact reranks, and exact reranks are the full-precision distances the filter existed to avoid. How much wider, the table does not say: at width 200 the popcount is still at 0.781 and the point where it catches 0.9695 lies outside what was measured.

Grade on recall at width, then check throughput. Doing it the other way round is what the repo caught itself doing.

One caveat that the asymmetry brings with it, and it is not small. On a posting scan the full-precision side is one query and the coded side is thousands of rows, so the partial-sums table is built once and amortised over all of them. In assignment the full-precision side is the row being placed, 8.7 million of them, and the coded side is about 200 seeds. Build the table per row and most of the saving goes into building it. Put the bits on the row instead and the table is built once per seed for the whole build, but then the row has to be coded against a centroid before it has been assigned to one, which is the thing being decided. The popcount has no such problem, because a symmetric kernel does not care which side it came from. Which side wears the bits is an open question, and it is the one place where the throughput ordering could still win.

## What the field does, which is not GEMM

I checked two other engines before committing, and both point the same way.

[Chroma](https://github.com/chroma-core/chroma) keeps its distance code in `rust/distance/src`, as four hand-written SIMD files: `distance_avx.rs`, `distance_avx512.rs`, `distance_neon.rs`, `distance_sse.rs`. Each one carries a header saying it is copied from [Qdrant](https://github.com/qdrant/qdrant) under Apache 2.0. They are pair-at-a-time distances, vectorised by hand, one file per instruction set. There is no GEMM and no BLAS dependency anywhere in the Rust tree. Their index is [SPANN](https://arxiv.org/abs/2111.08566), the same family as turbopuffer's and this one's, and the `spann_posting_list_*` values live in their blockstore.

[HelixDB](https://github.com/HelixDB/helix-db) is more interesting for what it has reserved than for what it runs. Its versioned on-disk encoding allocates stable identity numbers for `BinaryQuantizedCosineHalfV1`, `BinaryQuantizedSquaredEuclideanV1` and `BinaryQuantizedManhattanV1`, plus a `Hamming` metric and a `BinaryQuantizedV1` codec. Every one of them is marked reserved, and the codec currently returns `UnsupportedVectorCodec`. So this is not a shipped feature to point at, it is a format designer spending scarce numbers in a stable encoding on binary quantization before writing it, which is its own kind of vote.

Between them: nobody I looked at reaches for a matrix multiply, and the one making a long-term format bet is betting on bits.

## The comparison, and the decision

| | GEMM expansion | codes table lookup | XOR and popcount |
| --- | --- | --- | --- |
| new arithmetic to write | yes, with a cancellation failure mode | none | none |
| measured in this repo | no | yes, and running in production | kernels yes, no caller |
| throughput | estimated 5 to 10x | 14 to 49 Mrows/s | 88 to 276 Mrows/s |
| recall@10 at width 200 | not measured | 0.9970 | 0.7805 |
| how the winner is chosen | exact rerank | exact rerank | exact rerank |
| worst case if the filter is wrong | a mis-assigned row, so lost recall | a mis-assigned row, so lost recall | a mis-assigned row, so lost recall |

The bottom two rows are identical on purpose. All three are filters, all three hand the decision to the same exact kernel, and all three fail the same single way. Everything that separates them is above the line.

The GEMM goes last, and not because it is slow. It is the only one of the three that would need arithmetic written from scratch, its failure mode is precisely the one that costs recall, and it is the only one with no number next to it in this repo. Two candidates that have already been measured beat one that has to be built and then measured.

Between the other two, the popcount's seven times is real and so is its recall gap, and the gap is the one that matters for a filter. So the filter goes in with the table lookup that already runs on the query path, and the popcount stays as the fallback for the case where the shortlist can be widened cheaply enough to pay for it. That is a measurement, not a preference, and the next section is how it gets measured.

## The counter that makes it safe

There is exactly one failure mode left: the true nearest seed falls outside the cheap top-K, so the exact pass never sees it and cannot rescue it. Nothing else can be wrong, because nothing else is allowed to decide.

That failure is directly countable. Run the cheap ranking and the exact recompute side by side and count how often the exact numbers change the winner relative to the cheap order. If K is 8 and the answer moves zero times in a million rows, the margin is measured rather than argued. If it moves, raise K and measure again, and the cost of raising K is a handful of exact distances on a path that was doing 200.

The structure for that is already in the tree, because the same idea runs one level up. The code ranks the group centroids, takes the nearest two, and widens to four only when the row is ambiguous:

```rust:ann/spfresh.rs {5}
let two = ranked.len().min(2);
let (seed, best, second) = scan(&ranked[..two]);
// Near-tie between the two best seeds = a boundary row; re-check with
// twice the groups. ~10% of rows on real corpora, 2x their cost only.
if second.is_finite() && best.is_finite() && second <= best * 1.21 {
    let four = ranked.len().min(4);
    if four > two {
        return scan(&ranked[..four]).0;
    }
}
seed
```

A cheap ordering, a widened beam when the top is close, a decision taken after the widening. That is the filter pattern already, built for a different reason. What changes is only how the initial ordering is produced.

## Honest limits

Five of them, in the order they would bite.

- **The query-path numbers do not transfer.** That whole width table is top-10 out of a thousand-row posting. The build's assignment is top-1 out of about 200 candidates. Same kernels, different regime, and the curve has to be measured again where it will run. Quoting the query number as if it predicted the build number is exactly the substitution this project keeps catching itself making, and a top-1 decision may well forgive a looser estimator that a top-10 one does not.
- **The phase is small.** By my own attribution, `assign` is 32.9 seconds of a 1096-second build at 8.7 million rows, 3.0%. A twentyfold win there buys 31 seconds. The counting pass ahead of it is 56.9%, and 24.1% of the build is still unattributed, so even that ordering is provisional.
- **I have just been wrong about this exact thing.** The attempt before this one was an optimisation that came back 11% slower than the code it replaced, because I assumed where the bottleneck was instead of measuring it.
- **There is no stable recall baseline yet.** Changing how rows are assigned changes which cluster they land in, which changes recall. Without a baseline that holds still, the measurement after the change cannot separate a real regression from run-to-run noise, and a number that says nothing is worse than no number, because it looks like evidence.
- **Which side wears the bits is undecided.** The asymmetry that makes the table pay on a posting scan points the other way in assignment. That is a design question, not a tuning one, it is the first thing to settle after the baseline, and it is the one that could hand the slot back to the popcount.

So the order is: baseline first, then the codes filter, then the counter. The GEMM stays on the shelf, and if the codes filter ever stops being enough it is still there.

## The part worth keeping

Two things, and the second one is the one I keep relearning.

Precision is not a property of a formula, it is a property of where the formula sits. A formula that decides has to be right. A formula that only shortlists has to be right about the shortlist, which is a far weaker requirement and one you can put a counter on.

And a filter is not graded on how fast it runs. It is graded on how often the winner survives it at a given width, which is a different number, sometimes an inverted one, and the only one of the two that can be turned into a counter.

## References

- [FAISS implementation notes: the BLAS decomposition, its threshold, and IndexRefine](https://github.com/facebookresearch/faiss/wiki/Implementation-notes)
- [FAISS issue 297: where the expansion's instability was discussed](https://github.com/facebookresearch/faiss/issues/297)
- [Chroma's hand-written SIMD distance kernels](https://github.com/chroma-core/chroma/tree/main/rust/distance/src)
- [Qdrant, which those kernels are copied from](https://github.com/qdrant/qdrant)
- [HelixDB, whose vector encoding reserves binary-quantized score identities](https://github.com/HelixDB/helix-db)
- [SPANN: highly-efficient billion-scale approximate nearest neighbor search](https://arxiv.org/abs/2111.08566)
- [RaBitQ, the bound behind the one-bit codes](https://arxiv.org/abs/2405.12497)
- [Loss of significance (catastrophic cancellation)](https://en.wikipedia.org/wiki/Loss_of_significance)
- [The SIFT and GIST corpora these builds are measured against](http://corpus-texmex.irisa.fr/)
- [When the Vectors Live in a Bucket, the index this build feeds](/blog/vector-search-on-object-storage)