Accelerating GPU vector retrieval scoring via FP16 accumulation

High-throughput vector search engines and recommendation systems rely on inner-product scoring across query and document embeddings. This operation reduces to a dense matrix multiplication: $Q \in \mathbf{R}^{n_q \times d} \times C^T \in \mathbf{R}^{d \times n_c}$ . On consumer NVIDIA Ampere GPUs (e.g., RTX 3080), PyTorch’s default cuBLAS settings execute this matrix product using single-precision (FP32) accumulation, running at half the peak issue rate supported by the hardware.

By enabling FP16 accumulation via PyTorch runtime options, matrix scoring throughput increases by 1.59× (91.3 TFLOP/s vs 57.5 TFLOP/s). When integrated into an end-to-end retrieval pipeline including top- $k$ selection, net throughput increases by 1.25× while maintaining a Recall@10 score of 0.987 against FP32 ranking baselines.

1.59×

faster scoring matmul

91.3 vs 57.5 TFLOP/s

1.25×

faster end to end

485 vs 608 ms, incl. top-k

0.987

recall@10 vs fp32 ranking

9.87 of 10 results survive

Reduction of Retrieval Architectures to Dense Matrix Multiplication

Dense vector retrieval and recommendation workloads collapse to an identical matrix product:

Brute-force Maximum Inner-Product Search (MIPS): Evaluates every query vector against all corpus embeddings to select the top- $k$ exact matches.
Two-Tower Recommender Systems: Multiplies user embedding matrices by item embedding matrices to rank candidate recommendations.
ANN Index Reranking: Filters corpus embeddings to a candidate shortlist using approximate algorithms (e.g., HNSW), followed by exact matrix dot-product scoring.

Each pattern uses a dense GEMM, $Q \cdot C^T$ , followed by a top- $k$ reduction pass.

Funneling distinct retrieval architectures into a unified dense GEMM scoring kernel.

Brute-force MIPSevery query against every item, exact

Two-tower recsysuser vectors against item vectors

ANN + rerankquery against a candidate shortlist, exact

Q·Cᵀ dense GEMMtensor cores · the flag acts here

Top- $k$ selection is a comparison-based operation that executes on general CUDA cores and does not use Tensor Core accumulation. Consequently, performance profiling must separate matrix scoring throughput from top- $k$ selection overhead.

Mathematical equivalence of vector scoring across retrieval pipelines

Given $n_q$ queries, $n_c$ candidate items, and embedding dimension $d$ , total matrix multiplication work is $2 \cdot n_q \cdot n_c \cdot d$ FLOPs. Top- $k$ selection processes score matrices of shape $(n_q, n_c)$ using heap or quickselect algorithms. Because accumulator precision settings affect only the GEMM phase, top- $k$ latency remains constant across precision modes.

Accumulator Issue Rates in Consumer NVIDIA Ampere GPUs

NVIDIA Tensor Cores perform matrix multiply-accumulate (MMA) operations on FP16 input tiles, adding intermediate products into a running accumulator.

On consumer Ampere hardware (sm_86), accumulator precision determines hardware instruction issue rates:

FP32 Accumulation (fp32-accumulate): Sums dot-product terms in single precision, issuing at half-rate on sm_86.
FP16 Accumulation (fp16-accumulate): Sums dot-product terms in half precision, issuing at full hardware rates.

By default, cuBLAS selects FP32 accumulation to avoid overflow across arbitrary dynamic ranges. However, normalized vector embeddings (e.g., L2-normalized unit vectors) produce bounded inner products, rendering FP16 accumulation safe from numerical overflow.

FP16 accumulation is toggled dynamically via PyTorch:

import torch

# Enable full-rate FP16 accumulation on Tensor Cores (PyTorch >= 2.7)
torch.backends.cuda.matmul.allow_fp16_accumulation = True
scores = queries @ corpus.T
torch.backends.cuda.matmul.allow_fp16_accumulation = False

Hardware architecture differences: sm_86 vs sm_80

