Fast transforms and optimal coadds

In December 2015 Barak Zackay and Eran Ofek, at the Weizmann Institute, posted two papers with a title that sounds like a beginner’s question: How to coadd images? Astronomers had been adding images for a century. The answer was not obvious all the same. Exposures taken in different seeing are not repeated measurements of the same picture: each is the sky blurred by a different point-spread function (PSF), and averaging them pixel by pixel, however carefully weighted, mixes sharp information with blurred. Zackay and Ofek showed that for a static sky in background-limited Gaussian noise there is one coadded image, with its own PSF and flux scale, that preserves the entire likelihood of the data. On real Palomar Transient Factory exposures, popular weighting schemes had been throwing away between a few percent and about 25% of the information (Zackay & Ofek 2017a; Zackay & Ofek 2017b).

For planet searches the result matters twice. A static coadd is the reference against which moving things are found (Chapter 10). And the same mathematics says exactly what a moving-object search must keep: not a coadd, which has already discarded the time at which each photon arrived, but per-exposure score and information maps that can be summed along any trial trajectory. That sum is the inner loop of KBMOD, heliostack and every shift-and-stack code in Part V. It is a huge number of weighted additions, and many of them are repeated between neighbouring trajectories.

The second half of this chapter asks how much of that repetition can be removed. Mathematics offers many fast algorithms for structured sums: the fast Fourier transform (FFT) and its nonuniform relatives, recursive line and ray transforms, dedispersion trees, butterfly factorisations, fast multipole methods and random sketches. Each is fast because the operator has a particular structure. Each fails when the structure is absent, and most of them can quietly change the noise normalisation of the statistic they compute. Understanding which structure a planet search actually has is the point of the exercise.

One sum, many hypotheses

Let exposure $i$ contain background-subtracted pixels $d_{ip}$ with independent noise variance $\sigma_{ip}^2$ and pixel-integrated PSF $P_i$. A source of constant reference flux $f$, with known relative throughput or light-curve factor $a_i$, following trajectory $x_\theta(t)$, predicts

Here $\theta$ collects the initial position and whatever velocity, distance or acceleration parameters the baseline resolves. For Gaussian background-limited noise the log likelihood ratio is quadratic in $f$, and everything it needs is in two additive sums:

The flux estimate is $\hat f=A/B$ with variance $1/B$. The signed detection score is $\nu=A/\sqrt B$, a unit normal under noise for any fixed $\theta$, and the maximised log likelihood improvement is $A^2/2B$, or $\max(A,0)^2/2B$ if the flux must be positive. This is the KBMOD likelihood construction (Whidden et al. 2019). Masks are zero inverse variance. Correlated noise replaces the per-pixel weights with the inverse covariance.

Two rules follow, and most bugs in stacking codes break one of them. First, combine $A$ and $B$, then normalise. Adding per-image significances, or dividing a sum by a fixed $\sqrt N$, is wrong whenever seeing, transparency, variance or coverage changes between exposures. Second, the information map $\phi_i$ contains the square of the PSF correlated with the inverse variance. It is not the PSF applied to the variance map, and it is not the square of the filtered variance.

Write $E$ for the number of exposures, $P$ for the number of starting positions and $V$ for the number of motion hypotheses, so there are $M=PV$ trials. Once the maps are built, direct evaluation costs $O(EM)$ reads and additions. KBMOD’s original GPU kernel ran at $2.4\times10^{-11}$ seconds per image per trajectory, about 240 billion trajectories through 30 images in 180 seconds (Whidden et al. 2019). That is the baseline. Any fast transform has to beat this practical direct sum, not a deliberately naive dense matrix product. It also has to write or scan at least $M$ outputs if it produces every score.

Coadding a static sky

The measurement model

Write the background-subtracted image from exposure $j$ as

where $T$ is the static sky, $P_j$ the unit-integral effective PSF (including detector pixel response), $F_j$ the flux calibration (exposure time, throughput, area) and $\epsilon_j$ noise. The ideal assumptions are that the noise is Gaussian, independent between pixels and exposures, with known and locally constant variance $\sigma_j^2$, that astrometry and PSFs are known, and that the scene does not change (Zackay & Ofek 2017a). Source photon noise is neglected, which is right for faint objects and wrong near bright stars.

The general likelihood coadd

For an arbitrary linear model $d_j=A_jt+n_j$, where $A_j$ includes calibration, PSF, geometry, pixel integration and the selection of valid pixels, and $C_j$ is the noise covariance, the log likelihood of the sky vector $t$ is

The pair $(b,H)$ is sufficient: every likelihood ratio about $t$ depends on the data only through it. $H$ is the Fisher information matrix and also the normal-equation operator. Its off-diagonal elements describe how neighbouring sky pixels respond together, and they are not optional. James Bosch’s Rubin technical note calls this the likelihood coadd, and estimates that its nonzero entries for one $4{\rm k}\times4{\rm k}$ patch would take about 200 GB in single precision, a modelling estimate rather than a measured requirement (Bosch 2016). Storing $H$ is the hard part. Applying it is easy: forward-project a trial sky, weight the valid residual pixels, and apply the adjoint, which supports matrix-free iterative methods with any mask.

Point sources and the limiting magnitude

For a point source of flux $f$ at position $q$, let $a_{jq}$ be its predicted unit-flux pixel vector in exposure $j$. The log likelihood ratio is $f\Psi(q)-\tfrac12f^2\Phi(q)$ with

so $\hat f=\Psi/\Phi$, $\nu=\Psi/\sqrt\Phi$ and $2\Delta\ell_{\max}=\Psi^2/\Phi$, as above. For a circular Gaussian PSF of standard deviation $s_j$, well sampled, $\|P_j\|_2^2=1/(4\pi s_j^2)$, so

This is the physically useful measure of an exposure’s worth: transparency, background and PSF area all enter, and seeing enters squared. For $N$ identical exposures the SNR grows as $\sqrt N$ and the limiting magnitude improves by $1.25\log_{10}N$, about 2.5 mag for 100 exposures. That assumes a correctly aligned constant source and a fixed threshold. Searching more trajectories can force the threshold up (Chapter 38).

Scalar weights

