GPUs: where the bytes go

On 16 October 2026, on NVIDIA A100 GPUs at the National Energy Research Scientific Computing Center, a search that the benchmark’s fastest conventional kernel finished in 1.79 seconds was finished in 0.30 seconds. The workload was synthetic but realistic in shape: 1,024 exposures of a 256 × 256 pixel patch, searched over 4,225 trial velocities for a source too faint to see in any one frame. Both programs computed the same numbers to within half-precision rounding. The faster one did it by rewriting a sum over shifted images as thousands of small matrix products and handing them to the tensor cores, the matrix units that modern GPUs carry for machine learning.

A factor of 5.9 is a satisfying number, and it is the wrong one to remember. The slower program was not running anywhere near the limit of the hardware. A conventional kernel that simply streamed the needed data at the card’s memory bandwidth, with no clever reuse at all, would have taken about half a second. Measured against that, the tensor-core kernel is about 1.7 times faster. And the workload, in which every one of the 1,024 exposures covers every output pixel, is the single most favourable case for the method. In a sparse archival search the advantage disappears.

This chapter is about why both statements are true. A GPU accelerates a planet search only when the search supplies enough independent work, keeps data close to the arithmetic, and preserves the detection statistic. The useful question is therefore not which device advertises the most operations per second. It is which representation and schedule minimize the cost of a scientifically equivalent search. We start with the machine, then the roofline model that explains most of what follows, then tensor cores and precision, then the kernels astronomers actually use. The last part is the A100 case study, told with its caveats attached.

What a GPU is, for a stacking problem

Lanes, warps and multiprocessors

A GPU is a large number of simple arithmetic lanes, grouped so that one instruction drives many of them at once. NVIDIA’s terms are now standard. A thread is one lane’s sequence of work. Thirty-two threads that execute the same instruction together form a warp. Warps are scheduled on a streaming multiprocessor (SM), which holds the arithmetic units, a large register file and a small, fast, programmer-managed scratchpad called shared memory. The A100 has 108 SMs; each can keep up to 64 warps (2,048 threads) resident at once and owns 65,536 32-bit registers, with at most 255 registers per thread (NVIDIA 2020).

The design hides latency by oversubscription. A read from the GPU’s main memory takes hundreds of cycles. Rather than wait, the SM switches to another warp whose operands are ready. Occupancy, the fraction of the 64 warp slots actually filled, measures how much of this latency-hiding the kernel allows. It is limited by whichever resource runs out first: registers, shared memory or thread slots. A kernel that needs 255 registers per thread and launches blocks of 256 threads consumes $255\times256=65{,}280$ registers per block, essentially the whole register file, so only one block of eight warps fits on each SM: an occupancy of $8/64=12.5\%$. We will meet exactly this kernel in the case study.

The memory hierarchy

Data live at several distances from the arithmetic. Off-chip high-bandwidth memory (HBM) holds the full data set: 40 or 80 GB on an A100. A shared L2 cache (40 MB on the A100) sits between HBM and the SMs. Each SM has up to 164 KB of shared memory and 256 KB of registers (NVIDIA 2020, Table 4). Each step inward is roughly an order of magnitude faster and an order of magnitude smaller. Fu and colleagues measured, for their cost model on an A100-40GB, an achieved HBM bandwidth of 1.35 TB/s, shared-memory bandwidth of 9.5 TB/s, 234 TFLOP/s for tensor-core matrix multiplication and 17.6 TFLOP/s for general arithmetic (Fu et al. 2023, Table 19). Those four measured constants already contain the plot of this chapter: data must be reused many times on-chip before the matrix units can be kept busy.

The manufacturer peak specifications for the three generations most relevant here are below. They are ceilings, not measured astronomical throughput.

Sources: NVIDIA A100, NVIDIA H100, NVIDIA HGX and DGX B200 specifications. Sparse tensor rates are halved to give dense rates, and B200 values divide the eight-GPU system figures by eight, following the research for this book. The 40 GB A100 SXM4 has 1,555 GB/s of HBM bandwidth rather than 2,039 GB/s (NVIDIA 2020, Table 4), a difference that matters in the case study. The ridge point, defined in the next section, is peak arithmetic divided by bandwidth (estimate).

Coalescing

HBM delivers data in aligned chunks. When the 32 threads of a warp read 32 adjacent words, the hardware combines them into a few wide transactions: the access is coalesced. When they read 32 scattered words, each may cost its own transaction, and most of each transferred chunk is wasted. NVIDIA’s best-practices guide treats coalescing, staging through shared memory and occupancy as the first-order optimizations for memory-bound code (NVIDIA, CUDA C++ Best Practices Guide). For an image search the natural mapping is spatial: neighbouring threads handle neighbouring starting pixels, so that at each exposure they read neighbouring pixels. This single decision, made in KBMOD in 2019, is most of what separates a good stacking kernel from a bad one.

The roofline: two different bottlenecks

Let $F$ count the useful operations a kernel performs and $Q_H$ the bytes that actually cross HBM, after whatever reuse the caches and shared memory provide. The operational intensity is $I_H=F/Q_H$. A kernel cannot compute faster than the arithmetic units allow, nor faster than memory can feed it:

where $\beta_H$ is HBM bandwidth. Plotted on logarithmic axes, the bound is a sloped line (memory-bound) meeting a flat roof (compute-bound) at the ridge point $I^\ast=P_{\rm peak}/\beta_H$. Williams, Waterman and Patterson, who introduced the model, stress two details that are easy to forget: the bytes are traffic after caches, and the bandwidth should be the sustainable rate, not the pin rate (Williams et al. 2008, §3). Real kernels also hit lower ceilings from shared-memory bandwidth, instruction issue, type conversion, address arithmetic, register spills and latency.