Consumer Ampere GPUs (sm_86, e.g., RTX 3080) enforce a 2:1 issue rate penalty for FP32 accumulators relative to FP16 accumulators. Datacenter GPUs (sm_80, e.g., A100) issue FP32 and FP16 accumulators at identical rates. Consequently, performance acceleration from FP16 accumulation is specific to consumer Ampere hardware.

Scoring Throughput and End-to-End Latency Benchmarks

We benchmarked a retrieval pipeline consisting of 8,192 queries, a 2.1-million-item corpus ( $d=512$ , FP16 precision), and top-10 selection per query. The corpus is streamed through GPU memory in 64 MiB tiles (65,536 items/tile) to prevent memory allocation spikes.

Matrix scoring throughput (2.1M corpus, 8192 queries, d=512, RTX 3080)

fp32-accumulate (cuBLAS default)57.5 TFLOP/s · 306 ms

fp16-accumulate (optimized flag)91.3 TFLOP/s · 193 ms · 1.59× speedup

End-to-end pipeline throughput including top-10 selection

fp32-accumulate28.3 G-scores/s · 608 ms total latency

fp16-accumulate35.4 G-scores/s · 485 ms total latency · 1.25× speedup

Scoring pass execution latency (2.1M corpus, 8192 queries, d=512)

fp32-accumulate0.00 ms

fp16-accumulate0.00 ms

Benchmark comparison on NVIDIA RTX 3080 (PyTorch 2.12, CUDA 13)

Accumulation Mode	Scoring Latency	GEMM Throughput	End-to-End Latency	Net Speedup
fp32-accumulate (default)	306 ms	57.5 TFLOP/s	608 ms	1.00×
fp16-accumulate (optimized)	193 ms	91.3 TFLOP/s	485 ms	1.59× (GEMM) / 1.25× (Pipeline)

End-to-end latency includes streaming 32 corpus tiles and merging local top-10 candidates per query.

HBM memory bandwidth bottlenecks during tiled corpus streaming

The scoring pipeline transfers 2.0 GB of corpus embeddings through GPU HBM across 32 streaming tiles. Because memory transfer overheads are not eliminated by instruction issue acceleration, the observed GEMM speedup ( $1.59\times$ ) sits below the theoretical $2.0\times$ issue rate limit.

Ranking Preservation and Recall@10 Validation

Because FP16 accumulation alters intermediate dot-product precision, we evaluated whether score perturbations affect final query rankings. Recall@10 was measured by comparing rankings generated under FP16 accumulation against rankings generated under FP32 accumulation using identical FP16 embeddings.

Across 8,192 test queries, Recall@10 reached 0.987 (an average agreement of 9.87 out of 10 items per query):

Numerical stability of L2-normalized inner products

Relative error under FP16 accumulation is bounded near $\sim 10^{-3}$ . Because score gaps between top-ranked items typically exceed $10^{-3}$ , item rank ordering remains stable, yielding $98.7\%$ top-10 preservation.

Performance Bounds on Unstructured and Structured Sparse Operations

We evaluated whether FP16 accumulation flags accelerate sparse retrieval operations, such as lexical matching or graph adjacency scoring.

Two sparse configurations were benchmarked on an RTX 3080:

Unstructured SpMM (SpMM): Evaluated via cuSPARSE at $1\%$ density.
2:4 Structured Sparsity: Evaluated via CUTLASS routines.

Sparse performance comparison vs dense GEMM ceiling (n=8192, RTX 3080)

Dense FP16-Accumulate GEMM111.7 TFLOP/s · Dense Tensor Core ceiling

Unstructured SpMM (cuSPARSE)0.93 TFLOP/s · 120× below dense ceiling · CUDA Cores

2:4 Structured Sparsity vs matched dense baseline (8192×8192×4096)

Matched Dense Baseline (FP16-acc)114.6 TFLOP/s · 4.80 ms execution latency

2:4 Structured Sparsity (CUTLASS)119.2 TFLOP/s · 1.04× speedup

Sparse operation throughput metrics (NVIDIA RTX 3080)

