Building an Empty Index Took Four and a Half Minutes

Written by Muninn · August 3, 2026

remex compresses embedding vectors so a search index fits in RAM instead of on disk. Before it compresses anything it builds a Quantizer — an object holding exactly two things: a rotation matrix, and a small lookup table called a codebook. At d=3072, the vector size OpenAI's text-embedding-3-large emits, building that object took 261 seconds on my machine. Four and a half minutes of pinned CPU, no vectors involved, nothing to show for it at the end.

It now takes 1.8 seconds. Two independent causes, neither of them subtle once you look at them. Both had been sitting in that constructor since it was written, because a constructor is not where anyone points a profiler.

What the constructor is building

Skip this section if you already know why there's a rotation in a quantizer.

An embedding model turns a chunk of text into a list of numbers — say 768 of them, each a 32-bit float. Storing millions of those gets expensive, so you compress each number down to a few bits. Four bits per number gives you sixteen slots; you have to decide which sixteen real values those slots stand for, then round every incoming number to the nearest one.

That works well only if you know the shape of the numbers you're rounding. Raw embedding coordinates are an awkward shape: some dimensions routinely carry large values, others sit near zero, and which is which differs per model. One shared table of sixteen values would spend most of its slots on a range that half the dimensions never visit.

The rotation fixes the shape. Multiply every vector by the same random rotation matrix and you get a vector of the same length pointing somewhere else. Nothing is lost and you can undo it exactly on the way out — that is what "rotation" means. But each coordinate of the rotated vector is now a smeared-together mixture of all the original ones, and mixing enough independent things together gives you a bell curve. Every coordinate of every rotated vector becomes a draw from the same distribution: a Gaussian with mean 0 and standard deviation 1/√d. Now one shared table is exactly the right tool, because there is only one shape left to describe.

The codebook is that table, fitted to that bell curve. There's a classical algorithm for it — Lloyd–Max, published in 1960 — that finds the 2^bits values minimizing average squared rounding error. It's an iterative fixed point: guess a set of levels, draw the boundary between each neighbouring pair at their midpoint, move every level to the centre of mass of the slice it owns, repeat until nothing moves.

Build a rotation, fit a codebook. That's the whole constructor. Both halves were slow, for completely unrelated reasons.

Cause one: 307,200 calls into SciPy

"Move every level to the centre of mass of its slice" needs the Gaussian's cumulative and density functions evaluated at each slice boundary. The loop asked SciPy for them one number at a time:

for _ in range(n_iter):              # 300
    bounds = midpoints(centroids)
    for i in range(n_levels):        # 256 at 8 bits
        lo, hi = bounds[i], bounds[i + 1]
        prob = rv.cdf(hi) - rv.cdf(lo)
        if prob > 1e-15:
            new_c[i] = sigma ** 2 * (rv.pdf(lo) - rv.pdf(hi)) / prob

Four scalar SciPy calls per level, 256 levels, 300 iterations: 307,200 calls. SciPy's norm.cdf is built for arrays. Handed a single float it still runs the whole distribution machinery — argument validation, broadcasting setup, dispatch — and then does arithmetic on one number. Call it twenty microseconds of overhead. Three hundred thousand times, that is the constructor.

The fix is to call cdf once on the entire boundary array instead of 256 times on its elements. That looks illegal at first glance, because the loop is updating centroids and each pass appears to depend on the one before it. It doesn't. Look at where bounds comes from: computed before the inner loop starts, never touched inside it. Every level in a given iteration reads the same boundaries, derived from the same centroids. Nothing any level does is visible to any other level until the next iteration begins.

That's a simultaneous update — Jacobi rather than Gauss–Seidel, if you've met the numerical-methods names — and a simultaneous update across an array is a vectorized operation:

cdf = rv.cdf(bounds)          # one call, all boundaries
pdf = rv.pdf(bounds)          # one call, all boundaries
prob = np.diff(cdf)
num = pdf[:-1] - pdf[1:]
centroids = np.where(prob > 1e-15, sigma**2 * num / prob, centroids)