Now place the two operations that recur in this book. A trajectory sample in a likelihood stack reads two FP32 numbers, the signal statistic $\psi$ and the information statistic $\phi$, and adds each to a running sum. That is two useful floating-point operations for eight bytes: $I_H=0.25$ FLOP/byte before any reuse (estimate). On an A100 the memory roof at that intensity is $0.25\times2.04\approx0.51$ TFLOP/s, on an H100 about 0.84 TFLOP/s. The FP32 arithmetic peak is irrelevant. Calling this a low-FLOP kernel is misleading: additions are cheap, and supplying the selected values is the work. The sample also needs coordinates, a bounds check and perhaps interpolation, so the real instruction count is several times the useful count.

A square matrix product is the opposite. Multiplying two $n\times n$ FP16 matrices into an FP32 result performs about $2n^3$ operations while, ideally, reading $2\times2n^2$ bytes and writing $4n^2$: $I_H\approx n/4$ (estimate). At $n=1024$ the intensity is 256 FLOP/byte, above the A100 ridge of 153 and near the H100 ridge of 295. Tiling the product through shared memory and registers is what makes that ideal intensity attainable.

Two consequences follow immediately. First, moving an unreused gather onto tensor instructions does nothing: it stays three orders of magnitude to the left of either ridge. Second, the hardware improvement for a memory-bound search is the bandwidth ratio, not the arithmetic ratio. For $10^{12}$ paired samples at eight bytes each, the bandwidth floor is 2.39 s on an H100 and 1.00 s on a B200 (estimate): a factor of 2.4, not of ten. Cache reuse among neighbouring trajectories can push the actual HBM bytes below eight per sample, and coordinate work and uncoalesced transactions push the time above the floor.

Capacity can matter more than speed. A cube of 1,000 images of $4096^2$ pixels needs about 134 GB as FP32 $(\psi,\phi)$ pairs, 67 GB as FP16 pairs, or 34 GB as eight-bit pairs (estimate, excluding metadata and working memory). On an 80 GB card, the choice of storage format decides whether the search runs in one pass or must be tiled through host memory, a discontinuous change that can outweigh any twofold instruction speedup.

Tensor cores

A tensor core is a fixed-function unit that computes a small matrix product and accumulate, $D=AB+C$, in one instruction. NVIDIA introduced them with the Volta V100; each Volta tensor core performs a $4\times4\times4$ FP16 product with FP32 accumulation per clock. Each of the A100’s 432 tensor cores performs an $8\times4\times8$ product per clock, which is where its 312 TFLOP/s of dense FP16 comes from: sixteen times its 19.5 TFLOP/s of ordinary FP32 (NVIDIA 2020). Programmers do not address individual units. A warp collectively issues a matrix instruction on fragments held in its registers.

Three practical facts follow from this design.

The inputs are low precision, the sums need not be. FP16 has 10 explicit significand bits and a narrow exponent range. BF16 keeps the FP32 exponent range with only 7 significand bits. TF32 is a 19-bit format with FP32 range used internally by Ampere and later. INT8 and INT4 codes are exact integers with 32-bit integer accumulation, and Hopper adds FP8. The common pattern for science is FP16 or BF16 operands with FP32 accumulation, which keeps the running sum accurate while halving operand bytes.

Complex arithmetic has to be embedded in real arithmetic. For $A=A_r+iA_i$ and $B=B_r+iB_i$,

so one complex product costs four real tensor products, or one product of a real block matrix with duplicated and sign-flipped entries. The extra executed arithmetic is worthwhile because it runs on the fast unit. Honest accounting counts useful operations (eight real FLOPs per complex multiply-accumulate), not every padded multiply the hardware performs (Oostrum et al. 2025, §§III-IV).

Shape and interface matter. Matrix fragments come in fixed sizes, so small or ragged products are padded with zeros, and the padding is executed arithmetic that does no science. On Hopper, Oostrum and colleagues found that the older warp-level WMMA interface reaches only about 60-65% of the tensor ceiling; the larger ceiling needs the warp-group WGMMA interface. They also found that one-bit XOR is compiled inefficiently on Hopper (2,361 TOP/s measured on GH200, against 10,276 TOP/s for AND), and rewrote signed one-bit products in terms of AND (Oostrum et al. 2025, §III, Table I). Peak numbers presuppose that the kernel uses the right instruction on the right shape.

The statistic the hardware must preserve

Before looking at kernels, fix what they must compute. Chapter 1 derives the matched filter. In the form used by stacking codes, a background-subtracted exposure $e$ is modelled as

with pixel $p$, sampled PSF $P_e$, predicted position $q_e(\theta)$ for trajectory parameters $\theta$, known relative throughput $a_e$, unknown common flux $f$, and independent Gaussian noise of variance $V_e(p)$ (masked pixels get inverse variance zero). Precompute, once per exposure,

Then for any trajectory

The flux estimate is $\hat f=A/B$, and under the ideal null hypothesis $Z$ is standard normal (Whidden et al. 2019, §2.1). Everything a GPU does in a stacking search is in service of evaluating $A$ and $B$ for an enormous number of $\theta$. Two points matter for accelerators. The sums must stay signed: replacing them with powers or per-segment maxima changes the detector. And $B$ must be the true variance of $A$. A kernel that changes interpolation, masks or precision without changing $B$ accordingly is fast at computing a miscalibrated statistic.

For equal independent exposures the SNR grows as $\sqrt E$, so the ideal depth gain is $1.25\log_{10}E$ magnitudes: 3.75 mag for $E=1000$ (derived). No kernel can exceed the information in the photons. What a kernel can do is make a longer coherent integration, or a larger trajectory bank, affordable.

Four families of astronomical kernels

Likelihood gathers: KBMOD and DEEP