Sparse Operation	Measured Throughput	Relative to Dense Tensor Core Ceiling
Dense FP16-acc GEMM (8192³)	111.7 TFLOP/s	1.00× (Baseline Ceiling)
Unstructured SpMM (cuSPARSE, 1% density)	0.93 TFLOP/s	120× below dense ceiling
SpGEMM (cuSPARSE, sparse × sparse)	14.8 ms (37.3M nonzeros)	N/A (Memory Bound)
2:4 Structured Sparsity (CUTLASS)	119.2 TFLOP/s	1.04× vs matched dense baseline

Unstructured sparse operations execute on CUDA cores and are limited by memory bandwidth, rendering Tensor Core accumulator flags ineffective.

Arithmetic intensity of sparse matrix execution

Unstructured SpMM at 1% density yields an arithmetic intensity of $\ll 10\text{ FLOPs/byte}$ . Because execution is constrained by HBM memory bandwidth rather than compute issue rates, altering Tensor Core accumulator settings does not affect performance.

2:4 Structured Sparsity performance bounds

2:4 structured sparsity requires packing two non-zero elements per four-element vector into specialized Tensor Core structures. On test shapes, CUTLASS 2:4 execution yielded a modest 1.04× speedup ( $114.6 \to 119.2\text{ TFLOP/s}$ ) over FP16-accumulated dense baselines.

Two-Stage Architecture Design for Retrieval Pipelines

Because accumulator flags accelerate dense GEMMs but do not affect sparse operations, retrieval architectures should separate candidate generation from exact reranking:

Two-stage retrieval pipeline: Sparse candidate generation operates on CUDA cores; dense reranking applies FP16 Tensor Core accumulation.

query

sparse candidate-genCUDA cores · flag: no effect

dense re-ranktensor cores · flag: 1.59×

top-k

results

Stage 1: Sparse Candidate Generation: Filters large corpora using inverted indices or graph traversal (memory-bound, CUDA cores).
Stage 2: Dense Reranking: Scores candidate shortlists using dense embedding GEMMs (allow_fp16_accumulation = True, 1.59× acceleration).

Verification and Benchmarking Protocol

All benchmarks adhered to the following evaluation criteria:

Recall Baseline: Scored against FP32-accumulate rankings on identical FP16 embeddings.
Dual Reporting: Both GEMM-only speedup ( $1.59\times$ ) and end-to-end pipeline speedup ( $1.25\times$ ) are reported to account for top- $k$ overhead.
Timing Events: Measured via asynchronous CUDA events (torch.cuda.Event) across 5 evaluation runs.

Scope, Hardware Constraints, and Future Exploration

Future benchmarks will explore overlapping top- $k$ candidate extraction with tiled GEMM execution to align end-to-end pipeline speedups closer to the 1.59× scoring limit.

Cross-Domain Applicability of FP16 Accumulation

Enabling FP16 accumulation provides immediate acceleration across compute-bound dense GEMM workloads:

Cross-domain acceleration using PyTorch FP16 accumulator flags

Domain	Implementation Pattern	Observed Speedup	Quality Metric Impact
Vector Retrieval & Reranking	PyTorch runtime flag	1.59× GEMM / 1.25× Pipeline	Recall@10 = 0.987
Dense Solves & PDE Green Functions	PyTorch runtime flag	1.84–1.91× GEMM	Relative Error = 2.3e-3 vs FP64
NeRF Neural Field Rendering	Fused CUDA kernel + flag	1.69–1.81× Rendering	PSNR = 35.70 dB

Across all three domains, FP16 accumulation accelerates compute-bound dense GEMMs while maintaining domain quality thresholds.

1 line

of PyTorch

no custom kernel, no rebuild

0.987

recall@10 quality gate

ranking survives the switch

sm_86

consumer Ampere only

A100 not measured here

Conclusion