The simplest coadd converts exposures to common flux units $I_j=M_j/F_j$ and averages them with inverse-variance weights $w_j=F_j^2/\sigma_j^2$. If the weights are constant over an object, its effective PSF is the weighted mean of the input PSFs. This is the right estimate of a common pixel signal when all PSFs are the same. When they differ, the inputs are differently blurred signals, and no single scalar weight per exposure can treat every spatial frequency optimally: high frequencies are present only in the sharp exposures, low frequencies in all of them.

A common improvement favours sharp exposures. In the raw-image weighted-sum convention, Zackay and Ofek list the Annis weights as $w_j\propto F_j/(s_j^2\sigma_j^2)$ and the Jiang weights as $F_j/(s_j\sigma_j^2)$. Weight formulas cannot be compared without stating which units they multiply. Streaming a weighted mean costs $O(NP)$ for $N$ exposures of $P$ pixels, plus $O(NPk)$ for resampling with $k$ taps. One consequence of the exact likelihood is worth stating plainly: every correctly modelled exposure adds nonnegative information. Rejecting poor-seeing exposures can help an approximate pipeline, but under the exact model it never helps.

Matched-filter coaddition

For point-source detection the optimal operation is to correlate each exposure with its own PSF and add:

then divide by the square root of its null variance to get $\nu$ (Zackay & Ofek 2017a). Filtering a scalar-weighted coadd once with a single approximate PSF is different and worse. For spatially varying noise, the same recipe computes per-exposure fields

for exposure j:
    M, F, P = calibrated image, throughput, local PSF
    ivar = inverse background variance, zero where masked
    Psi[j] = F    * correlate(M * ivar, P)
    Phi[j] = F**2 * correlate(ivar, P**2)
Psi_total = sum over j of Psi[j];  Phi_total = sum over j of Phi[j]
nu = Psi_total / sqrt(Phi_total)  where Phi_total > 0

Filtering costs $O(NPk)$ for a PSF of $k$ pixels or $O(NP\log P)$ by FFT, and Gaussian PSFs are separable into two one-dimensional passes. Dustin Lang and David Hogg report that Gaussian approximations retained above 97% matched-filter efficiency for the DESI Legacy Imaging Surveys PSFs, an empirical figure for those data (Lang & Hogg 2026). $S$ is an optimal detection score, not an image. Its pixels are correlated and its response to a point source is a sum of PSF autocorrelations, so it must not be fed to software that applies another PSF filter.

Proper coaddition

Paper II turns the detection score into an image (Zackay & Ofek 2017b). With hats for Fourier transforms, define

$A$ is the Fourier transform of the matched-filter sum $S$, and $B$ is its noise power: the information available at each spatial frequency. The proper image, its PSF and its flux scale are

Dividing by $\sqrt B$ whitens the noise, so

The sufficiency proof is short. In Fourier space the original log likelihood is a sum over frequencies, and the only data-dependent term involving the sky is $\mathrm{Re}(\hat T^*A)$, with quadratic coefficient $B$. Since $A=F_R\hat P_R\hat R$ and $B=F_R^2|\hat P_R|^2$, fitting $R$ with PSF $P_R$, flux scale $F_R$ and unit white noise reproduces the same likelihood up to a constant that does not depend on the sky. And applying the matched filter of $P_R$ to $R$ gives back $F_R\hat P_R^*\hat R=A$, the optimal detection score. $R$ is not a deconvolved sky. It is an image-shaped representation of the normal equations.

A = complex zeros on the Fourier grid;  B = real zeros
for exposure j:
    Mhat = FFT(registered background-subtracted image j)
    Phat = FFT(unit-integral PSF j, origin at pixel 0)
    A += F[j] * conj(Phat) * Mhat / sigma[j]**2
    B += F[j]**2 * abs(Phat)**2 / sigma[j]**2
Rhat = A / sqrt(B) on measured modes only
FR = sqrt(sum of F[j]**2 / sigma[j]**2);  PRhat = sqrt(B) / FR
R = inverse FFT(Rhat);  PR = inverse FFT(PRhat)

The cost is $O(NP\log P)$, streaming, with $O(P)$ working memory. The accumulators $(A,B)$ are additive, so proper coadds of disjoint subsets combine exactly. Modes where $B$ is zero carry no information and must be dropped, not filled; regularising them changes the likelihood. Padding prevents periodic FFT boundaries from mixing opposite image edges. The proper PSF often has a narrow core and can have negative wings, and a narrow FWHM is not itself a measure of information: $B(k)$ is.

The real test used PTF field 100031 near M51, with four input sets of 45, 48, 39 and 17 exposures and deep references built from 425 and 263 exposures. Measured as the median squared SNR of detected stars relative to the optimal score, equal weighting retained 0.53–0.83, the Annis weights 0.82–0.93, the Jiang weights 0.81–0.92 and proper coaddition 0.99–1. In set 3, proper coaddition carries $1/0.84\approx1.19$ times the information of the Annis coadd (derived from the table). At the nominal $4\sigma$ threshold that set gave 1,976 reference-matched and 230 unmatched detections with proper coaddition against 1,796 and 159 with Annis weights. More detections of both kinds do not by themselves mean a worse false-alarm trade-off: the two images have different resolution and correlated pixels (Zackay & Ofek 2017b).

The idea has antecedents. Zackay and Ofek acknowledge a related method attributed to a 2004 communication by Nick Kaiser, and Bosch’s note identifies Kaiser’s unpublished Pan-STARRS document Addition of images with varying seeing (Bosch 2016). The proper coadd’s distinctive contributions are the explicit sufficiency theorem and its white-noise image form.

What resampling does to noise

Every coadd starts by warping exposures onto a common grid, and warping correlates noise. A linear warp $W$ transforms data and covariance together: $d'=Wd$, $A'=WA$, $C'=WCW^T$. Even if $C$ is diagonal, $C'$ is not. For one-dimensional stationary interpolation with kernel $h$ and white input noise,

The variance plane stored with a resampled image keeps only the $\Delta=0$ term. For any filter or aperture $u$, the true variance is $u^TC'u$, which includes all the cross terms, not $\sum_pu_p^2C'_{pp}$.

Andrew Fruchter and Richard Hook’s Drizzle, developed for the Hubble Deep Field North, made this concrete (Fruchter & Hook 2002). Drizzle shrinks each input pixel to a “drop” of linear size pixfrac $p$, maps it through the geometric distortion onto a finer output grid of relative scale $s$, and distributes its flux by overlap area. It reconstructs undersampled, dithered images with little extra blur, and it correlates the output noise whenever one input pixel lands in several output pixels.