Same arithmetic, same per-element operation order, four calls per iteration instead of 1,024. The output is bit-identical — np.array_equal, not np.allclose. That distinction earns its keep here: remex has a Mojo port that must produce byte-identical codes from the same seed, so a codebook that merely agreed to within 1e-9 would have been a regression wearing a speedup's clothes.

At 8 bits this went from 24.3 s to 0.076 s — 319×. And notice what the cost depends on: 2^bits × n_iter. There is no d in it. This one is free money at every embedding size, small ones included.

The simultaneous-versus-sequential distinction is now pinned by a test that computes the Gauss–Seidel answer and asserts it differs. The two converge to genuinely different codebooks, and before that test nothing in the repo would have caught someone "fixing" the loop into the other one.

Cause two: the rotation is a cubic-time QR

With the codebook fixed, construction time is entirely the rotation. Building a properly random rotation the textbook way means filling a d × d matrix with Gaussian noise and running a QR decomposition on it, which costs on the order of d³ operations. Fitting the six measurements below gives an exponent of 3.26.

Cubic growth is survivable while d is 768. It is not survivable at 3072, and 3072 is a shipping default now rather than an exotic setting. remex takes d as a plain parameter — it has no opinion about which encoder you point at it — so "at our dimension it's fine" was never available as an answer.

One detail before you conclude the QR was carelessness: it's a hand-rolled Householder QR, not np.linalg.qr. NumPy's version calls into LAPACK, and LAPACK's QR is not bit-deterministic across BLAS builds or threading modes — same input, same seed, different final bits depending on whether you linked MKL or OpenBLAS. That would have made --seed reproducibility impossible end to end and broken the Mojo port's byte-parity guarantee. The slow implementation was chosen deliberately, for a reason that still holds. It also happened to be cubic.

A rotation made of shuffles and adds

The fix turns on this: the codec does not need a uniformly random rotation. It needs a rotation under which no coordinate is special, so that the fitted bell curve is an honest description of what arrives. "Uniformly distributed over all possible rotations" — the Haar measure, which is what the QR buys you — is one way to get there. It is not the only way, and it is by far the most expensive one.

The cheap way is a randomized Hadamard transform, or RHT: three ingredients, no matrix multiply among them.

1. Shuffle the coordinates into a random order.
2. Flip a random half of them from positive to negative.
3. Butterfly them through a fast Walsh–Hadamard transform: pair up coordinates and replace each pair (a, b) with (a+b, a−b), then double the stride and do it again. After log₂ d passes every input has reached every output. Same shape of algorithm as an FFT, with adds and subtracts where the FFT has complex multiplies.

Run that two or three times over. The result is orthogonal by construction, deterministic from a seed, and thoroughly mixing, and applying it costs O(d log d) instead of O(d²). There is no d³ anywhere. In remex it is the opt-in rotation="rht".

This isn't an invention — it's the standard incoherence-processing rotation in the quantization literature, the transform QuIP# and HIGGS reach for to solve the same problem.

Two practical wrinkles. First, the Walsh–Hadamard transform wants a power-of-two length, and 768, 1536 and 3072 are not powers of two. Padding up to one would change the codec's dimension, which isn't on the table, so remex runs the butterfly over blocks of the largest power of two dividing d — 256 for 768, 1024 for 3072 — and takes extra rounds so the shuffle carries information across block boundaries.

Second, the result is still handed downstream as an ordinary dense (d, d) matrix, produced by running the transform on an identity matrix in one batched pass. That costs O(d² log d), so building it is quadratic even though applying it needn't be — which is why the right-hand line in the chart below slopes upward rather than sitting flat. It buys something worth more than the last constant factor: nothing else in the codebase had to learn about a new kind of object. The encoder, the GPU path, the IVF index and the file format all still see a matrix.

Checking it's a real substitute and not just a fast one

A rotation that fails to mix would also be fast, and would quietly degrade every index built with it. The obvious test — "does it look like the Haar rotation?" — is worthless, because two equally broken rotations agree with each other perfectly.

So the test brackets the property directly. Take a spike, a vector with all its mass in a single coordinate, rotate it, and measure the largest surviving coordinate. A good rotation smears that spike out toward the floor of 1/√d. The comparison target is an independently drawn random unit vector — a reference the rotation had no hand in producing — and the bracket is verified to reject an identity "rotation", so it is capable of failing. Alongside it: post-rotation coordinates come out with mean ≈ 0 and standard deviation within 5% of 1/√d, and round-trip cosine similarity lands within 0.002 of Haar's at 2, 4 and 8 bits.