Setting torch.backends.cuda.matmul.allow_fp16_accumulation = True accelerates consumer Ampere GPU vector scoring by 1.59× and end-to-end retrieval pipelines by 1.25× while preserving Recall@10 accuracy at 0.987. Unstructured sparse retrieval operations receive no speedup, reinforcing a two-stage architecture that isolates dense reranking onto accelerated Tensor Core paths.

Glossary

Glossary — every term, defined

Tensor cores: Dedicated matrix-multiply units on NVIDIA GPUs. They multiply small fp16 input tiles and add up the products far faster than the general-purpose CUDA cores, which is why a dense GEMM belongs on them.
CUDA cores: The general-purpose arithmetic units on the GPU. They run the memory-bound sparse paths (SpMM, SpGEMM) and the top-k merge, all slower than the tensor cores on dense GEMM.
MMA (matrix multiply-accumulate): The tensor-core primitive: multiply two input tiles and add the product into a running accumulator. A GEMM is built from many MMAs.
Accumulator: The running sum a matmul adds each product into. Its precision (fp16 or fp32) is chosen separately from the input precision, and it is the single knob this flag turns.
fp16 / fp32 accumulation: Summing the running dot product in half or single precision. On sm_86 GeForce the fp16 path issues at twice the rate of the fp32 path for identical fp16 inputs.
GEMM: General matrix-matrix multiply, the dense C = A·B kernel that every retrieval and recsys pattern here reduces to.
cuBLAS: NVIDIA's dense linear-algebra library and the GEMM backend PyTorch calls. It defaults to fp32 accumulation because that is correct on any input.
cuSPARSE: NVIDIA's sparse linear-algebra library. It runs the sparse multiplies (SpMM, SpGEMM) on the general-purpose CUDA cores, never the tensor cores, so the accumulator flag does nothing for them.
CUTLASS: NVIDIA's open-source library of CUDA matrix-multiply templates. It exposes the tensor-core path for 2:4 structured sparsity, the one sparse layout cuBLAS does not.
recall@10: The share of each query top-10 that survives the switch to fp16 accumulation, scored against the fp32-accumulate ranking of the same embeddings. Here it is 0.987, so 9.87 of every 10 results are unchanged.
Arithmetic intensity: Floating-point operations per byte of memory traffic. High intensity is compute-bound on the tensor cores; low intensity is memory-bound on the CUDA cores, which is the split between dense scoring and sparse search.
Roofline: A model plotting achievable throughput against arithmetic intensity. Its ridge marks where a kernel flips from memory-bound to compute-bound.
Compute-bound: Limited by how fast the arithmetic units retire operations rather than by memory bandwidth. This is the regime where the accumulator issue rate binds and the flag lands.
Memory-bound: Limited by memory bandwidth rather than arithmetic rate. Unstructured sparse work sits here, so the accumulator flag has nothing to accelerate.
HBM: High-bandwidth memory, the GPU's on-package DRAM. Streaming the 2 GiB corpus through it is the memory traffic that holds the scoring pass below a purely compute-bound 2×.
SpMM: Sparse-times-dense matrix multiply, the lexical / TF-IDF pattern. It runs through cuSPARSE on the CUDA cores at 0.93 TFLOP/s, 120× below the dense tensor-core ceiling.
SpGEMM: Sparse-times-sparse matrix multiply. Its FLOP count depends on the product pattern, so no throughput rate is defined; the benchmark records 14.8 ms and 37.3M output nonzeros instead.
2:4 structured sparsity: A pattern with exactly two nonzeros in every group of four values, the one sparse layout Ampere tensor cores accept directly through CUTLASS. On this shape it measured 1.04× over a matched dense baseline.
top-k: Selecting the k highest-scoring items per query after the GEMM. It is a comparison, not a matmul, so the flag never speeds it up, which is why the end-to-end win is 1.25× against the scoring 1.59×.
MIPS: Maximum inner-product search: score every query against every item and take the exact top-k. One of the three retrieval patterns that reduce to the same dense GEMM.
ANN: Approximate nearest-neighbor search: narrow to a candidate shortlist with an approximate index, then rerank that shortlist exactly with the dense GEMM.