For square drops and uniform, dense dithers, the ratio by which large-aperture noise exceeds the variance-plane prediction is $\mathcal R=r/(1-1/3r)$ for $r=p/s\ge1$ and $1/(1-r/3)$ for $r\le1$. Fruchter and Hook’s example, $p=0.6$ and $s=0.5$, gives $\mathcal R=1.662$. Real dither patterns need their own calibration.

Dustin Lang’s unWISE coadds made a related point about resolution (Lang 2014). The public WISE Atlas images had been PSF-filtered for point-source detection, which made them about 8.5 arcseconds wide in the first three bands. Lang’s coadds, resampled with Lanczos kernels, kept about 6 arcseconds, better for forward-modelling crowded scenes. A typical tile combined 432 W1 exposures, about $3\times10^8$ input pixels, in 630 seconds on one CPU core with 2.3 GB of memory. The north-ecliptic-pole tile touched about $10^{10}$ input pixels from 20,000 exposures and took about 12 hours with 100 GB of memory. The unWISE images reappear in the infrared searches of Chapter 24.

Controlled reconstruction: Imcom and Fast Imcom

Barnaby Rowe, Christopher Hirata and Jason Rhodes’s Imcom (2011) solves a different problem: reconstruct an image with a requested output PSF from undersampled, irregularly dithered inputs, and report both the mismatch to that PSF and the output noise (Rowe, Hirata & Rhodes 2011). Each output pixel is a linear combination $h_\alpha=\sum_iT_{\alpha i}d_i$ of input pixels, with weights that minimise a sum of PSF leakage and noise variance, and the output covariance $TCT^T$ is computed exactly. The solve costs $O(n^3+n^2m)$ for $n$ input and $m$ output pixels per region. For a hypothetical space weak-lensing survey the authors estimated about $10^{21}$ operations, roughly 12 days on a sustained petaflop machine.

Kaili Cao’s Fast Imcom, submitted in October 2025 and revised in January 2026, changes the mathematics before accelerating it (Cao 2026). It precomputes continuous interpolation weight fields from oversampled PSFs, regrids each exposure independently, and then combines the regridded images with a small image-count-sized set of weights, replacing the large pixel-level solve. For the Roman High Latitude Imaging Survey its Table 1 contrasts about 100 million core-hours and 1.5 PB for the existing approach with a projected fewer than 10 million core-hours and about 0.2 PB. The paper calls the speed prediction educated speculation: its demonstrations are one-dimensional and the two-dimensional production code is future work.

Coadds in a survey pipeline

Production pipelines rarely use the proper Fourier formula, for practical reasons. Rubin’s first data preview (DP1) built per-band coadds by inverse-variance mean stacking of warped visits, on 2.8-square-degree tracts divided into $10\times10$ patches of about 128 square arcminutes. Deep-coadd inputs had to have PSF FWHM below 1.7 arcseconds, and difference-imaging templates used the best-seeing third of visits, or the best 12 when fewer than 36 were available (Rubin Observatory 2025). A direct mean lets the pipeline compute the coadd PSF lazily at any position from the contributing exposures’ PSFs and weights, though that PSF jumps at input edges.

Artifact rejection is the other reason. Sigma-clipping differently blurred images preferentially clips the cores of bright stars in the sharp exposures, so the effective PSF becomes flux dependent, and Rachel Mandelbaum and colleagues show that even mild clipping or median combination can break the coadd PSF (Mandelbaum et al. 2023). Rubin’s CompareWarp therefore separates the nonlinear decision from the linear mean. It compares PSF-matched warps against a static model, finds residual footprints, uses their temporal recurrence to tell transient artifacts from variable sources, masks them in the direct warps and then takes a plain weighted mean. In its example the stellar-locus width was 14.35 mmag against 18.48 mmag with sigma-clipping (AlSayyad 2018). Asteroids are deliberately removed from a static reference by exactly this logic, so a moving-object search must keep the visit-level data.

The 2026 pipeline paper for the Large Array Survey Telescope (LAST) is a reminder that proper subtraction and proper coaddition are separate choices. LAST uses the ZOGY subtraction algorithm, but its references were still weighted coadds, with proper coaddition planned (Konno et al. 2026).

The motion boundary

A static proper coadd is sufficient for a static sky. It is not sufficient for an unknown trajectory, because it has summed away the epoch labels.

For trajectory $q_j(\theta)$ the moving-frame statistic is

The expensive filtering happens once per exposure. The search only reads the two fields. A trajectory through fewer valid epochs automatically gets less information, and setting a masked numerator to zero while keeping the unmasked denominator is a bug. If brightness follows a known pattern $f_j=fa_j$, the sums carry $a_j$ and $a_j^2$. An unknown rotational light curve can be expanded in a small temporal basis, giving a score vector and a small information matrix instead of scalars (derived extension).

Hayden Smotherman and colleagues’ 2021 KBMOD survey shows the practicalities (Smotherman et al. 2021). A reference image that contains the moving source subtracts part of it. They split their observations into groups at least twelve hours apart and differenced each against the opposite group’s reference, so that objects moving faster than their minimum searched rate were never at their reference positions. Their search over about 132 square degrees reported 75 unlinked detections and 24 recoveries of known objects, the faintest at $VR=25.02\pm0.93$. They also quote the KBMOD kernel’s capability of $10^{10}$ trajectories through 10–15 $4{\rm k}\times4{\rm k}$ images in under a minute on a GeForce GTX 1080 Ti. Difference imaging reduced the number of individually detected sources but raised high-score trajectory contamination, which needed further filtering.

The per-epoch maps are large. In single precision a score and information cube costs $8NP$ bytes, about 13.4 GB for $N=100$ exposures of $4096^2$ pixels (estimate, excluding masks). That can exceed a GPU’s memory and make host-to-device transfers the bottleneck, which is why tiling and compact encodings matter (Chapter 40).

Approximation is a template change

Every fast algorithm below computes, exactly or approximately, some modified sum. The right way to judge it is to ask what template it actually applies. Whiten the data, $y=C^{-1/2}d$, and the source template, $h_\theta=C^{-1/2}s_\theta$. If an algorithm effectively applies $\tilde h_\theta$ instead of $h_\theta$, the fraction of the ideal expected SNR it retains is

