Knuth's Root of All Evil Has a Break-Even Point

Don Knuth said premature optimization is the root of all evil: about 97% of the time, small efficiencies don't matter. In the same paragraph he said not to pass up the critical 3%. The hard part is telling which is which.

Finding every pair of near-identical items in a big collection is how you line up candidate pairs for duplicate review, or build a similarity graph, a web that links each item to the ones most like it. I built an exact similarity join, which finds every such pair and misses none, that ran 3.041x faster than FAISS, the widely used library for this kind of search, on a million real arXiv title embeddings, end to end, and returned exactly the same pairs. An embedding turns each title (these come from arXiv, a big public archive of research papers) into a long list of numbers, so titles that mean similar things get similar lists; end to end means the stopwatch ran for the whole job. Then I built a GPU version (a GPU is the graphics chip, built to do thousands of small sums at once) with a fast path that answers queries, single complete runs of the search, on those titles 2.5x faster, and taught it to skip that path when the setup, the prep work done once before the first query, won't pay for itself. On arXiv titles, the fast path comes out ahead from the ninth query.

I set out to beat tiled FAISS by at least 2x on the same Apple M2 Max, return every expected pair, and keep memory flat no matter how many results came out. I wrote that down before the deciding run. On a million 768-dimensional vectors (lists of 768 numbers each) with 500,000 planted duplicate pairs (copies put in on purpose), the MLX engine, my Mac version built on Apple's MLX toolkit, took 191.423 seconds to FAISS's 518.203. That's 2.707x, and all 18 confirmation runs came back exact. The arXiv run also counted the time to write every scored pair to disk.

On NVIDIA GPUs I added the fast path, which skips most of the arithmetic when most pairs can't match, and a selector that picks which way each job runs. Building the selector made setup cost part of the algorithm.