Peter Whidden, Bryce Kalmbach, Andrew Connolly and collaborators built Kernel-Based Moving Object Detection, KBMOD, around the equations above (Chapter 17 tells its story). Its decisive implementation choice was to precompute $\psi$ and $\phi$ for each exposure, interleave them in memory, and inspect trajectories without ever building a shifted image. Each GPU thread owned one starting pixel and one velocity, and adjacent threads owned adjacent pixels in $32\times16$ blocks, so reads coalesced. The paper reports a runtime coefficient near $2.4\times10^{-11}$ s per trajectory per image and a search of 240 billion trajectories through 30 images in about 180 s, on a GPU whose model that passage does not name (Whidden et al. 2019, §2.3). That is about $4.2\times10^{10}$ paired samples per second (derived). Applied to the 2015 HiTS DECam data, the method found 39 previously unknown Kuiper belt objects over 150 square degrees.

The minimal exhaustive gather is short:

for each spatial tile and motion packet, in parallel:
    neighbouring origin pixels -> neighbouring lanes
    for each motion in the packet:
        A, B = 0
        for each exposure e:
            q = predicted_position(origin, motion, e)
            if q is inside the image and valid:
                psi, phi = decode(stats[e, q])
                A += a[e] * psi
                B += a[e]^2 * phi
        if B > 0: score = A / sqrt(B); keep if above threshold

For $P_0$ origins, $H$ motions and $E$ exposures this costs $O(P_0HE)$ samples plus coordinate work, against $O(EP)$ storage. The pinned current KBMOD source uses $128\times2$ thread blocks, supports FP32 and scaled UINT8/UINT16 encodings with a reserved invalid code, computes trajectory times in double precision, and has a compile-time limit of 448 images (KBMOD source, rev. f4b16e4). The natural improvement is to batch nearby motions as well as nearby origins: load the union of their source pixels for one exposure into shared memory once, decode once, and update several sums per thread. Packets whose trajectories diverge overfetch, and each extra motion costs registers, so the best packet size is an empirical question.

The search kernel is only part of a discovery pipeline. The DECam Ecliptic Exploration Project (DEEP) used about 100 exposures per CCD and a $50\times30=1{,}500$ motion grid. Its median time per CCD was 576 s on an NVIDIA A40, including disk loading, the GPU search, false-positive filtering and stamp construction, and total search runtime was about 282 hours. Reverse-direction searches measured the artifact background, and dense or slow-motion fields required raising the threshold from 7 to 10 (Smotherman et al. 2023, §3). A new kernel must be timed alongside those stages, and its sensitivity compared with their selection effects intact.

Whole-image warps: heliostack

Kevin Napier, Matthew Holman and collaborators’ heliostack organizes the work differently (Chapter 20). It asserts orbital parameters, maps each image’s corners to a common barycentric frame, approximates the mapping by a perspective transform (errors below 0.2 pixels), and warp-adds whole images of signal and information onto a canvas with nearest-neighbour sampling. Geometry is shared across a whole image instead of being solved per pixel per trajectory. Reprocessing 2003 Hubble ACS data, it searched 1,118 image planes over 1,148,817 asserted initial conditions on a roughly $14{,}000\times13{,}000$ pixel canvas, and reported two recovered and two new faint objects at measured ST magnitudes $28.99\pm0.10$ and $29.17\pm0.11$. The whole analysis took a few hundred H200 GPU-hours plus a few thousand CPU-hours (Napier et al. 2026, §§2-4). That aggregate covers every stage and should not be divided by a nominal trial count to make a kernel benchmark.

Correlators and beamformers on tensor cores

Radio astronomy reached dense linear algebra first, because its central operation already is one. An interferometer correlates the complex voltages of $R$ receivers over $T$ samples in each frequency channel:

Each voltage is reused in $R$ products, which is precisely the reuse a gather lacks. Michael Clark, Paul La Plante and Lincoln Greenhill showed in 2011 that explicitly tiling voltages through the memory hierarchy, rather than trusting the cache, lifted a GTX 480 correlator to 1,058 GFLOP/s, 79% of single-precision peak, where earlier GPU correlators had reached roughly 10-30% (Clark et al. 2013, §§1, 6.1).

John Romein’s Tensor-Core Correlator (TCC), presented at an ALMA workshop in February 2020 and published in 2021, embedded the complex Hermitian product in tensor-core instructions, with separate kernels for square and diagonal-triangle blocks of $64\times64$ dual-polarization receivers (Romein 2020; Romein 2021). The maintained library reports five- to tenfold efficiency gains over regular-core correlators, an author-reported comparison (TCC repository). Its README also documents the unglamorous requirements. Input must arrive in a specific layout with eight time samples per block for FP16, sixteen for INT8 and FP8, and 32 for INT4, and transposing the stream just before correlation can consume the saving. The most negative integer code is excluded because conjugation cannot negate it. A change of 10 September 2026, “Avoid visibility overflows with i8 input”, added saturating arithmetic: saturation prevents wraparound but cannot preserve an out-of-range sum, so accumulator bounds remain part of scientific correctness.

Leon Oostrum, John Romein and colleagues’ Tensor-Core Beamformer (2025) generalizes the approach to beamforming, $Y=WX$, with a portable complex GEMM library for NVIDIA and AMD hardware.

These are kernel measurements on large tuning matrices ($M=N=K=8192$ for FP16), with useful operations counted as $8MNK$ (Oostrum et al. 2025, Table III). The A100 figure is 55% of its dense FP16 peak. For LOFAR-like beamforming the library was up to 20 times faster and 10 times more energy-efficient on an A100 than the existing FP32 kernel. The roofline study showed small matrices to be bandwidth-limited and large ones reaching about 50-85% of measured tensor ceilings. The lesson for stacking is direct: a large synthetic GEMM is not the performance of a skinny search contraction, so measure the real shapes. In September 2026 the EDD backend paper described a tensor-core correlator-beamformer deployed at the ARGOS prototype (Barr et al. 2026, §5).

Dedispersion and FFT convolution