This is the cosine of an angle in noise-weighted template space. An algorithm can reproduce a template’s total flux and still get its noise norm wrong. Smith’s practical audit is to apply the algorithm’s adjoint to a unit impulse at one output, recover the effective input template, and measure its overlap and norm (Smith 2016).

For a well-sampled Gaussian PSF of width $\sigma_i$, a positional error $\delta x_i$ in exposure $i$ gives $R\simeq\sum_ic_i\exp[-\|\delta x_i\|^2/4\sigma_i^2]$, with $c_i$ the information weights (Geringer-Sameth et al. 2025). A common displacement costing 2% of the SNR is $0.284\sigma$, or $0.121$ FWHM (derived). That is the scale to which any approximate path must be accurate.

Interpolation needs particular care. If a fast algorithm interpolates score-map values with weights $w$, the variance of the interpolated score is $w^T\mathrm{Cov}(\psi)w$, including the covariance between neighbouring score pixels. Interpolating the variance map with the same weights gives a different, wrong answer. Near masks and changing PSFs, a fast numerator with a convenient denominator produces a miscalibrated tail. The project’s design studies (unpublished, 2026) found this the most common failure across many proposed accelerations: interpolated or Fourier-shifted statistics under masks and variable seeing distorted the null variance by up to tens of percent. The arithmetic shows why it matters. If the noise width is underestimated by 1%, the false-alarm rate at a $7.65\sigma$ threshold rises by a factor of about 1.8. The exact correction in those studies cost about 2.5 times the numerator work, and it has to be booked.

A single template’s false-alarm probability is not the search’s. The union bound $P(\max_\theta\nu_\theta>u)\le M\bar\Phi(u)$ holds for correlated templates, but a transform that changes the bank or its correlations needs its own empirical calibration of the search maximum. Lang and Hogg’s illustrative calculation shows how fragile the tail is: a residual background of only $0.05\sigma$ per pixel can double the false-positive rate at a $5\sigma$ threshold in good seeing, because coherent errors survive averaging much better than independent noise (Lang & Hogg 2026).

Fast transforms

Fourier identities and the NUFFT

For straight-line motion $x_\theta(t_i)=x_0+vt_i$, the Fourier shift theorem turns a moving sum into phases. With $\hat\psi_i(k)$ the spatial Fourier transform of exposure $i$‘s score map,

At each spatial frequency $k$, the sum over exposures is an exponential sum with irregular nodes $t_i$ and target frequencies $k\cdot v$. That is a nonuniform FFT (NUFFT) problem of type 3, or type 1 if the target frequencies can be put on a regular grid. The identities are a derived formulation of optical trajectory stacking. The NUFFT algorithms are established. The three types are

Type 2 is the adjoint of type 1, not its inverse. Inverting irregular samples is a separate linear inverse problem, and confusing the two is a real hazard when borrowing interferometric imaging code for detection.

A type-1 transform sorts the nonuniform points into bins, spreads each onto an oversampled regular grid with a compact smooth kernel, takes an FFT, and divides out the kernel’s transform. In $d$ dimensions with $J$ points, $N$ modes and kernel width $w=O(\log1/\epsilon)$ for accuracy $\epsilon$, the cost is $O(Jw^d+N\log N)$. Alex Barnett, Jeremy Magland and Ludvig af Klinteberg’s FINUFFT uses an “exponential of semicircle” kernel with cache-aware sorting (Barnett, Magland & af Klinteberg 2019), and cuFINUFFT carries the design to CUDA (Shih et al. 2021).

On a Tesla V100, a three-dimensional type-1 transform of $1.34\times10^8$ random points onto $256^3$ modes took 0.18 s at tolerance $10^{-2}$ and 0.87 s at $10^{-5}$ with shared-memory spreading, 16.1 and 3.9 times faster than FINUFFT on 28 CPU cores, using about 6.1 GB. Spreading took 92–98% of the execution time. The bottleneck of a GPU NUFFT is a scattered accumulation, not the FFT.

For a planet search three costs remain after the temporal sums. Producing a position image for every velocity costs $O(VP\log P)$ inverse FFTs. The type-3 transform’s cost depends on the space–bandwidth product, here proportional to the time span times the range of $k\cdot v$, so long gaps inflate it without adding photons; the FINUFFT documentation warns that type 3 can be slow even for small inputs when that product is large. And curved motion needs more dimensions: a local model $x(t)=x_0+\sum_rb_rg_r(t)$ lifts the sum into an $R$-dimensional type-3 transform, with spreading cost growing as $w^R$ (derived extension). The transform does not make extra motion dimensions free.

Radio interferometry is where astronomy already runs production NUFFTs. Philipp Arras and colleagues’ wide-field operator combines gridding, FFTs and the non-coplanar $w$ term with errors down to about $10^{-12}$, and is used in WSClean and the DUCC library (Arras et al. 2021).

The Fourier-slice transform

In a single image, the Fourier transform of a projection along parallel lines is a radial slice through the image’s two-dimensional Fourier transform. Paul Hickson used this to integrate long, faint orbital-debris streaks: FFT the image, interpolate radial slices, inverse transform each, and apply a transverse matched filter (Hickson 2018). He reports about 3, 10 and 30 seconds for $1024^2$, $2048^2$ and $4096^2$ images on a 2.2-GHz processor.

Interpolation in the Fourier plane is not cosmetic. Linear interpolation gave flux errors up to about 15%, a seven-point sinc interpolator under 5%, and about 50 points under 1% at four times the runtime. FFT periodicity also produced displaced ghost detections. These are photometric errors in one implementation, not universal SNR losses, but they show that a transform’s accuracy is a parameter with a cost. Streak detection is treated fully in Chapter 11.

Recursive line transforms

A different family avoids Fourier space altogether. The fast discrete Radon transform of Götz and Druckmüller (1996) and Brady (1998), explained clearly by William Press in 2006, sums an image along digital lines by recursion (Press 2006). Let $T_m(h,s)$ be the sum along a digital line of length $2^m$ starting at offset $h$ with total displacement $s$. A line twice as long is two shorter lines joined end to end:

For an $n\times n$ image there are $O(n^2)$ states per level and $\log_2n$ levels, so the whole transform costs $O(n^2\log n)$ instead of $O(n^3)$ for direct summation, with four orientations covering all slopes. The sum is exact along its own digital paths, not along straight lines. Press bounds the deviation by $(\log_2n)/6$ pixels in his convention, which grows with image size at fixed pixel scale. Kazimirov and Nikolaev’s 2025 FHT2SP algorithm bounds the error directly, at $O(wh\ln^3w)$ arithmetic for a $w\times h$ image, but its example shows the price: for a $1024^2$ image and a 5/8-pixel error bound it uses $29\times29$ superpixels, an effective $29696^2$ image (Kazimirov & Nikolaev 2025).

Guy Nir, Barak Zackay and Eran Ofek built a PSF-filtered streak detector on this recursion (Nir, Zackay & Ofek 2018). They observed satellite streaks on 20 July and 8 August 2017 with a 40-cm Meade telescope and an Andor Zyla 5.5 sCMOS camera at the Kraar Observatory, and detected low-Earth-orbit and GPS satellites. In simulations with a Gaussian PSF of $\sigma=2$ pixels the method recovered on average more than 95% of the theoretical SNR for long streaks. Their extreme-value fit to real images suggested a threshold near $\nu\simeq10.25$ for about one false positive per $10^6$ images, a fitted extrapolation rather than a million-image test.

The space–time cube: X-ray transforms and FaXT

Moving point sources across many exposures are lines through a space–time cube, not lines in an image. In three dimensions the Radon transform integrates over planes. The transform that integrates over lines is the John or X-ray transform, and Marichal-Hernández and colleagues give recursions for both (Marichal-Hernández et al. 2020). Their ray transform returns $n\times n\times2n\times2n$ values per orientation sector: an $n^3$ input has an $O(n^4)$ output, one value per starting position and velocity.

FaXT, by Nguyen, Woods, Ruprecht and Birge, applies this recursion to moving-object searches (Nguyen et al. 2024). When the number of exposures, the image side and each side of the velocity grid all grow as $n$, direct shift-and-add costs $O(n^5)$ and FaXT $O(n^4)$, proportional to the output size. Its abstract reports more than a thousand TESS asteroids below the single-frame detection limit; it does not establish how many were new. Matthews and Nguyen’s 2026 assessment found about a 100-fold speed-up over nearest-pixel summation at 256 frames, but a mean photometric loss rising from 32% at 64 frames to 37% at 256, against about 27–28% for nearest-pixel summation, for a simulated PRF with 0.5-pixel standard deviation (Matthews & Nguyen 2026). Some simulated objects were lost before refinement could rescue them. The digital-path geometry is drawn in Chapter 21.

Endpoint recursions: FDMT, Smith’s tree and ESAM

The Radon recursion works because a line is fixed by its endpoints. Barak Zackay and Eran Ofek’s fast dispersion measure transform (FDMT) showed that the same trick works for curves (Zackay & Ofek 2017c). A radio burst arrives at frequency $f$ at time $t(f)=t_0+K\,\mathrm{DM}\,f^{-2}$. Over any frequency interval the curve is fixed by its arrival time at one edge and its total delay, so a sum over a band is the sum over two sub-bands joined at the midpoint crossing:

The cost is bounded by $2N_tN_f+N_tN_\Delta\log_2N_f$ for $N_t$ time samples, $N_f$ channels and $N_\Delta$ delays, against $N_tN_fN_\Delta$ directly (Chapter 36). The paper states the condition for other curve families: a unique curve must pass through any two endpoints. A path $x(t)=x_0+b\,g(t)$ with known $g$ satisfies it. A free quadratic $x_0+vt+\tfrac12at^2$ does not, because two endpoints leave the curvature undetermined. A projected Keplerian orbit has even more state. Carrying extra derivatives through the recursion can close it, but the extra states must be counted.

Kendrick Smith’s pulsar acceleration tree shows how to do the counting (Smith 2016). Writing the quadratic phase in a rescaled Legendre basis around each interval’s midpoint, halving the interval coarsens the required acceleration and frequency resolution, so the number of intermediate states forms a convergent geometric series and the full coherent search costs $O(1)$ per final model. The default implementation was more than 90% optimal in the overlap sense.

ESAM, by Gupta, Bannister, Flynn and James, removes the need for an analytic curve law (Gupta et al. 2024). It takes an arbitrary bank of two-dimensional masks, splits each recursively into halves, and gives identical sub-products the same identifier, so evaluation becomes a bottom-up traversal of a graph of unique partial sums. It is common-subexpression elimination for a template bank, and its cost is proportional to the number of unique products, which approaches the brute-force cost when masks share nothing. Its example reproduced a brute-force bank’s SNR with about 10 times fewer operations. The authors did not provide a hardware timing and noted that brute force on a GPU might outrun their Python code.

Butterfly factorisations

Some operators are full rank globally but low rank in pairs of complementary boxes. Write an oscillatory transform as $u(x)=\sum_\xi a(x,\xi)e^{2\pi i\Phi(x,\xi)}g(\xi)$. Emmanuel Candès, Laurent Demanet and Lexing Ying showed that for a small output box $A$ and a large input box $B$ with widths $w_Aw_B\lesssim1/n$, the residual phase

varies slowly enough to interpolate with a few terms (Candès, Demanet & Ying 2009). The butterfly algorithm walks down one tree and up the other, transferring short local expansions, and evaluates an $n\times n$ to $n\times n$ Fourier integral operator in $O(n^2\log n)$ instead of $O(n^4)$, with separation rank $O(\log^4(1/\epsilon))$ under its assumptions. Those assumptions include phase homogeneity and careful treatment of the frequency origin; “the orbit is smooth” is not a substitute. Li, Yang, Martin, Ho and Ying’s butterfly factorisation writes such an operator as $O(\log N)$ sparse factors applied in $O(N\log N)$, but building them costs about $O(N^{3/2})$ (Li et al. 2015). A one-pass sky search may not amortise that.

Curved-trajectory sums with butterflies already exist outside astronomy. Hu, Fomel, Demanet and Ying accelerated the hyperbolic seismic Radon transform $Rd(\tau,p)=\int d(\sqrt{\tau^2+p^2h^2},h)\,dh$ this way (Hu et al. 2013), and Demanet and colleagues did the same for synthetic-aperture radar imaging (Demanet et al. 2012). Kielstra and colleagues’ tensor butterfly reaches $O(n^d)$ construction and application for $2d$-mode tensors at fixed rank, but its storage carries a factor $r_t^{2d}$ in the rank, which can dominate (Kielstra et al. 2024).