The engine has two ways to answer the same question. The direct way compares every pair of rows in full (a row is one item's list of numbers). The screening way first builds a shorter summary of every row, uses the summaries to rule pairs out cheaply, and compares only the rest in full, the way you'd skim a thumbnail before opening the full photo. The summaries cost time up front, and screening pays only when the pairs it rules out save enough on each run to earn that back over the runs you plan. So you tell the engine how many runs you plan. It times a small sample of the work and picks the way it estimates will finish the rest sooner. Told the job will run once, it skips the sample and takes the direct way.

A threshold join owes you every qualifying pair

A threshold join returns every unordered pair of rows whose cosine similarity, the dot product of the two rows after each is scaled to unit length, is at or above a threshold, the minimum score that counts as a match. (A dot product multiplies two lists number by number and adds up the results; scaling each list to length one, so the score measures only how closely two rows point the same way.) Each row can have zero qualifying neighbors or thousands, and the join has to be fine with either. That is the shape you need to build a similarity graph or to produce candidate pairs for duplicate review. A top-k search returns a fixed number of neighbors per row, the k nearest, so a row with no real neighbors still gets its k nearest even when they fall below the threshold, and a row with many gets cut off at k. Wrong for this job in both directions at once.

With 50,000 rows, the join must consider 1,249,975,000 unordered pairs, and the full FP32 score matrix would take 9.3 GiB. (FP32 is the standard way a computer stores a number with a decimal point, in 32 bits.) So the engine never builds it. It visits the pair space in tiles, square patches of the grid of all pairs, handled one patch at a time like mowing a lawn in strips, and streams each tile's qualifying pairs out as a batch of indices (row numbers) and scores. A 50,000-row join with 8,192-row tiles takes 28 initial tiles; the Notes list their shapes.

The CUDA engine has two routes through those tiles. (CUDA is NVIDIA's toolkit for programming its graphics chips.) The dense route does the direct thing and computes every similarity with a full FP32 cuBLAS matrix product (cuBLAS is NVIDIA's ready-made math library; the matrix product does all of a tile's multiply-and-add in one go). The screening route starts with dot products between lower-dimensional projections of the rows, the shorter summaries from earlier, uses those to bound each pair's full similarity, meaning it works out the highest that similarity could possibly be, and computes full dot products only for the pairs it cannot reject. On 50,000 arXiv title embeddings at threshold 0.9, on an RTX 4090, a high-end NVIDIA graphics card, with the data already on the GPU and the projection built, screening answered a query in about 26 ms (thousandths of a second), including copying the output to the host, the computer's main memory. Dense took about 64 ms. Both returned the same complete set of pairs. Note the clause "with the projection built." The rest of this essay is about that clause.

Every rounding step rounds the safe way, or none of them count

Screening rejects a pair when an upper bound on its similarity is below the threshold. A rejected pair never reaches the exact path, and nothing later can bring it back. So an upper bound that came out even slightly too low would discard a qualifying pair for good.

That requires saying exactly what is bounded. The prepared API, the set of calls a program makes to the engine once setup is done, takes host-normalized FP32 vectors, decides membership by their exact real dot product (the true answer, as if computed with no rounding at all) against an FP32 threshold, and returns FP32 approximations as scores.

The projection is scaled so it cannot increase a vector's length. For each pair, the screen combines the projected dot product with an upper bound on the contribution the projection missed, and together they bound the full dot product from above. In real arithmetic, that's the whole argument. In floating point, of course, every operation rounds, and one rounding in the wrong direction is exactly the too-low bound the screen can't afford. (Floating point is how computers do decimals: each result gets rounded to the nearest number the computer can store.) So the part of each row's squared length that the projection captured has to be rounded down, and the part it missed rounded up. This is the core of the CUDA bound kernel, screen_bounds.cu (a kernel is a small program that runs on the graphics chip):

double lower=fmax(0.,__dsub_rd(__dsqrt_rd(sum),error));double captured=__ddiv_rd(__dmul_rd(lower,lower),beta);double remaining=fmax(0.,__dsub_ru(norm_squared,captured));out[row]=bound_normal_up(__double2float_ru(__dsqrt_ru(remaining)));

sum is a downward bound on the projection's squared norm, error bounds the projection error, and beta bounds the projection operator's squared norm. Each intrinsic, a built-in chip instruction, names its rounding direction: _rd rounds toward negative infinity and _ru toward positive infinity. Every operation rounds in the direction that keeps the bound safe. The code never adds one fudge tolerance at the end and hopes. The full kernel also handles FP32 subnormal conversion and flushing (subnormals are the tiniest numbers near zero, which need special care), which these four lines don't show.

Accepting a pair is a separate procedure, and both routes share it. Screening passes on the pairs its bound could not reject. Dense passes on every pair whose FP32 score clears the threshold minus that score's error bound, rounded down. Each candidate's full FP32 score then gets an error interval, a range the true score is guaranteed to sit inside. If the interval's lower endpoint is at or above the threshold, the pair is emitted. If its upper endpoint is below the threshold, the pair is discarded. Otherwise the pair is re-evaluated with FP64 intervals (FP64 is the same format with twice the bits, so far finer rounding) under the same rule. If that still leaves it undecided, the engine stops approximating altogether: integer arithmetic computes the exact dot product of the FP32 operands and emits the pair when that product is at or above the threshold. Three tiers of arithmetic for one yes-or-no question, and the last tier doesn't round at all. The derivation and the operation-by-operation error ledger are in the project's NUMERICS.md. CPU checks, run on the computer's main processor with Python's Fraction, which does arithmetic on exact fractions, test the inequalities exactly, and independent CUDA fixtures exercise the device arithmetic.

Comparing two routes follows the same reasoning: they must return the same multiset of pairs (the same pairs, each the same number of times) under the same predicate. A score tolerance can't excuse a missing pair, and equal counts can hide a missing pair and an extra one that cancel perfectly.

Bound the memory you own, and split a tile before it says a word

On the Mac, emitting 49.5 million clustered pairs kept peak memory, the most in use at any one moment, close to a run that emitted none. On the GPU, a score tile is only the first allocation of a query, the first chunk of memory it borrows. After scoring come per-block counts, scan scratch (temporary working space), candidate indices and scores, classification flags, and the compacted output, and during compaction the old candidate arrays and the new compacted arrays are live at the same time. A budget that counted only the score tile could be exceeded by a dense result even when the score tile fit, which makes it a budget in name only. So the query planner bounds each allocation phase and takes the largest; Appendix A has the function and what each term counts. The default budget is 512 MiB, a cap on temporary graphics-chip memory, and it covers the temporary device arrays a query owns.

The executor, the part that actually runs the tiles, checks the score and count bound before scoring a tile, and the remaining phases after counting, before allocating candidate arrays. If a tile won't fit, it bisects the tile along its longer side and puts both halves back on the queue. The split always happens before that tile emits any output. That ordering is the whole trick. The halves are disjoint and nothing from the original tile has gone out, so no pair can be emitted twice, and the engine never keeps a global record of what it has sent.

In a separate diagnostic on an RTX A6000, another NVIDIA card, at revision d1a7f1e (a saved version of the code), measured owned peaks were 258.0 to 258.5 MiB, all within the computed bounds, and every completed query left zero live bytes in a dedicated allocation pool.

Sparse scoring pays only when screening actually rejects something

Where screening spends its time depends on the data. An earlier A6000 profile at d1a7f1e measured it, using Nsight Systems kernel-time sums inside the measured query ranges, excluding preparation and warmup. (A profile times each step to see where the time goes; Nsight Systems is NVIDIA's tool for it; warmup is the untimed practice runs.) On arXiv at 0.9 with screening forced, bound counting and scattering took 31.32 ms, the projected matrix product 12.01 ms, and the sparse full dot products, full-dimensional scores for just the surviving pairs, 3.29 ms. On an isotropic control at 0.8, where the bound rejects very little, the sparse full dot products alone took 5,171.23 ms. (A control is a made-up test set whose shape you know in advance; isotropic means random points spread evenly in every direction.) That is not a typo. I checked. On selective data the work left after rejection is small and the bound itself is the main cost. On unselective data the survivors are the main cost. Change the data and the bottleneck moves.

The engine already counts each tile's survivors to size the candidate arrays, so the same count also decides whether sparse scoring still makes sense for that tile. From _screen_gpu.py:

if use_dense(total, bounds, self.dense_fallback_ratio):    # No prefix arrays remain while the full score tile is allocated.    del scores, ends    stats["dense_fallback_tiles"] += 1    return self._dense_candidates(bounds, threshold, stats, stage)

The cutoff started at 25% of the tile's eligible pairs. In the first installed campaign (a campaign is a planned batch of benchmark runs: timed, repeatable tests), eight-cluster controls, data built as eight tight clumps, kept about 12.5% of their pairs as candidates. That is under 25%, so they stayed on the expensive sparse path, faithfully obeying a cutoff set too late. I changed the cutoff to 1% and froze it before the confirmation run. On the A6000 at d1a7f1e, that cut those controls' query medians, the middle value of the runs, from about 715 ms to about 167 ms at thresholds 0.8 and 0.9, and arXiv 0.8 from about 1,009 ms to about 191 ms. But dense was still faster on those cells (a cell is one dataset at one threshold): a tile that falls back arrives having already paid for the projected product and the count.

One measurement goes the other way. In the RTX 4090 run at d1a7f1e, clustered data at 0.95 took 54.28 ms with screening forced, 59.88 ms with dense, and 67.77 ms with the 1% fallback. The cutoff that rescued the clustered controls gave away this small win. The policy was frozen, so it stayed frozen, and that number went into the report beside the big ones.

Give the baseline equal engineering attention, and watch what it does with it

Before screening existed, the GPU comparison that mattered was my custom compaction kernel, which packs the matching pairs into a tight list, against fused CuPy masking on the same GPU, with both sharing the matrix product and the streaming. CuPy, a popular Python library for GPU math, was the baseline: the standard approach you have to beat for the new one to count. The CUDA stage had its own written gate, a pass-or-fail bar set in advance, separate from the Mac one: at least a 2x median speedup over an equally tuned standard GPU control, on real and synthetic workloads. At 2,048-row tiles the custom kernel was 1.20 to 1.38x faster. Then I tuned both paths the same way. Both chose 8,192-row tiles, and the advantage fell to 1.01 to 1.06x. I reran the comparison before I believed it. The lead had lived at a tile size neither path picks once it's tuned. The first screening prototype didn't beat tuned CuPy on 100,000 arXiv rows either. That gate failed, and that optimization campaign stopped right there. The CPU baseline told the same story from the other side: correcting its BLAS dispatch, which math routine it called, made it about 3.85x faster.

Later, profiling showed overhead in the bound and scan steps and in classifying every pair unconditionally. Two kernel fusions, each merging separate steps into one pass, roughly halved screening time on the real workloads. Against the strongest dense route tested, the fused screening prototype's prepared query was 2.88x faster on arXiv (51.4 ms against a 148.3 ms dense control) and 1.78x faster on COCO (photo embeddings). GloVe (word embeddings) had an apparent screening advantage right up until the dense control improved, and then the two were essentially tied. Those runs used 50,000 vectors at threshold 0.9 and five repetitions, with setup and calibration, the small timed sample from earlier, excluded, on the prototype benchmark code. That was enough to reopen screening for another round. The rule I took from both rounds: an improvement the baseline could also use goes into the baseline too.

Estimate the work the executor actually performs, all 28 tiles of it

The selector uses calibration measurements to choose dense or a screening rank, weighing a one-time cost, preparing the projection for the full dataset, against a per-query saving times the number of queries the caller plans. (The rank is how short the summaries are: rank 64 keeps 64 numbers per row.) Whether to calibrate is a separate decision. Once calibration has run, its cost is the same whichever route wins, so it can't make either route cheaper from there. It counts only when the whole job is compared with dense run without calibration. This sounds like accounting trivia. It is half of what went wrong next.

The first version of the selector chose dense in every cell of the confirmation campaign at d1a7f1e, on both GPUs, and missed all five measured screening wins. Its cost model had two faults. First, it timed one diagonal sample tile and compared that tile's saving, unscaled, with preparation costs for the full dataset, though a 50,000-row join has 28 initial tiles. And it charged calibration against screening, though calibration had already been paid by the time the choice was made. For arXiv at 0.9 on the RTX 4090, that model put screening's break-even at 413 queries, far past the 100 it was told to plan for. So: dense, everywhere. Well, of course.

The fixed selector, at revision 7d46ea5, times each of the four tile shapes through the complete query pipeline: candidate counting, fallback, subdivision, refinement and transfers to the host. It samples and shuffles 16,384 rows, enough for two full tiles. From that sample it constructs all four tile shapes at their actual execution sizes, including the 848-row strips at the edges, left over because 50,000 rows don't split evenly into 8,192-row tiles. It multiplies each shape's median time by the number of tiles of that shape in the full join. It times the projection's preparation on the sample and scales it to the full row count, because that preparation is still ahead when the choice is made. Then it chooses the route with the lowest future cost, from experimental.py:

future = {    name: expected_queries * value + (        0.0 if name == "dense" else preparation_seconds.get(name, 0.0)    )    for name, value in estimates.items()}

Calibration is already spent by then, so it stays out of the choice, though the reported total still includes it. A screening route also has to beat dense's estimated query time by at least 10% to be considered, and a tie goes to dense.

For one planned query, the selector (strategy "auto" in the code, "Auto" in the tables) skips calibration and chooses dense:

if strategy == "auto" and expected_queries > 1:    resident, selection = calibrate(

Otherwise it builds the dense resident data and goes dense. No sample, no calibration bill.

Before the hardware run, I added 18 holdout workloads to the 18 existing ones: thresholds 0.85, 0.925 and 0.975 instead of 0.8, 0.9 and 0.95, the next 50,000 rows of each real corpus, and a fresh synthetic seed. Holdouts are fresh tests the policy never saw while it was being built, so it can't have been tuned to them. I froze the policy before running the existing and holdout workloads, and no corpus name or threshold selects a route.

The campaign ran all 36 workloads on the RTX 4090 at 7d46ea5: arXiv (768 dimensions), COCO (512), GloVe (25), and isotropic, eight-cluster and dense-output synthetic controls. Each workload measured dense, forced screening, adaptive screening and the selector, with three first-use observations and nine prepared-query observations per route: 1,728 main timings, every output fully consumed. The campaign used expected_queries=100 to test the selector's repeated-query decisions.

The selector that had chosen dense everywhere now chose screening in 8 of 18 existing workloads, with prepared-query speedups of 1.45 to 3.17x over dense, and in 10 of 18 holdout workloads, with speedups of 1.20 to 3.17x. It chose dense in the rest. Across all 36 workloads, the selector's prepared-query time was never more than 2.5% slower than the fastest of the three explicit routes tested: dense, forced screening and adaptive screening at rank 64. Because the selector can also pick rank 128, it sometimes beat all three of the routes it was being graded against. For arXiv at 0.9 it chose rank-64 screening, and its six selection reports for that cell put the break-even at 9.0 to 11.1 queries counting calibration, and 1.8 to 2.3 queries counting only the projection's preparation. The old model had said 413.

Its query-time estimates had a median absolute error of 3.04% across the 36 cells, and 14 of 18 cells in each phase were within 7%. The six dense-output controls were underestimated by 33.54 to 36.07%, and the selector chose dense for all six. Its biggest misses all landed on cells where it picked dense anyway.

Charge setup to useful reuse, because setup always sends the bill

The campaign measured two timing scopes. Prepared-query time covers one query on an object that already exists, the engine already set up and sitting in memory, including copying batches to the host and concatenating them. First-use time starts from the host rows and covers normalization, device setup, calibration, preparation, the first complete query and closing the object, with every route warmed before measurement.

The two scopes rank the routes differently. At first-use scope, all 36 of the selector's medians were slower than dense. First use is where setup lands, and one fast query doesn't pay for it. These are three of the cells, in milliseconds:

The last column is the break-even for a caller who keeps running the same prepared join; Appendix C has the arithmetic for existing arXiv 0.9.

After the first query, auto is 288.07 ms behind dense. Each additional query claws back 38.23 ms, which puts auto ahead from the ninth query. The selector had estimated 9.0 to 11.1 before any of these queries ran. The measurement said nine. I enjoyed that one. The fast path loses the first query and wins the job.

The third scope is the whole job. The project's example program builds a similarity graph: it loads the rows, sets up, runs one query, and writes the edges, the links between similar titles, to a CSV file, a plain spreadsheet-style text file. I ran it on the existing 50,000 arXiv rows at 0.9 once per strategy, each in a fresh process after the campaign. The total starts before input loading and ends after the CSV is written. It is timed on its own clock, so it also includes about 1 ms of closing the stream and prepared object and of output-file bookkeeping outside the three phases.

The 340 edge IDs match across all three. Forced screening did its job and cut the query and CSV phase from about 97 ms to 28 ms. It also spent about 145 ms more on setup, so its export took about 77 ms longer than dense. Faster query, slower job.

For the one-query export, auto chose dense and finished about 60 ms ahead of forced screening. Given the same rows and 100 planned queries, it chooses screening; the measured first-use and query medians put the crossover at query nine. Same engine, same rows, opposite answers, and both of them right. Which of Knuth's two sentences applies depends on how many times you plan to run the job. For arXiv titles at threshold 0.9, the answer is nine.

Appendix A: The memory accounting

The query planner bounds each allocation phase and takes the largest. This is the candidate-dependent part, from _workspace.py:

def query_peak_bytes(cells, candidates):    blocks = (cells + 255) // 256    candidate_bytes = 3 * aligned_bytes(4 * candidates)    scatter = aligned_bytes(4 * cells) + aligned_bytes(8 * blocks) + candidate_bytes    cblocks = (candidates + 255) // 256    classification = candidate_bytes + aligned_bytes(candidates) + aligned_bytes(16 * cblocks)    return max(        count_peak_bytes(cells),        scatter,        classification + scan_scratch_bytes(cblocks),        classification + candidate_bytes,    )

The three candidate arrays hold two int32 indices and one FP32 score per candidate, and every allocation rounds up to the 512-byte pool quantum. The first term is the peak while scores are counted. The second is the scatter phase, where the score tile, the block offsets and the candidate arrays coexist. The third is classification plus the scan that compacts its results. The fourth is the overlap of old candidate arrays and new compacted arrays during compaction. The scan is written out explicitly so its temporary arrays can be counted.

The budget's minimum is 4 KiB. It bounds the library-owned temporary device arrays a query allocates; resident inputs and projections, preparation, the CUDA context and cuBLAS internals, the allocator's cache, host memory and output the caller keeps are accounted separately. The computed bound is conservative: in the A6000 diagnostic at d1a7f1e, forced screening on an isotropic control had a bound of about 401 MiB against a measured peak of about 258 MiB.

Freeing memory had its own twist. When a query raised an exception, the saved traceback's frame locals (Python's record of where the error happened, along with the variables it was holding) could keep device arrays alive after the error. The exception was holding the GPU memory hostage. The cleanup code now clears unwound frames across chained exceptions, and regression tests hold on to the exception while checking that the buffers have been released.

Appendix B: How to check it

On a rented RTX A6000, the first CUDA build passed 111 tests; 53 tests also passed under each of the four Compute Sanitizer tools, NVIDIA's checkers for memory and timing mistakes on the chip. The installed revision 7d46ea5 passed 153 CUDA tests without skips on the RTX 4090. A 42-case subset covering the prepared path, the workspace accounting and the directed-rounding arithmetic passed under each of Compute Sanitizer's memcheck, racecheck, initcheck and synccheck. All 36 campaign workloads passed complete-output checks, and 224 collected artifacts were checked against their recorded hashes. Small fixtures check membership against exact rational arithmetic; on the full corpora, routes are checked for complete agreement with dense. Targeted mutation checks seeded nine faults into the package source, and the tests caught all nine.

On the Mac, the gate's other profiles and dimensions ran 2.448 to 2.924x faster than FAISS.

The code, the preserved per-trial CSVs and the reports are in the repository, under an MIT license.

Appendix C: The break-even arithmetic

For the existing arXiv 0.9 cell, first use already includes the first query, so each further query adds its prepared-query time to the running total:

import math# Existing arXiv 0.9 cell, RTX 4090, revision 7d46ea5. Medians in ms.dense_first, dense_query = 180.30, 64.04auto_first,  auto_query  = 468.37, 25.81first_use_gap = auto_first - dense_first          # first use includes the first querysaving = dense_query - auto_queryrepeats = math.floor(first_use_gap / saving) + 1  # later queries needed to pass denseprint(f"first-use gap {first_use_gap:.2f} ms, saving {saving:.2f} ms per query")print(f"auto is ahead from query {repeats + 1}")for q in (1, 8, 9, 100):    d = dense_first + (q - 1) * dense_query    a = auto_first + (q - 1) * auto_query    print(f"{q:>3} queries: dense {d:8.2f} ms  auto {a:8.2f} ms")

first-use gap 288.07 ms, saving 38.23 ms per queryauto is ahead from query 9  1 queries: dense   180.30 ms  auto   468.37 ms  8 queries: dense   628.58 ms  auto   649.04 ms  9 queries: dense   692.62 ms  auto   674.85 ms100 queries: dense  6520.26 ms  auto  3023.56 ms

Notes

Why not FAISS range_search? FAISS has Index::range_search(n, x, radius, result, params), documented as "return all vectors with distance < radius" and writing into a RangeSearchResult. It is not required of every index type; the documentation says "many indexes do not implement the range_search (only the k-NN search is mandatory)." A normalized flat inner-product index supports batched similarity range searches. What this engine adds is the bounded query scratch memory, the certified early rejection, and the streamed tile-by-tile output.

Why FP32 and not FP16? The project's contract fixes FP32 input, normalization, multiplication and accumulation by default. Reduced-precision multiplication, including TF32 on CUDA, is not part of that contract, and the certificate's error ledger is written for FP32 operands.

Why is the work quadratic? Quadratic means double the rows and you get four times the pairs. An exhaustive join in this engine visits every pair, and screening makes each visit cheaper rather than skipping visits. On the historical million-row arXiv run at 0.9, screening sent about 0.09% of pairs to full scoring, and the prefix comparisons still touched all of them.

Why medians and observed ranges? Each cell has three first-use and nine prepared-query observations per route. The reports give medians, and the charts draw the observed minimum-to-maximum range rather than a confidence interval.

What are the 28 tiles? For 50,000 rows and 8,192-row tiles, the upper triangle is covered by 28 initial tiles of four shapes: six full diagonal tiles, fifteen full off-diagonal tiles, six 8,192-by-848 edge rectangles, and one 848-by-848 edge diagonal.