Incoherent dedispersion (Chapter 36) is the closest radio analogue of a trajectory gather. It sums a dynamic spectrum along curves $\Delta t(\nu,d)=k_{\rm DM}d(\nu^{-2}-\nu_0^{-2})$, at a cost $O(N_dN_tN_\nu)$ for $N_d$ trial dispersion measures, $N_t$ times and $N_\nu$ channels. Novotný and colleagues stage time-frequency tiles in shared memory and keep several outputs per thread in registers. On an A100 they report about 8.5 TFLOP/s useful for four-bit input and 5 TFLOP/s for eight-bit, at least ten times faster than Heimdall in the studied low-bit cases. Their limiting resources were shared-memory service at low bit depth and about 95% conversion-unit utilization in the sixteen-bit path, not HBM (Novotný et al. 2023, §§5-6). These are ordinary GPU kernels, not tensor-core results.

Fourier-domain dedispersion turns the shifts into phases, $\widetilde D(d,f)=\sum_\nu e^{2\pi i f\Delta t(\nu,d)}\widetilde I(\nu,f)$. At each Fourier frequency this is a matrix-vector product, not automatically a GEMM. Bassa and colleagues’ Titan RTX implementation, in ordinary FP32, was about 20% faster than their optimized time-domain code once roughly 512 or more trial DMs were needed. They left tensor cores to future work (Bassa et al. 2021, §§3-5). This investigation found no primary measured tensor-core dedispersion result suitable for quoting an equal-sensitivity gain.

Fused FFT convolution is the third route. In 2018 Sofia Dimoudi, Karel Adámek and colleagues kept short FFTs and filtering inside shared memory for pulsar acceleration searches and measured 1.5-3.2 times their cuFFT-based application on Maxwell GPUs (Dimoudi et al. 2018, §5). FlashFFTConv (2023) adds tensor cores: it factorizes the discrete Fourier transform into small dense matrix products, $F_N=P(I_{N_2}\otimes F_{N_1})D\,P^{-1}(I_{N_1}\otimes F_{N_2})P$ for $N=N_1N_2$, and fuses the whole convolution on-chip up to length 32K. A fixed-order decomposition does more arithmetic than a radix FFT, $O(N^{1+1/p})$ against $O(N\log N)$, yet runs faster because the arithmetic lands on tensor cores. On an H100 its plain forward convolution beat PyTorch by up to 6.54 times; the often quoted 7.93 times is for gated convolution, where adjacent elementwise operations are fused too (Fu et al. 2023, Tables 3-4). These measurements are of machine-learning convolutions, not PSF filtering or orbit stacking.

Precision is a detection-theory question

Lower precision halves or quarters the bytes per sample, and the GPU rewards that. The question is what it costs in SNR. Suppose each statistic is stored with independent, unbiased rounding of step $\Delta_e$. The rounding adds variance $\Delta_e^2/12$, and if the signal response is unchanged the retained SNR fraction is

Bernstein and colleagues derived this codec-noise relation for image compression, and Pence, White and Seaman showed why the quantization step should follow the noise rather than the raw dynamic range (Bernstein et al. 2010, §III; Pence et al. 2010, §II). The fractional penalty does not shrink with more exposures, because original noise and independent rounding noise average down together.

Real codes also saturate. For a unit-variance statistic $X$ and a symmetric quantizer $Q$, the response to a small mean signal is $E[XQ(X)]$ and the output noise is $\sqrt{E[Q(X)^2]}$, so

This is exact in the weak-signal limit and includes clipping. Evaluated for signed codes $-m,\dots,m$ with $m=2^{b-1}-1$, saturating at $\pm c\sigma$, it gives these analytic estimates (derived for this book):

Depth loss is $-2.5\log_{10}\eta$. Eight bits scaled to the noise cost almost nothing. The same eight bits scaled to a bright residual at $100\sigma$ cost 2.5% of the SNR, because the faint-noise steps become coarse. Keeping only the sign gives $\eta=\sqrt{2/\pi}=0.798$, a 20% loss (derived). Radio correlators process one-bit voltages efficiently, but that does not license sign-only storage of optical likelihood maps.

Storage precision and accumulation precision are separate choices. FP16 represents normalized noisy maps adequately; a running FP16 sum over a thousand exposures does not, because late additions fall below the accumulator’s spacing. A positive information sum $B$ is especially exposed. The cuFFT guide warns that the zero-frequency output of a half-precision transform sums all inputs and can overflow (NVIDIA, cuFFT guide, §2.3). White and colleagues moved the convolution stage of a pulsar acceleration search to BF16 while keeping the first long FFT in FP32 for exactly this reason. In about 20,000 synthetic files, 14,602 of 14,971 bright fundamental peaks (97.53%) stayed in the same bin and the other 369 moved to an adjacent acceleration bin; peak SNR changed by at most about 3%, and the speedup was about 1.6 (White et al. 2022, §§3-5). That is careful empirical evidence, but their SNR is peak power over plane standard deviation, so it does not certify a matched-filter efficiency or an unchanged false-alarm rate.

This project’s benchmark kit measured precision directly on its synthetic stacks, on CPU with operand rounding emulated. Over 512 independent paired noise and injection cubes, INT8/UINT8 storage of $(\psi,\phi)$ with per-image noise scaling retained 0.99949 of the FP32 gather’s effective SNR (95% interval 0.99815-1.00084), FP16 storage 0.999998, and the FP16 separable tensor-core operator of the case study 1.000002, indistinguishable from exact. These are precision results on Gaussian fixtures. They are not survey completeness measurements and not calibrations of an eight-sigma tail.

Case study: tensor-core velocity stacking on an A100

The kernels above give a clear rule: tensor cores pay when the computation has dense reuse. A shift-and-stack search over a grid of velocities seems to have none, because each trajectory touches one pixel per exposure. This section describes a reformulation that creates the reuse, a benchmark of it on real hardware, and what the numbers do and do not show. The measurements were made for this book on NERSC Perlmutter on 16 October 2026. They are unpublished and have not been independently reproduced.