For an optical search, the relevant kernel is $K(\theta;(i,k))=e^{ik\cdot x_\theta(t_i)}$, mapping exposure–frequency data to a high-dimensional orbit bank. Its complementary ranks have to be measured on real cadences, ephemerides and PSF bandwidths before any factorisation is built. The project’s design studies (unpublished, 2026) measured Fourier ranks of 27–82 at $10^{-3}$ tolerance on local patches, with larger patches reaching 159, and modelled stored-factor implementations at only 0.1–2.5 times a direct gather even with 1,118 fully contributing epochs. Their conclusion was that a stored-factor butterfly is not a route to tenfold, and that an analytically generated, nested version is the missing algorithm.

Multipoles and hierarchical matrices

Leslie Greengard and Vladimir Rokhlin’s 1987 fast multipole method (FMM) replaced $O(N^2)$ pairwise interactions by multipole expansions of distant clusters (Beatson & Greengard). Hierarchical ($\mathcal H$ and $\mathcal H^2$) matrices express the same idea as algebra on block trees, with near-linear storage and application for bounded-rank blocks. These methods need kernels that are smooth and low rank between well-separated regions, like the Newtonian potential.

The optical matched filter is the opposite: a sharply concentrated tube around each trajectory, with essentially zero weight on distant pixels. A literal FMM does not remove the sum over exposure-dependent positions. Its plausible uses here are narrower: smooth interpolation, kernel density accumulation in catalogue or orbit coordinates (Chapter 21), or compression of structured background covariance. For a nonoscillatory kernel the FMM structure applies; for an oscillatory phase, the butterfly does.

Sparse FFTs and random sketches

Haitham Hassanieh, Piotr Indyk, Dina Katabi and Eric Price’s sparse FFT finds the $k$ largest Fourier coefficients of a length-$n$ signal in $O(k\log n)$ time when the spectrum is exactly $k$-sparse, and beat FFTW for $n=2^{22}$ and $k\le2^{17}$ in a preliminary implementation (Hassanieh et al. 2012). Its guarantee is an $\ell_2$ approximation relative to the spectrum’s tail. A whitened astronomical image’s noise is dense in every orthonormal basis, so the observed spectrum is never sparse, and an error bounded by the total tail can exceed a near-threshold source. A moving point’s Fourier signal is also a phase-organised ridge, not a few isolated bins. Sparse FFTs suit strong periodic lines and interference, not faint-object detection.

Random projections fail for a related, sharper reason, which a short derivation shows. Let the whitened data be $y=fh+n$ with $n\sim\mathcal N(0,I_D)$, and keep only $Sy$ for an $m\times D$ matrix $S$ with $m<D$. With the correct compressed covariance $SS^T$, the Fisher information about $f$ is

where $P_S$ projects onto the row space of $S$. For a randomly oriented row space this keeps on average a fraction $m/D$ of the original information $\|h\|^2$. Keeping one tenth of the dimensions leaves a typical SNR ratio of $\sqrt{0.1}\approx0.32$ (derived). The Johnson–Lindenstrauss lemma, which preserves the norms of a finite set of vectors with $O(\epsilon^{-2}\log(M/\delta))$ dimensions (Woodruff 2014), does not contradict this: it preserves large norms, not a small signal buried in noise. A sketch can preserve a low-dimensional template subspace that it was designed to contain, but a generic random projection cannot preserve every direction. Randomised methods are best used to find structure in an operator, as in Nathan Halko, Per-Gunnar Martinsson and Joel Tropp’s range finder $Y=K\Omega$, $K\approx QQ^*K$ (Halko, Martinsson & Tropp 2011), for example to discover the local template basis of an orbit cell. Sketch-based voting can also generate proposals, provided exact rescoring follows and the miss probability before rescoring is measured.

Landmark results

Tensor cores and the separable velocity sum

One reformulation gets particular attention in the project’s own design studies, because it maps the moving-object sum onto the dense matrix multiplication that GPU tensor cores are built for. On a Cartesian grid of velocities $(v_{x,a},v_{y,b})$, the Fourier-domain temporal sum factorises:

with $c_{ik}=a_i\hat\psi_i(k)$, $L_k[a,i]=e^{ik_xv_{x,a}t_i}$ and $R_k[b,i]=e^{ik_yv_{y,b}t_i}$. The transpose is not conjugated. For each spatial frequency, the sums over all velocities become one complex matrix product of size $V_x\times E$ by $E\times V_y$, and the phase tables grow as $V_x+V_y$ rather than $V_xV_y$. Irregular exposure times are allowed. The arithmetic is still $O(EV)$ per frequency: the gain is in data reuse and sustained matrix throughput, not in operation count. It competes with the NUFFT for the same temporal sum, so their gains do not multiply.

The design studies are unpublished models and CPU tests, not GPU measurements, and should be read as such. On CPU, FP16 phase tables retained 0.999999976 of the normalised SNR, an eight-pixel halo test retained at least 0.99957, and a separable construction of the squared-PSF information norm agreed with the exact masked norm to $1.7\times10^{-15}$ relative in a periodic test. Ordinary shifted information maps were generally not sufficient for that norm. Roofline models of an H100 projected 7.3–24 times the throughput of a strong direct-gather kernel at 1,000 fully overlapping epochs, and 2.7–10 times at 100, falling to 2.9–9.6 times at 1,000 epochs once a general mask-normalisation pass at 2.5 times the cost was charged. Small fused blocks underfilled the tensor cores and modelled at only about 0.4–1.6 times. At the footprint occupancy of the real heliostack search, about 51 epochs per output pixel, the advantage disappeared. Reaching tenfold in the 1,000-epoch case would need about 99 TFLOP/s of sustained complex contraction throughput for this matrix shape, which no measurement yet supports.

The same studies modelled several of the other structures in this chapter. Unmodified FaXT-style integer paths added up to 4.35% extra SNR loss at 256 frames with a PSF of $\sigma=1.2$ pixels and 8.35% at $\sigma=0.8$. A naive six-coordinate interpolating tree projected at only 0.55–0.94 times a direct gather at 1,024 epochs. These numbers are design-study models, and they make the same point from the other side: the fast operator must beat a strong, footprint-aware direct sum, not a weak one.