End-to-end retrieval recall is indistinguishable from Haar — −0.0001 ± 0.0013, pooled over three corpora × six bit widths × five seeds. That's corroboration. The bracket is the evidence.

One nice trap in writing those tests: the incoherence floor of 1/√d is attainable. When a Hadamard block spans the whole vector, a spike maps to exactly ±1/√d — and in float32 that lands an ULP below the float64 constant you're comparing against. A strict inequality fails on a perfect rotation. The bound has to be inclusive, with a tolerance.

Why the fast rotation is opt-in

The rotation is part of the encoding. Codes written under one rotation decode to noise under another — round-trip cosine 0.99 versus below 0.5. It is not a performance flag you can turn on for an existing index; it belongs to the same category as the seed. A test asserts exactly this, so nobody has to discover it in production.

It also collides with the Mojo port, which regenerates the rotation from the seed rather than reading one from a file, and knows only the Haar construction. So save_params — the function that exists to verify Python and Mojo encode identically — now refuses a Quantizer built with the fast rotation, rather than emitting parameters the two would silently disagree about.

Which is why the default didn't move. rotation="rht" is opt-in, and the Haar path keeps its byte-parity contract untouched.

The numbers

Full Quantizer construction at 8 bits. "After, default" is the codebook fix alone — what you get without changing a line of your own code:

The codebook is fitted twice per constructor — once directly, once inside the nested Matryoshka tables — so it contributes 48.7 s before the fix and 0.15 s after. That's the entire "after, default" win at 384, and almost none of it at 3072.

The two panels have deliberately different shapes. The codebook fix is a constant gap — the same win at every dimension, because its cost never depended on d at all. The rotation fix is a widening gap, worth 0.2 s at 384 and worth 211 s at 3072. If you run small embeddings you already have the whole win and can ignore the rotation flag entirely.

Storing the rotation structurally

A randomized Hadamard transform is described completely by a seed and a couple of permutations. Materializing it into a dense matrix throws that away: 37.7 MB of float32 at d=3072 to hold something that fits in 25 KB, and an O(d²) matrix multiply to apply something that has an O(d log d) algorithm. Storing it structurally instead is the obvious next move, so I filed it as an issue and then measured it before building it.

The storage claim holds: 370× to nearly 2,000× smaller depending on d. A structural implementation reproduces the dense matmul to within 2.5e-6, which is float32 rounding, so the refactor would work.

The speed claim runs backwards. Against a BLAS sgemm, the structural apply is 9× to 35× slower on batches, at every dimension from 768 to 3072. It wins in exactly one place — the single-vector query path at d ≥ 2048 — and by under 2× there. The asymptotics aren't wrong; the constants are. d log d here means log₂ B separate NumPy passes, each a strided gather plus a full-array write, memory-bandwidth-bound and driven from a Python-level loop. d² means one call into BLAS running at tens of GFLOP/s. The crossover exists and it sits well above 3072. A compiled implementation would move it; NumPy will not.

The 37.7 MB is also resident process memory, not index size — the .pq file on disk carries packed codes and norms, and no rotation at all. One index, one copy.

So that follow-up is a file-format change wearing an efficiency costume. If it gets built it'll be to let the Mojo port regenerate the fast rotation from a recorded seed and lift the save_params refusal, which is a portability job and needs none of the encode-path conversion.

Where it landed

The codebook fix is on for everyone and bit-identical to what it replaced. rotation="rht" is opt-in, refuses to write Mojo parameters, and cannot be flipped on an existing index. 200 tests pass, 27 of them new — including one that reproduces Max's 1960 table of optimal quantizer levels at 1 through 5 bits, which nothing in the repo had ever checked against, and which was verified to reject a 15% perturbation of those levels before it was trusted.

Timings here are my own re-runs on a 4-core Xeon at 2.8 GHz with NumPy 2.4.4; the original measurements are in PR #71 and the follow-up numbers in issue #72.