The separable-phase identity

Take linear motion plus any known common curve (parallax at an assumed distance, for instance), on a Cartesian grid of velocities: $V_a$ values $v_{x,a}$ along one axis and $V_b$ values $v_{y,b}$ along the other. The displacement of exposure $i$ is then

so the $x$-displacement depends only on the $x$-velocity and the $y$-displacement only on the $y$-velocity. The stack for velocity $(a,b)$ is $A_{ab}(\mathbf x)=\sum_i\psi_i(\mathbf x+\mathbf s_{ab,i})$. Take the spatial Fourier transform with $\hat\psi(\mathbf k)=\sum_{\mathbf x}\psi(\mathbf x)e^{-i\mathbf k\cdot\mathbf x}$. A shift becomes a phase, and because $\mathbf k\cdot\mathbf s=k_xs_x+k_ys_y$ the phase factorizes:

where $(\mathsf A_{k_x})_{ai}=e^{ik_xs_x[a,i]}$ is a $V_a\times E$ phase matrix, $(\mathsf B_{k_y})_{bi}=e^{ik_ys_y[b,i]}$ is $V_b\times E$, and the transpose is a plain transpose, not a conjugate. For each spatial frequency, the whole velocity plane of the stack is one complex matrix product whose inner dimension is the number of exposures. The same operator applied to $\hat\phi$ gives the information. An inverse FFT per velocity, a crop of the padding, and $Z=A/\sqrt B$ finish the job.

Two prior bodies of work frame the novelty. Fourier velocity filtering of image sequences dates at least to Reed, Gagliardi and Shao in 1983 (Reed et al. 1983), and tensor-core complex GEMM is established in radio astronomy. The project’s design study searched for this exact Cartesian factorization (each exposure contributing a separable rank-one matrix, equivalently Khatri-Rao columns in the full phase matrix) applied to optical asteroid stacking, and did not find it. That is a bounded search result, not a claim of priority.

What the identity saves, and what it does not

The identity saves phase generation and storage. The full phase tensor has $E\cdot F\cdot V_aV_b$ entries for $F$ spatial frequencies; the two axis tables have only $E(V_aW_f+V_bH_f)$, where $W_f$ and $H_f$ are the padded transform widths. In the design study’s example of a $128^2$ transform and $32\times32$ velocity blocks, that is a thousandfold reduction in phase work.

It does not reduce the sum. Each frequency still performs $E\cdot V_aV_b$ complex multiply-accumulates per channel, so the useful work is about $16FEV$ real FLOPs for both channels, with $V=V_aV_b$. A direct gather does $P\cdot V\cdot E$ paired samples for $P$ retained origins. Since $F$ exceeds $P$ by the padding needed for the largest displacement, the Fourier route does more useful arithmetic per sample than the gather’s two additions. It can win only by moving the work from a memory-bound gather onto units sixteen times faster than ordinary FP32. Then there are the costs outside the GEMM: a forward FFT of every exposure (amortized over all velocities), an inverse FFT per velocity and channel, and the crop and scoring.

The method therefore needs density. Every exposure must overlap every output pixel of the tile, or the GEMM multiplies zeros. The project’s design model put the gain for 1,000 fully overlapping epochs at 7.3-24 times a strong modelled H100 gather, falling to 2.9-9.6 times if general mask normalization costs 2.5 times more. Applied to heliostack’s real Hubble footprint, where an average canvas output is covered by about 51.5 useful epochs rather than 1,118, the same model gave only 0.30-1.11 times. These are projections, not measurements, and they define the regime in which the experiment below was run.

The benchmark

The benchmark kit generates synthetic stacks with Gaussian noise, spatially varying variance, random masked holes, a moving bad column, variable seeing ($\sigma$ from 1.0 to 1.65 pixels), irregular exposure times and a common quadratic-plus-periodic track resembling parallax. One source is injected at stack SNR 8, below single-exposure detectability. For timing, tracks are rounded to integer pixels and the Fourier methods are given exactly the same rounded shifts, so every method evaluates the same discrete operator and their outputs can be compared number by number. Each method returns a digest: the sums of $A$ and $B$, the maximum score, and the count of scores above a common threshold set by a Bonferroni correction over the full bank. Every score is computed and tested; nothing is dropped.

The kit’s conventional kernels are an auditable serial gather in KBMOD’s style ($128\times2$ blocks, one origin per thread, serial velocity and epoch loops, FP64 coordinates, FP32 sums), FP16 and INT8 variants of it, and a tiled gather that owns a $32\times8$ block of origins and up to 32 neighbouring velocities, stages the union of their source pixels in shared memory, decodes once and updates all sums from there. It also has a Fourier path (cuFFT plus a dense phase sum) and a separable GEMM path built on cuBLAS.

First run: the gather is not bandwidth-bound

The first run used the kit as designed on an A100-SXM4-80GB. All eleven paths agreed with the CPU reference (FP32 to about $10^{-7}$, the FP16 GEMM to about $10^{-3}$ RMS). The tiled gather was fastest, at 42, 106 and 135 billion paired samples per second for 100, 448 and 1,024 exposures, 1.1-2.0 times the serial KBMOD-style kernel. For scale only, that is about three times the 2019 KBMOD coefficient, on unlike hardware.

The surprise was precision. FP16 and INT8 storage changed gather speed by at most 1.1 times. Halving or quartering the bytes did nothing because bytes were not the bottleneck: the serial FP32 kernel moved a logical payload of only 300-660 GB/s against about 2 TB/s available. Its time went to double-precision coordinate arithmetic for every sample, plus branching. On an A100 FP64 runs at half the FP32 rate. A gather that recomputed coordinates in FP32 from precomputed offsets would very likely be faster, and building it is the first item on the project’s list.