Limitations and open problems

Static sufficiency stops at the static model. A proper coadd preserves the likelihood of a constant scene with known PSFs and background-limited Gaussian noise. It removes the time coordinate needed for variability, proper motion, parallax and planets. Bright sources make the noise depend on the scene. Masks couple Fourier frequencies, so a zero-filled hole breaks the scalar Fourier formula, and stationary coloured noise needs a noise power spectrum in place of $\sigma_j^2$. Interpolating across artifacts repairs the picture but not the statistics. A perfect compression of the wrong forward model preserves the wrong likelihood, and proper coaddition is not self-calibration.

Shared references correlate epochs. If a common noisy reference $r$ is subtracted from several exposures, $d_i=s_i-K_ir$, the residuals are correlated across epochs: $\mathrm{Cov}(d_i,d_j)=\delta_{ij}C_{s_i}+K_iC_rK_j^T$. Whitening each difference image separately does not make a trajectory’s samples independent. Deep references, a joint covariance or a partition of the exposures that controls reference reuse each trade sensitivity against computation.

Digital exactness is not physical accuracy. Recursive transforms are exact along their own paths. Broader smoothing tolerates path errors but changes the noise and the source confusion. FaXT’s photometric loss and FHT2SP’s expanded mesh are two faces of the same tension.

State closure. A tree works only if a parent trajectory restricts unambiguously to its children. Constant velocity closes. Free acceleration needs more state, and Keplerian motion with observer parallax needs a physically justified local representation. A recursion whose hidden state grows has not solved the high-dimensional bank problem.

Output dimensionality survives. A space–time FFT does not deliver an arbitrary position–velocity–acceleration search, and a plane transform does not deliver ray integrals. The output bank and its resolution survive most reformulations: the estimated $10^{19}$ trials for a six-year ZTF search beyond 100 au remain prohibitive even at constant work per trial (Geringer-Sameth et al. 2025).

Footprints and bandwidth. Dense transforms spend work on empty padding when exposures overlap only partly, irregular cadence inflates temporal bandwidth, and sharp PSFs need more spatial modes. A transform with a favourable dense operation count can lose to footprint-aware direct accumulation.

What has not been shown. No source examined for this chapter demonstrates a NUFFT, butterfly or multipole replacement for a curved optical orbit bank, at matched completeness and matched false-alarm rate, on a modern GPU. Fast transforms are common in astronomy. That particular combination is not yet.

Where the room is

Separate the physics from the engineering. For a correct Gaussian model, the information matrix $H$ determines what can be learned. No kernel can create information at spatial frequencies where every PSF has zero transfer, remove photon noise, or recover the time coordinate from a static coadd. Exact proper coaddition already saturates the static likelihood, so a “better” coadd under the same assumptions must either exploit an assumption its comparator mishandled or change the objective. The $\sqrt N$ gain is physics. Repeated arithmetic, padding, factor construction, poor parameterisation and uncoalesced memory traffic are engineering, and that is where the room is (Chapter 42).

Fix the normalisation first. The denominator $B$ is cheap to get wrong and expensive to get right under masks, interpolation and variable seeing. Any proposed acceleration should report its effective template’s overlap and noise norm by an adjoint test, and its null distribution at the search threshold, before reporting a speed. An exact, cheap information norm for realistic PSFs and masks would itself be a useful result, because every faster numerator needs one.

Reuse score and information, and move fewer bytes. For moving-point searches the per-exposure maps already avoid re-filtering. Where coverage and variance are uniform, the information map collapses to a scalar per exposure with sparse corrections at masks, removing half the reads. A 13-GB cube makes tiling, epoch blocking and compact encodings at least as important as arithmetic (Chapter 40).

Match the transform to the regime. Four directions are worth testing, in order of evidence.

1. Local trees in physical coordinates. Extend FaXT and ESAM to signed score and information states in local orbital coordinates, compile identical child products once, and measure how many unique states survive and how often each is reused. The risk is the compact representation of curved, irregularly sampled paths with an affordable noise norm.
2. NUFFT or separable GEMM inside well-conditioned dense patches. Where hundreds to a thousand exposures share a footprint and motion is locally linear, the mathematics is transparent and GPU libraries exist. Batch fixed-node transforms, reuse plans, and include inverse FFTs, padding, the denominator and exact refinement in the timing. A numerator microbenchmark is not a detector.
3. Butterflies as a rank experiment first. Sample blocks of the orbit-to-Fourier kernel, strip separable phases, and measure singular-value decay at masks, field edges, narrow PSFs and strong parallax. If ranks stay bounded, try analytic or tensor transfers. If they rise, stop.
4. Randomisation to build operators, not to compress photons. Randomised range finders can find local template bases cheaply. Generic random compression of the data loses information as $m/D$.

For coadds themselves, the open engineering problem is an economical representation of a spatially varying likelihood that keeps the covariance near defects and boundaries: a stationary convolution in the interior plus low-rank or sparse corrections at masks, applied matrix-free. Fast Imcom’s change of mathematics, from pixel-sized solves to image-sized combinations, is the right kind of idea, and its two-dimensional benchmark with real masks and end-to-end timing is the experiment that would settle it.

An order-of-magnitude saving at unchanged limiting magnitude would still matter a great deal. It would pay for more trajectories, more nights, richer nuisance models and better completeness calibration. “Equal sensitivity” should mean equal recovery and false alarms over the same physical hypothesis space, measured with the same injections through the same pipeline.

References