The kit’s GEMM path, with $8\times8$ velocity blocks, reached only 0.2-0.34 TFLOP/s and was 5-12 times slower than the gather. Small blocks underfill the tensor cores.

Why the library GEMM lost

The second run enlarged the problem to 1,024 exposures, $256^2$ pixels and a $65\times65$ grid of 4,225 velocities, so that the velocity axes could fill tensor-core tiles. The tiled gather took 1.79 s per scan (about 145 billion samples per second) on the A100-SXM4-40GB the job received. The kit’s cuBLAS GEMM path took 5.60 s with 32-velocity blocks and 4.63 s with 64-velocity blocks, and its Fourier path 5.2 s.

The GEMM itself was not the problem. A library GEMM multiplies matrices that already exist in memory, and the operand here does not: it is the phase table multiplied elementwise by each frequency’s spectrum, $\mathsf B_{k_y}\operatorname{diag}(\hat{\boldsymbol\psi}(\mathbf k))$, different for every frequency. Materializing it means writing and re-reading a tensor of size $F\cdot E\cdot V_b$, which turns a compute-bound product back into a memory-bound one. A version that materialized the operand once per chunk and called batched cuBLAS took 1.73 s, with the contraction at 3.7 TFLOP/s. The diagonal scaling has to happen in registers, as a prologue to the product, and no off-the-shelf GEMM call does that.

A fused kernel in Triton

The fused kernel (tc_stack.py) was written in Triton, a Python-embedded language that compiles tiled array programs to GPU code (Tillet et al. 2019). One program instance handles one spatial frequency and one block of each velocity axis. It loops over exposures in chunks; for each chunk it loads the four real components of the two spectra (signal and information), loads or computes the corresponding phase-table tiles, forms the scaled operand $\mathsf B\circ\hat\psi$ in FP32 registers and casts it to FP16, and issues four real tensor-core dot products per channel with FP32 accumulation. To keep FP16 in range, each frequency’s spectrum is divided by its largest magnitude before the product and multiplied back afterwards. Then PyTorch’s inverse FFT, a crop and the same scoring and digest as the kit complete the scan.

The kernel went through five versions on an A100-SXM4-80GB, each timed per complete scan:

Kernel TFLOP/s counts only useful operations, $16FEV$, over kernel time, and kernel time includes building the per-tile phase tables. The first step is an old lesson: transcendental functions are slow, and a phase table that is reused across all frequencies sharing a $k_x$ is cheap to precompute. The last step came from the profiler. Before it, Nsight Compute showed 255 registers per thread with spills to local memory, 12.5% occupancy, the tensor pipe active only 20% of the time, and stalls dominated by long-scoreboard waits (warps waiting on memory loads). The cause was four $64\times64$ FP32 accumulator tiles per program, real and imaginary parts for two channels, alongside the operand fragments. Splitting the two channels into separate programs halved the accumulators, smaller $64\times32$ output tiles halved them again, and four warps per program let more programs fit on each SM.

The best kernel sustains about 27 TFLOP/s, roughly 9% of the A100’s 312 TFLOP/s dense FP16 peak. For comparison, the Tensor-Core Beamformer reaches 173 TOP/s (55% of peak) on the same card with $8192^3$ matrices. Our products are much smaller, with an inner dimension of 1,024 and output tiles of 64 by 32. With the kernel this fast, post-processing (inverse FFT, crop, score) takes 0.051 s per scan, about 17% of the total. The remaining few hundredths of a second are scoring reductions and host overhead outside the two timed regions.

One derived number shows what kind of win this is. The best kernel performs about $6\times10^{12}$ useful operations per scan, about 24 for each of the roughly $2.6\times10^{11}$ paired samples a gather would read (derived from the measured rates). That is twelve times the gather’s two useful additions. The tensor route wins by doing much more arithmetic on much faster units, which is the opposite of algorithmic savings.

Correctness

The fused kernel’s digest gave a maximum score of 8.6313 with 9,331 trials above threshold. The FP32 paths gave 8.6311 and 9,329. Most of those threshold crossings are neighbouring trajectories of the injected source, and the small differences come from rounding the phase tables to FP16. The project’s CPU emulations of the same operator found FP16 score errors of a few $10^{-4}$ RMS and template SNR retention indistinguishable from one. The SNR-retention test of the GPU kernel itself, with paired injections and null banks, has not yet been run. Until it has, “matches FP32 to the fourth decimal on one digest” is the correct description.

How much faster, honestly

The experiment does establish three things. The separable-phase reformulation runs on real tensor cores, at a useful fraction of their throughput, and reproduces the exact gather’s scores to half-precision accuracy. The decisive implementation detail is fusing the per-frequency diagonal into the operand load, which no library call does. And on dense long stacks the tensor-core route is at least competitive with the best gather one could reasonably build, and probably somewhat faster.

It does not establish an order-of-magnitude gain at equal sensitivity, or anything about a real survey. The project’s own next steps are, in order: an optimized gather (FP32 coordinates, precomputed offsets, velocity-batch reuse) as the fair denominator; a CUTLASS-class kernel with a fused post-processing stage; the SNR-retention and null-calibration test of the FP16 path; and real footprints with masks, where normalization (below) becomes the hard part. The project’s design study also pre-registered an H100 test: the contraction must sustain at least 99 TFLOP/s at 1,000 exposures to clear a tenfold margin over the strongest modelled gather. That is about 10% of H100 dense FP16 peak, close to the fraction the A100 kernel reached. Whether the fraction survives on Hopper, which needs WGMMA instructions to reach its ceiling, is untested.

Landmark results

The rows below have different timing scopes. They are evidence for mechanisms, not a cross-hardware ranking.

Limitations and open problems

A fast statistic can be a different statistic

The integer-shift benchmark avoids the hardest problem by construction. With fractional shifts, the Fourier method samples a band-limited interpolant of $\psi$, and the variance of an interpolated $\psi$ is a quadratic form in the covariance of neighbouring samples:

Interpolating $\phi$ with the same weights does not compute this, because squaring a band-limited PSF doubles its bandwidth. Near masks and with sharp PSFs the error is not small. The kit measured Fourier null-variance ratios of 1.31, 1.17 and 1.04 at a mask boundary for PSF widths $\sigma=0.55$, 0.65 and 0.85 pixels, and the separable prototype found stack null-variance errors of $\pm4\%$ at $\sigma=0.65$. The exact fixes exist: a doubled-bandwidth information spectrum, or, for separable PSFs, modified axis tables that keep the same GEMM count. They have to be implemented and charged. This failure follows from the likelihood, not from GPU arithmetic, and it applies to any shortcut that interpolates a filtered map.

Trials and tails

For $M$ independent unit-normal scores the expected number of false excursions above $z$ is about $M\,\overline\Phi(z)$, where $\overline\Phi$ is the normal upper-tail probability. Heliostack discusses roughly $2\times10^{14}$ trials; at $7\sigma$ that is about 256 raw excursions under the Gaussian model (estimate), which is why a nominally overwhelming threshold still yields a candidate list. Small calibration errors matter enormously in this tail. If the true noise standard deviation is 1.01 times the assumed one, the single-trial false-alarm probability at a nominal $8\sigma$ rises by about 1.9 times (estimate). A precision or interpolation change that looks harmless on a mean SNR plot can therefore change the survey’s false-alarm rate. Equal sensitivity has to be shown as completeness at matched false alarms on identical injections.

Amdahl’s law

If a kernel takes a fraction $f$ of the pipeline and becomes $s$ times faster, the whole pipeline speeds up by

A tenfold kernel gain gives 1.82 times overall at $f=0.5$ and 5.26 times at $f=0.9$. The case study already shows the effect inside one scan: post-processing was 3% of the first fused kernel’s time and 17% of the best one’s. In real pipelines the unchanged stages are larger: image loading, PSF filtering, candidate filtering, stamps, vetting. DEEP’s 576 s per CCD includes all of them. The October 2026 Allen Telescope Array prototype is a clear example in radio. Of its 98 ms integrations it spends 33 ms on GPU imaging and 20 ms on CPU dispersion combination, but writing the dedispersed output does not fit in one integration and candidate ranking is offline (Perez et al. 2026, §§6-7).

Representations that do not fit

Dense reformulations fail through padding and materialization. A one-hot trajectory selector written as a dense matrix is mostly zeros. A Fourier phase matrix pays only if it is reused enough times. An FFT helps only if the search family has translation or separability structure. Arbitrary long-arc Keplerian motion, irregular cadence, varying WCS, chip gaps and changing PSFs all weaken those symmetries. Exposure trailing breaks exact axis separability: at 5″/h a 120 s DECam exposure trails by 0.64 pixels. A gather also supports early rejection, abandoning a trajectory as soon as partial sums make detection impossible, while a GEMM computes everything.

Where the room is

The physical limits are photon statistics, PSF bandwidth, sampling, cadence and the degeneracies of short arcs. No GPU infers motion parameters the observations do not constrain. Everything else in this chapter is an engineering or algorithmic limit: HBM bandwidth and capacity, coordinate arithmetic, padding, transposition, launch overhead, candidate transfer. Those can be reduced. The size of the template bank lies between the two: physics sets a minimum coverage, and placement, hierarchy and reuse decide how much redundant work is done (Chapter 38).

Build the fair gather first. The cheapest large gains are in the conventional kernel. The A100 measurements show a gather at 135-145 billion samples per second, held back by FP64 coordinates, against a streaming bound near 500 billion. FP32 coordinates from precomputed per-exposure offsets, velocity packets staged through shared memory, decode-once compressed storage, fused thresholding and compact candidate output are known techniques, and together they could plausibly give two to three times on the same card (estimate). Every claimed tensor-core or algorithmic gain should be quoted against that kernel, in paired samples per second as well as FLOP/s, with measured HBM traffic.

Finish the tensor-core kernel, then stop expecting miracles from it. At 9% of peak, the fused contraction has headroom. If it reached the Tensor-Core Beamformer’s 55%, kernel time on this workload would fall from 0.225 s to about 0.04 s, and with post-processing and overhead unchanged the scan would take about 0.11 s: roughly 16 times the kit gather and roughly 4.5 times the ideal no-reuse gather (estimate). Post-processing would then dominate, so the inverse FFT, crop and scoring must be fused too. That is a real but bounded prize, confined to dense long stacks: Rubin deep drilling fields, HSC or New Horizons multi-night sequences, and Hubble or JWST stares with hundreds to thousands of fully overlapping exposures.

Make many hypotheses share a representation. Order-of-magnitude gains at equal sensitivity require reducing the number of operations, not executing them faster. Signed partial sums over nightly segments, screened through coarse orbit cells and refined exactly on the original pixels, are the route the project’s frontier study ranked first for dense multi-night fields. Its CPU tests on real New Horizons/Subaru timestamps were encouraging but are not GPU measurements (Chapter 21, Chapter 42). Local derivative or low-rank bases within an orbit cell, nonuniform FFTs for irregular cadence (Chapter 39), and hierarchical trees in the spirit of FDMT (Zackay & Ofek 2017) are the other candidates. Each must keep the statistic signed, compute its normalization exactly under masks, include the cost of building the representation, and fall back where the representation fails.

Spend precision deliberately. Noise-scaled INT8 storage of $(\psi,\phi)$ with FP32 accumulation is nearly free in SNR under the Gaussian model. It matters most where it changes capacity, letting a thousand-image cube fit in HBM, rather than where it changes bandwidth for a kernel that is not bandwidth-bound. Use per-image or per-tile scales, an exception path for bright residuals, and a validated invalid code.