1. AlSayyad, Y. (2018). Coaddition artifact rejection and CompareWarp. Rubin Observatory DMTN-080. dmtn-080.lsst.io
2. Arras, P., Reinecke, M., Westermann, R., Enßlin, T. A. (2021). Efficient wide-field radio interferometry response. Astronomy & Astrophysics. arXiv:2010.10122
3. Barnett, A. H., Magland, J. F., af Klinteberg, L. (2019). A parallel non-uniform fast Fourier transform library based on an "exponential of semicircle" kernel. SIAM Journal on Scientific Computing 41, C479. arXiv:1808.06736
4. Beatson, R., Greengard, L. A short course on fast multipole methods. Lecture notes. math.nyu.edu
5. Bosch, J. F. (2016, revised 2025). Flavors of coadds. Rubin Observatory DMTN-015. dmtn-015.lsst.io
6. Candès, E., Demanet, L., Ying, L. (2009). A fast butterfly algorithm for the computation of Fourier integral operators. Multiscale Modeling & Simulation 7, 1727. arXiv:0809.0719
7. Cao, K. (2026). Linear image regridding and coaddition with oversampled point spread functions: lessons from 1D. Astronomical Journal 171, 140. arXiv:2510.16110
8. Demanet, L., Ferrara, M., Maxwell, N., Poulson, J., Ying, L. (2012). A butterfly algorithm for synthetic aperture radar imaging. SIAM Journal on Imaging Sciences 5, 203. doi:10.1137/100811593
9. Fruchter, A. S., Hook, R. N. (2002). Drizzle: a method for the linear reconstruction of undersampled images. PASP 114, 144. arXiv:astro-ph/9808087
10. Geringer-Sameth, A., Golovich, N., Iwabuchi, K. (2025). Multi-year stacking searches for solar system bodies. Preprint. arXiv:2509.25428
11. Gupta, V., Bannister, K. W., Flynn, C., James, C. (2024). Efficient summation of arbitrary masks: ESAM. Preprint. arXiv:2412.10678
12. Halko, N., Martinsson, P.-G., Tropp, J. A. (2011). Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review. arXiv:0909.4061
13. Hassanieh, H., Indyk, P., Katabi, D., Price, E. (2012). Nearly optimal sparse Fourier transform. STOC. arXiv:1201.2501
14. Hickson, P. (2018). A fast algorithm for the detection of faint orbital debris tracks in optical images. Advances in Space Research. arXiv:1809.00239
15. Hu, J., Fomel, S., Demanet, L., Ying, L. (2013). A fast butterfly algorithm for generalized Radon transforms. Geophysics 78, U41. doi:10.1190/geo2012-0240.1
16. Kazimirov, D., Nikolaev, D. (2025). A high-accuracy fast Hough transform with linear-log-cubed computational complexity for arbitrary-shaped images. Preprint. arXiv:2509.00231
17. Kielstra, P. M., Shi, T., Luo, H., Qian, J., Liu, Y. (2024). A linear-complexity tensor butterfly algorithm for compressing high-dimensional oscillatory integral operators. Preprint. arXiv:2411.03029
18. Konno, R., Ofek, E. O., Krassilchtchikov, A. et al. (2026). The Large Array Survey Telescope: pipeline. II. Image subtraction and transient detection. Astronomy & Astrophysics. arXiv:2604.27921
19. Lang, D. (2014). unWISE: unblurred coadds of the WISE imaging. Astronomical Journal 147, 108. arXiv:1405.0308
20. Lang, D., Hogg, D. W. (2026; preprint 2020). Principled point-source detection in collections of astronomical images. Astronomical Journal. arXiv:2012.15836
21. Li, Y., Yang, H., Martin, E. R., Ho, K. L., Ying, L. (2015). Butterfly factorization. Multiscale Modeling & Simulation 13, 714. arXiv:1502.01379
22. Mandelbaum, R., Jarvis, M., Lupton, R. H. et al. (2023). PSFs of coadded images. Open Journal of Astrophysics. arXiv:2209.09253
23. Marichal-Hernández, J. G., Gómez-Cárdenes, Ó., Rosa, F., Kim, D. H., Rodríguez-Ramos, J. M. (2020). Three-dimensional multiscale discrete Radon and John transforms. Optical Engineering 59, 093104. arXiv:2501.13664
24. Matthews, S., Nguyen, T. (2026). Performance assessment of track-before-detect algorithms and post-search refinement methods. AMOS Conference. amostech.com
25. Nguyen, T., Woods, D. F., Ruprecht, J., Birge, J. (2024). Efficient search and detection of faint moving objects in image data. Astronomical Journal 167, 113. doi:10.3847/1538-3881/ad20e0
26. Nir, G., Zackay, B., Ofek, E. O. (2018). Optimal and efficient streak detection in astronomical images. Astronomical Journal. arXiv:1806.04204
27. Press, W. H. (2006). Discrete Radon transform has an exact, fast inverse and generalizes to operations other than sums along lines. PNAS 103, 19249. doi:10.1073/pnas.0609228103
28. Rowe, B., Hirata, C., Rhodes, J. (2011). Optimal linear image combination. Astrophysical Journal 741, 46. arXiv:1105.2852
29. Rubin Observatory (2025). Image coaddition. Data Preview 1 documentation. dp1.lsst.io
30. Shih, Y.-h., Wright, G., Andén, J., Blaschke, J., Barnett, A. H. (2021). cuFINUFFT: a load-balanced GPU library for general-purpose nonuniform FFTs. IPDPS workshops. arXiv:2102.08463
31. Smith, K. M. (2016). New algorithms for radio pulsar search. Preprint. arXiv:1610.06831
32. Smotherman, H., Connolly, A. J., Kalmbach, J. B. et al. (2021). Sifting through the static: moving object detection in difference images. Astronomical Journal. arXiv:2109.03296
33. Whidden, P. J. et al. (2019). Fast algorithms for slow moving asteroids: constraints on the distribution of Kuiper Belt Objects. Astronomical Journal 157, 119. arXiv:1901.02492
34. Woodruff, D. P. (2014). Sketching as a tool for numerical linear algebra. Foundations and Trends in Theoretical Computer Science 10, 1. arXiv:1411.4357
35. Zackay, B., Ofek, E. O. (2017a). How to coadd images? I. Optimal source detection and photometry using ensembles of images. Astrophysical Journal 836, 187. arXiv:1512.06872
36. Zackay, B., Ofek, E. O. (2017b). How to coadd images? II. A coaddition image that is optimal for any purpose in the background dominated noise limit. Astrophysical Journal 836, 188. arXiv:1512.06879
37. Zackay, B., Ofek, E. O. (2017c). An accurate and efficient algorithm for detection of radio bursts with an unknown dispersion measure, for single-dish telescopes and interferometers. Astrophysical Journal 835, 11. arXiv:1411.5373
38. Finding Planets project (2026). Faint-mover search frontier: separable tensor-core, butterfly, tree and normalisation design studies. Unpublished design studies; roofline models and CPU tests, not GPU measurements.