Benchmark the experiment, not the kernel. A fair acceptance plot puts runtime or energy on one axis and completeness at matched false alarms on the other, with identical injected sources, masks and bank. It separates four different things: accelerating the exact bank, coarsening the bank, shrinking the physical search volume, and approximate screening. All four are useful, but they answer different scientific questions. The biggest evidence gap in this field is still a public, sensitivity-controlled benchmark of optical gather kernels, compressed and tiled variants, and genuine tensor reformulations on the same H100 or B200. This chapter’s A100 measurement is a first, partial entry in that table, not the table.

References

1. Bassa, C. G., Romein, J. W., Veenboer, B., van der Vlugt, S., Wijnholds, S. J. (2021). Fourier-domain dedispersion. Astronomy & Astrophysics. arXiv:2110.03482
2. Barr, E. D. et al. (2026). The EDD radio astronomy backend framework. arXiv preprint. arXiv:2609.11177
3. Bernstein, G. M., Bebek, C., Rhodes, J. et al. (2010). Noise and bias in square-root compression schemes. PASP. arXiv:0910.4571
4. Clark, M. A., La Plante, P. C., Greenhill, L. J. (2013). Accelerating radio astronomy cross-correlation with graphics processing units. International Journal of High Performance Computing Applications. arXiv:1107.4264
5. Dimoudi, S., Adámek, K., Thiagaraj, P. et al. (2018). A GPU implementation of the correlation technique for real-time Fourier domain pulsar acceleration searches. ApJS. arXiv:1804.05335
6. DIRAC Institute / KBMOD contributors (accessed 2026). KBMOD search and PSF convolution source, revision f4b16e4. GitHub. github.com/dirac-institute/kbmod
7. Fu, D. Y., Kumbong, H., Nguyen, E., Ré, C. (2023). FlashFFTConv: efficient convolutions for long sequences with tensor cores. arXiv preprint. arXiv:2311.05908
8. Napier, K. J., Holman, M. J., Lin, H.-W., Gerdes, D. W., Ruch, T. R. (2026). heliostack: a novel approach to minor planet discovery. Planetary Science Journal (accepted). arXiv:2604.16731
9. Novotný, J., Adámek, K., Clark, M. A., Giles, M., Armour, W. (2023). Accelerating dedispersion using many-core architectures. arXiv preprint. arXiv:2311.05341
10. NVIDIA (2020). NVIDIA A100 Tensor Core GPU architecture (whitepaper). nvidia.com
11. NVIDIA (accessed 2026). A100 Tensor Core GPU specifications. nvidia.com/a100
12. NVIDIA (accessed 2026). H100 Tensor Core GPU specifications. nvidia.com/h100
13. NVIDIA (accessed 2026). HGX B200 and DGX B200 specifications. nvidia.com/hgx, nvidia.com/dgx-b200
14. NVIDIA (accessed 2026). CUDA C++ Best Practices Guide. docs.nvidia.com
15. NVIDIA (accessed 2026). cuFFT Library User Guide. docs.nvidia.com
16. Oostrum, L., Veenboer, B., Rook, R., Brown, M., Kruizinga, P., Romein, J. W. (2025). The Tensor-Core Beamformer: a high-speed signal-processing library for multidisciplinary use. arXiv preprint. arXiv:2505.03269
17. Pence, W. D., White, R. L., Seaman, R. (2010). Optimal compression of floating-point astronomical images without significant loss of information. PASP. arXiv:1007.1179
18. Perez, K. I., Bright, J. S., Gajjar, V., Cruz, L. F. (2026). Towards a GPU-accelerated real-time imaging pipeline for SETI at the Allen Telescope Array. arXiv preprint. arXiv:2610.06410
19. Reed, I. S., Gagliardi, R. M., Shao, H. M. (1983). Application of three-dimensional filtering to moving target detection. IEEE Transactions on Aerospace and Electronic Systems. doi:10.1109/TAES.1983.309401
20. Romein, J. W. (2020). The Tensor-Core Correlator (ALMA workshop presentation). osf.io/vzrcw
21. Romein, J. W. (2021). The Tensor-Core Correlator. Astronomy & Astrophysics 656, A52. doi:10.1051/0004-6361/202141896
22. Romein, J. W. and contributors (accessed 2026). Tensor-Core Correlator source and README. git.astron.nl
23. Smotherman, H. et al. (2023). The DECam Ecliptic Exploration Project (DEEP) VI: first multi-year observations of trans-Neptunian objects. arXiv preprint. arXiv:2310.03678
24. Tillet, P., Kung, H. T., Cox, D. (2019). Triton: an intermediate language and compiler for tiled neural network computations. Proceedings of MAPL 2019. doi:10.1145/3315508.3329973
25. White, J. et al. (2022). Bits missing: finding exotic pulsars using bfloat16 on NVIDIA GPUs. arXiv preprint. arXiv:2206.12205
26. Whidden, P. J., Kalmbach, J. B., Connolly, A. J. et al. (2019). Fast algorithms for slow moving asteroids: constraints on the distribution of Kuiper Belt objects. Astronomical Journal. arXiv:1901.02492
27. Williams, S., Waterman, A., Patterson, D. (2008). Roofline: an insightful visual performance model for floating-point programs and multicore architectures. UC Berkeley Technical Report UCB/EECS-2008-134. EECS-2008-134
28. Zackay, B., Ofek, E. O. (2017). An accurate and efficient algorithm for detection of radio bursts with an unknown dispersion measure. ApJ. arXiv:1411.5373
29. Finding Planets project (2026). A100 benchmark of separable-phase tensor-core velocity stacking, NERSC Perlmutter, 16 October 2026 (benchmark kit and tc_stack.py). Unpublished measurement.