Reusing the sums: trees and catalogues

In 2010 Alex Parker and JJ Kavelaars worked out what it would cost to stack LSST images along trial orbits for a 56-hour baseline. Their answer was $6.4\times10^{16}$ pixel additions over 7,743 motion vectors, to reach a forecast depth of $r'=28.6$ for one selected population of distant bodies (Parker & Kavelaars 2010). The number is large, but most of it is redundant. Two trial orbits that end a few pixels apart after 56 hours are indistinguishable over any single night. Each night’s stack along “both” of them is the same stack, computed twice. Parker and Kavelaars saw this and proposed computing each nightly stack once and reusing it. They called the idea a Grid Tree, credited it to a 2002 thesis by R. L. Allen, and stated plainly that they had not implemented it.

This chapter is about the two families of shortcut that follow from that observation, and about what each one costs. The first family keeps the pixels and reuses partial sums: trees, recursive transforms, and hierarchies that prune. The second family throws most of the pixels away. It reduces each image to a catalogue of faint peaks, perhaps a hundred times smaller than the image, and stacks the catalogue instead. Both are real, published, and demonstrated on data. Both deliver large speedups in the papers that introduce them. And in both cases the speedup arrives together with a loss of sensitivity that has to be read off the same page. A tree whose digital paths drift away from the physical track, or a catalogue whose threshold discards the negative fluctuations, sees less signal than a full matched filter. Whether the trade is good depends on numbers that are often reported in different units, on different hardware, against different baselines. Much of the work of this chapter is lining those numbers up.

The background is Chapter 17 (shift-and-stack and KBMOD) and Chapter 20 (orbit-space banks). The template-bank geometry that sets every trial count here is developed in Chapter 38, and the radio-astronomy transforms that inspired the trees in Chapter 36.

Three ideas that share one name

The word “hierarchical” covers three different things, and most confusion about tree stacking comes from mixing them up.

A full-sum tree keeps every partial sum it needs and only changes the order in which additions are done. Like the fast Fourier transform, it can be exactly equivalent to the direct calculation: same templates, same scores, fewer operations.

A grid tree goes further. It uses the coarser motion resolution that suits a shorter time interval. A night needs fewer distinguishable velocities than a month, so the nightly stacks are computed on a coarser grid and each long trajectory is assigned to its nearest nightly representative. That is an approximation, and it needs a mismatch budget: a stated maximum loss of signal-to-noise ratio (SNR) from using a nearby template instead of the true one.

A pruning hierarchy discards whole regions of parameter space that look unpromising at an early level and refines only the survivors. On top of a mismatch budget it needs a false-dismissal budget: a stated probability of throwing away a source that the full search would have found.

The catalogue methods in the second half of the chapter belong to a fourth category: lossy compression of the input. They can be combined with any of the three.

The statistic that has to survive

Every method here tries to compute, or approximate, one quantity: the matched-filter likelihood of a moving point source. It is worth writing down once, because it tells us which shortcuts are free and which are not.

Take registered difference images $d_i$ at times $t_i$, with pixel covariance $C_i$, point-spread function (PSF) template $p_i(q)$ centred at position $q$, and relative photometric response $a_i$. A source of flux $f$ on trajectory $q_i(\theta)$ produces

Define for each exposure the PSF-weighted image and its information,

and sum them along the trajectory: $A(\theta)=\sum_i\psi_i(q_i(\theta))$, $B(\theta)=\sum_i\phi_i(q_i(\theta))$. Then, up to a constant, $\log\mathcal L=fA-\tfrac12f^2B$. The best-fit flux is $\hat f=A/B$, and the detection statistic is $Z=A/\sqrt B$, which has unit variance under the null for a fixed template. This is the sufficient-statistic form that KBMOD implements (Whidden et al. 2019). With diagonal noise, $\psi_i$ is the inverse-variance-weighted PSF correlation of the image and $\phi_i$ is the inverse variance correlated with the squared PSF. Masked pixels simply have zero inverse variance.

The property that makes trees possible is additivity. For any two disjoint sets of exposures $L$ and $R$,

This identity is exact. It says that a partial stack over one night is a legitimate piece of the full stack, provided we store the right thing. We must store the signed numerator $A$ and the information $B$, not a per-night significance, not a clipped positive flux, not a maximum over motions, and not a squared SNR. A night in which the faint source happened to land on a negative noise fluctuation still carries evidence, and that evidence has to reach the final sum with its sign intact. A weighted mean is nestable only if its total weight travels with it. A median is not nestable at all. Parker and Kavelaars made exactly this distinction between weighted-average coadds and median coadds (2010).

For $N$ identical exposures and constant flux the ideal SNR grows as $\sqrt N$. In general it is $f\sqrt{\sum_i\phi_i}$. A tree changes how this sum is computed. It cannot change the ceiling.

Why shorter intervals need fewer velocities

The redundancy that trees exploit comes from how motion resolution depends on time. For a circular Gaussian PSF of dispersion $\sigma$, a template displaced by $\delta$ from the true source has normalised overlap $R=\exp[-|\delta|^2/(4\sigma^2)]$. Over many exposures the overlap is approximately the information-weighted average of this factor, and for small parameter errors $\Delta\theta$ it becomes a quadratic form,

where $J_i$ is the Jacobian of position with respect to parameters at exposure $i$ and the weights $w_i$ sum to one. This is the template-bank metric of Chapter 38, applied to orbits by Geringer-Sameth et al. (2025).

For linear motion over a duration $T$, the velocity spacing at fixed loss scales as $\Delta v\propto\sigma/T$. A fixed two-dimensional region of physical velocities therefore needs $V\propto T^2$ trial velocities. Halve the interval and you need a quarter as many. With acceleration as a free parameter, the acceleration spacing scales as $T^{-2}$ and the count grows much faster still. This is the source of all reuse: short segments have few distinguishable motions, so many long trajectories collapse onto the same short pieces.

A history in three threads

Allen, Parker and Kavelaars

The astronomical lineage begins with deep “pencil-beam” searches for trans-Neptunian objects (TNOs): many short exposures of one small field, stacked along guessed motions to reach below the single-exposure limit (Chapter 17, Chapter 19). Once those searches spanned several nights, the motion grid needed for the full baseline was far finer than any single night could resolve.

Radio precedents

Radio astronomers had been summing along curved tracks recursively for decades. A pulse from a distant source arrives later at low frequencies because of dispersion in the interstellar plasma, so a search must sum a dynamic spectrum along many candidate delay curves. Taylor’s 1974 tree dedispersion combined frequency channels recursively along digital approximations to those curves. Its straight-line pieces mismatch the true $f^{-2}$ dispersion law. In 2014 Barak Zackay and Eran Ofek published the fast dispersion measure transform (FDMT), a recurrence tied to the actual dispersion curve, with an explicit treatment of discretisation (Zackay & Ofek 2014). It is discussed in Chapter 36.

Kendrick Smith’s 2016 coherent pulsar-search tree is the closest mathematical precedent for optical trees that interpolate. It combines statistics from shorter intervals, tabulated on coarser local grids, and reaches constant asymptotic work per final template when the template bank dominates other costs (Smith 2016). His Appendix A shows how transposing the search operator reveals the effective template, and so the accumulated interpolation error, that the tree actually applies. In 2018 Guy Nir, Zackay and Ofek applied recursive line summation to streaks inside a single astronomical image, with $O(L^2\log L)$ cost against $O(L^3)$ for direct line sums on an $L\times L$ image (Nir, Zackay & Ofek 2018); see Chapter 11.

FaXT and its scrutiny

The earliest recursive stacker for faint moving objects that this book’s sources show being built and applied to survey data is the fast discrete X-ray transform, FaXT. Tam Nguyen, Deborah Woods, Jessica Ruprecht and Jonathan Birge published it in the Astronomical Journal in February 2024 (Nguyen et al. 2024). Its abstract describes staged pairwise shift-and-add with reuse of short tracks, and reports more than a thousand asteroids detected in TESS images below the single-frame detection limit, using “moderate computational resources.” That establishes a real-data application. It is a count of detections, not of new discoveries. The full text could not be retrieved for this book, so the TESS sector, limiting magnitude, hardware and run time are not given here.

The clearest look at FaXT’s speed–sensitivity trade comes from Sophia Matthews and Nguyen at the AMOS conference in 2026 (Matthews & Nguyen 2026). They find it about 100 times faster than nearest-pixel shift-and-add at 256 frames, and they also find that it loses more signal. Those two results belong together, and the sections below explain why.

Three other recent works supply ingredients. Vivek Gupta and colleagues’ ESAM (2024) compiles an arbitrary set of radio masks into reusable exact partial sums (Gupta et al. 2024). Geringer-Sameth, Golovich and Iwabuchi (2025) quantify orbit-bank mismatch for multi-year optical stacks and propose segmentation (2025). Pravir Kumar and Zackay’s LOKI (2026) develops calibrated, staged pruning for coherent pulsar searches on GPUs (Kumar & Zackay 2026).

The Grid Tree, and what its savings count

Parker and Kavelaars’s construction has two separate savings, and it is worth keeping them apart.

The first is choosing the grid. On a two-dimensional plane of endpoint displacements, covering every point within a maximum error $\epsilon$ with square cells needs one grid point per area $2\epsilon^2$. A triangular lattice needs one per $(3\sqrt3/2)\epsilon^2$. The ratio of point counts is $4/(3\sqrt3)\simeq0.770$, a 23% saving. The paper reports about this improvement, and stresses a related point about coordinates: grids laid out in fixed steps of speed and angle oversample slow motions and undersample fast ones.

The second saving is the tree. Here is the paper’s cost argument, rewritten. Let $P$ be the number of starting pixels, $n$ the number of nights, $e$ the number of exposures per night, and $F$ the number of motion vectors needed for the full baseline. Direct stacking costs about $PFne$ additions. The paper’s illustrative schedule has nights and days of equal length, so the total baseline is $(2n-1)$ night-lengths. The tree introduces two levels of approximation, and the worst-case position errors at the two levels add. To keep the same total budget, each level gets half the error, which multiplies the density of each two-dimensional grid by four. So the top-level grid needs about $4F$ vectors, and each night, being $(2n-1)$ times shorter, needs $(2n-1)^2$ times fewer:

The total cost is the nightly stacking plus the top-level merge:

When there are many exposures per night, $e\gg(2n-1)^2$, the first term dominates and

For two nights this is $9/4=2.25$; for three, $25/4=6.25$. These are the published figures. Three things deserve notice. They count additions only. The final grid has four times as many trajectories as the flat one, so there are four times as many stack images to inspect and four times as many statistical trials. And the large-$e$ condition is essential: dropping the merge term without it exaggerates the gain. For more than two levels, splitting the error equally between levels quickly becomes expensive, and the right hierarchy has to be optimised jointly over segmentation, per-level resolution and memory, not by applying the two-level gain repeatedly.

The LSST simulation that motivated all this is plotted below. It excludes the tree. It is still useful, because it shows how expensive depth is: going from an 8-hour to a 56-hour baseline at opposition gained 0.6 magnitudes and cost about 120 times more arithmetic.

FaXT: recursive digital straight tracks

FaXT takes reuse to its logical end. Instead of two levels, it uses $\log_2N$ levels, halving the time interval at each step. The cost of doing so is that its tracks are not quite straight.

The recurrence

Consider $N=2^m$ equally spaced frames and, in one coordinate, an integer endpoint displacement $d$ between $0$ and $N-1$ pixels. Nearest-pixel shift-and-add follows the path $n_N(k,d)=\operatorname{round}(kd/(N-1))$ for frames $k=0,\dots,N-1$. FaXT instead defines its path recursively:

In words: the path of length $N$ with displacement $d$ is the path of length $N/2$ with displacement $\lfloor d/2\rfloor$, followed by the same half-length path shifted by $\lceil d/2\rceil$. The path ends at the right place, and the two halves are the same child, so its partial sum can be computed once and used by every parent that needs it. In two dimensions the stack obeys

with $a=\lfloor d/2\rfloor$ and $b=\lceil d/2\rceil$ componentwise; the right child is the later half in time. (This zero-based form follows Matthews and Nguyen’s Eqs. 3–5.) For a likelihood search the same recurrence would be applied separately to the $\psi$ and $\phi$ maps. The published comparison uses aperture sums, not this full noise-weighted model.

The catch is that the recursive path is not the nearest digital approximation of the straight line. Matthews and Nguyen find that the two coincide for every allowed displacement up to eight frames and diverge beyond that. The figure shows a 16-frame example with a 3-pixel displacement. The nearest-pixel path steps at frames 3, 8 and 13. The FaXT path steps at 4, 8 and 12, because each half inherits the same child. At frame 3 the physical position is 0.6 pixels, the nearest path is at 1, and FaXT is at 0.

Counting the work

Let $P$ be the number of starting positions and $V=N^2$ the number of final joint displacements, so the search has $M=PV$ trials. A level with blocks of length $k$ has $N/k$ blocks, each with $k^2$ motions, so it does $Nk$ merges per origin. Summing over the dyadic levels $k=2,4,\dots,N$,

Direct stacking visits $NM$ epoch–trial pairs. FaXT does about $2M$ merges. The arithmetic ratio is $N/c_F$, about 128.5 at $N=256$. If image size also scales as $P\sim N^2$, the familiar statements follow: direct stacking is $O(N^5)$ and FaXT $O(N^4)$. Those exponents assume that image side, frame count and both velocity ranges all grow together. They are not complexity bounds for a fixed image and bank.

Memory traffic gives a less flattering ratio. Carrying two FP32 channels ($A$ and $B$), one merge reads two 8-byte pairs and writes one: 24 bytes. A streamed direct gather reads 8 bytes per epoch–trial pair and keeps its accumulators in registers. Ignoring caches and coordinate arithmetic,

At 256 frames that is 42.8, a third of the arithmetic ratio. On a GPU, where such a low-arithmetic-intensity kernel is limited by memory bandwidth rather than by additions (Chapter 40), the traffic ratio is the more relevant screening number. Neither ratio is a measured speedup.

What the 2026 assessment measured

Matthews and Nguyen simulated 1,000 image stacks each at 64, 128 and 256 frames, with a per-frame SNR of 4 and random velocities. Because the nearest-pixel search was slow, only 240 stacks were run with it at 256 frames. FaXT was about 100 times faster at 256 frames. The processor is not stated. The mean photometric loss of nearest-pixel stacking stayed between 27% and 28%, while FaXT’s grew from 32% at 64 frames to 37% at 256 (Matthews & Nguyen 2026, §6.2).

These loss percentages need a careful reading. The paper’s simulations use a one-pixel detection aperture, and with its modelled Gaussian pixel response (dispersion 0.5 pixel) that aperture can hold at most about 47% of the source’s light, so a large part of both methods’ loss is the aperture, not the tracking. The difference between the methods, five to ten percentage points, is the part due to FaXT’s path. In a separate faint-source test at 64 frames, with a per-frame SNR of 0.875 chosen so that the ideal stacked SNR is 7, the detection rate at threshold 7 was about 6% for nearest-pixel and 3% for FaXT. These are aperture-based losses on undersampled images. They cannot be converted directly into matched-filter overlap.

The paper also tests refinement: fitting a physical moving PSF to the retained detections. Refinement recovers almost everything for the detections that survive. In a geosynchronous-orbit case study, retained SNR before refinement was 0.50 for nearest-pixel and 0.48 for FaXT, and 0.94 for both afterwards. But refinement has a price, and it cannot help detections that never crossed the threshold.

In that example, search plus refinement for FaXT took 26 units, more than the unrefined nearest-pixel search. In the paper’s low-Earth-orbit example, two FaXT scores were 6.74 and 6.98 where nearest-pixel gave 7.37 and 7.22. At a threshold of 7, FaXT lost both before refinement could see them. In the 256-frame runtime test, on the other hand, FaXT plus refinement ran well within the nearest-pixel time. The honest summary is that FaXT’s search is much faster, its threshold sees less signal, and whether the whole pipeline is faster depends on how many candidates go to refinement.

The book project’s own design studies looked at the same question with a matched-filter measure instead of aperture flux. They are unpublished CPU experiments and models, and are reported here as such. With pixel-integrated Gaussian PSFs of dispersion 1.2 pixels at 256 epochs, the worst selected FaXT template responded at 0.956 of the nearest-pixel template, an extra 4.35% SNR loss, or 0.048 mag. At dispersion 0.8 pixel the ratio was 0.917. Searching neighbouring origins over the same displacement bank did not remove the worst case. Against a sensitivity budget of 2% SNR, which is $-2.5\log_{10}0.98=0.0219$ mag, unmodified FaXT does not qualify at these frame counts.

Assumptions hidden in the frame index

FaXT’s recurrence runs over frame indices. That assumes equally spaced exposures and straight-line motion with displacements in a fixed non-negative range. Negative motions need reflected copies of the transform. Motions outside the range need other banks. Padding the frame count to a power of two does not make irregular timestamps regular. A missing image can contribute zero information, but the path is still defined by frame index, not by the true time. Real survey cadences, with gaps, unequal exposures, parallax and curvature, therefore need either a different recurrence or a different tree.

Exact trees: compile the paths you actually need

One alternative keeps a discrete trajectory bank exactly and reuses only subtracks that really are identical. ESAM showed how to do this for arbitrary radio masks: each internal product is identified by two child identifiers and an offset, and duplicate products are stored once (Gupta et al. 2024). Its Python reference implementation needed more than ten times fewer operations than direct dedispersion on its example bank, while reproducing the supplied masks exactly.

For images, the analogue is to remove the common translation from each temporal subtrack, hash the remaining relative path, and intern child pairs recursively:

compile(track, interval):
    remove the track's common integer translation over this interval
    if interval is one exposure: return leaf id
    left  = compile(track restricted to first half)
    right = compile(track restricted to second half)
    key   = (left, right, relative translation, weight metadata)
    return existing product with this key, or create it

The result is a directed acyclic graph rather than a tree, because many parents share a child. Evaluating $H$ unique products over $P$ origins costs $O(PH)$. There is no guarantee that $H$ stays close to the number of final motions. Arbitrary masks, weights and irregular cadence can sharply reduce the number of identical subtracks.

The project’s exact-path prototype measured this. With uniform cadence at 256 frames, the bank has 65,536 final joint motions but needs 1,476,320 internal states, 22.5 per final motion, against FaXT’s 1.99. The logical traffic ratio falls from FaXT’s 42.8 to 3.79. Gaps do not always hurt. A schedule with two observing seasons needed only 7.66 states per motion, a traffic ratio of 11.1, because many within-season paths become identical. These are prototype counts and derived ratios, not a lower bound for all exact algorithms.

The two figures make the central trade visible. FaXT’s economy comes from forcing every path into a shape that shares children. Keeping the physical path costs most of that economy back.

Interpolating trees in local coordinates

A second alternative keeps a physical family of curves and interpolates the shorter-interval statistics. For a segment $s$, map the parent’s parameters into each child’s local coordinates with transformations $F_{sL}$ and $F_{sR}$:

This identity is exact for continuous child surfaces. Storing the children on a finite grid and interpolating introduces the error. Smith’s pulsar tree uses exactly this structure, with local Legendre-polynomial coordinates (Smith 2016).

The project’s optical study writes each spatial coordinate in a segment centred at $c_s$ with half-span $h_s$ as $x(t)=u_0+u_1z+u_2(z^2-1/3)$, where $z=(t-c_s)/h_s$. For a child with centre offset $d=(c_c-c_s)/h_s$ and half-span ratio $r=h_c/h_s$, the child’s coefficients are

These identities are exact for a quadratic and work for irregular exposure times. What is approximate is the interpolation, and whether a quadratic is an adequate model of the motion. A single quadratic over years is not a Keplerian orbit. A useful optical tree needs local orbit Jacobians or an explicit parallax basis.

Why the geometric series is not the whole story

In $d$ spatial dimensions with independent motion derivatives up to order $q$, the order-$k$ coefficient count at fixed mismatch scales as $T^{dk}$, so the bank scales as $T^p$ with

If the full-span bank has $M$ states including position, each halving of the interval multiplies the number of nodes by two and shrinks each node’s bank by $2^{-p}$. Summed over depths, the ideal total for $p>1$ is

Linear motion in two dimensions has $p=2$ and $I\simeq2M$, matching FaXT’s $c_F\to2$. Quadratic motion has $p=6$ and $I\simeq32M/31$, barely more than the final bank. That looks almost free. It is not, because there are now six coordinates including position, and interpolating in six dimensions is expensive. Multilinear interpolation reads $2^6=64$ table entries per child. A four-point cubic reads $4^6=4{,}096$. Oversampling each intermediate axis by a factor $o$ to control interpolation error multiplies the intermediate density by $o^6$:

which is 3.07 at $o=2$ and 133 at $o=4$. The geometric series is real, but so is the curse of dimensionality in the interpolation stencil.

The project’s CPU experiments put numbers on the sensitivity side, for a two-dimensional Gaussian PSF. Native multilinear interpolation lost 2.5% to 4.2% of expected SNR across the tested cadences and frame counts. Twice-oversampled linear interpolation lost at most 0.63%. Native cubic interpolation reached about 0.42% worst-case loss but at many more taps. With a narrower PSF (dispersion 0.8 pixel), native cubic loss rose to 2.8%, a reminder that results on well-sampled PSFs do not transfer to undersampled instruments. On speed, a straightforward table implementation was modelled at only 0.55–0.94 times the speed of an ideal direct gather at 1,024 exposures, and a proposed tiled schedule at 5.6–18.9 times, before calibration and support costs. No GPU was used for any of these numbers.

The error budget

A hierarchical approximation has at least four separate sources of loss: finite template spacing, digital path quantisation, interpolation of partial tables, and numerical precision. Pruning adds a fifth, the chance of discarding a detectable source. Each should be allocated and measured separately.

Worst-case geometric errors add across levels, $\epsilon_{\rm total}\le\sum_\ell\epsilon_\ell$. They cannot be added in quadrature as if they were independent random errors. That is why Parker and Kavelaars split the budget between levels, and why their final grid got four times denser.

Interpolation raises a subtler problem: normalisation. If the tree estimates $\widetilde A=\sum_jc_jA(q_j)$ from neighbouring table entries, its null variance is

which is generally neither $\sum_jc_jB(q_j)$ nor $\sum_jc_j^2B(q_j)$, because neighbouring matched-filter values share pixels. The correct normalisation uses the norm of the effective template $h$ that the whole tree applies. For whitened true template $e$, the response is

Applying the transpose of the full linear search operator to one output cell reconstructs $h$ without building every output. This is Smith’s method, which the project’s full-sum study applied to optical trees. The study found that interpolating the information $B$ as if it were an ordinary image gave null standard deviations of 0.75–0.80 for native linear interpolation and 0.90–0.93 for twice-sampled linear interpolation, instead of 1. A score whose null is that far from unit variance has badly miscalibrated tails. Exact refinement at the end fixes the scores of retained candidates. It cannot fix a misnormalised early rejection.

Trials do not shrink

For $M$ tests with unit-normal one-sided scores, the union bound $P_0(\max Z>u)\le M\,\bar\Phi(u)$ holds even when the tests are correlated. For $M=10^{12}$ and a global false-alarm probability of 0.01, it requires $u\simeq7.65$. A tree does not reduce $M$: it evaluates the same templates faster. But a tree with smoother digital paths can change how correlated neighbouring scores are, so two methods with the same nominal threshold need not have the same actual false-alarm rate. Matthews and Nguyen quote a Gaussian false-positive probability of $3.58\times10^{-5}$ per $84\times84\times64$ stack at threshold 7. That is a nominal calculation, not a measured artefact rate. A real comparison needs injected sources and null runs through the whole pipeline. Chapter 1 covers the statistics.

Pruning hierarchies

A pruning hierarchy combines reusable sums with selection. Start from short segments or coarse parameter cells, accumulate their evidence, refine the grid where it looks promising, and drop the rest. The design rule that matters most is that early stages must not require an individual detection. A faint object’s first night may well be consistent with noise. Requiring each night to cross a threshold, or taking the maximum over motions within each night before combining, changes the statistic and gives up the $\sqrt N$ gain. The gates have to act on signed, accumulated evidence, and be permissive enough that their false-dismissal rate is small and measured.

Smith distinguishes the fully coherent tree from semicoherent and hierarchical searches. LOKI extends coherent accumulation with optimised thresholds and GPU batching; its implementation is C++20 with OpenMP and CUDA, and its published benchmarks ran on an NVIDIA L40S (Kumar & Zackay 2026). These are pulsar results. They establish an architecture, not an optical speed factor.

The project’s optical transfer of the idea keeps signed nightly statistics, calibrates permissive gates, and rescores survivors on the original maps. In a CPU image experiment with 480 epochs over 20 irregular nights, 37,632 fine trajectories and 192 injected cubes, it made 19.2 times fewer statistic accesses than exhaustive stacking and 8.6 times fewer than a sequential-rejection control, with identical blind recoveries in the sample. The measured ratio of local nightly table work to global gather work was 0.027. A more favourable production assumption of 0.005 led to a modelled H100 gain of 10.3 times over one comparator, falling to 6.8 after charging an extra normalisation pass and 6.0 against a stronger sequential comparator. The measured and favourable geometries are different numbers and should not be conflated.

A later project study counted that ratio on real cadences: actual Subaru Hyper Suprime-Cam timestamps from New Horizons search campaigns in 2020 and 2024, with topocentric Earth geometry and propagated Keplerian orbits. For month-scale banks of TNOs at 30–100 au the counted ratios were tiny, between $1.3\times10^{-9}$ and $3.7\times10^{-5}$, because a month resolves angular rates and curvature that a four-hour night cannot. That removes one objection. It does not settle the speed question: the modelled factor against a strong sequential baseline was 15–23 times in central cases, 3–4 under adverse hardware assumptions, and 7.3 for the hardest central case once a masked-normalisation pass was charged. All of these are models. The tables that are cheap to build have to be queried, compacted and normalised, and those costs now dominate.

Catalogue-domain stacking

The second family takes a different view of the redundancy. Most pixels in a difference image are noise, and most trial trajectories pass through nothing. Instead of summing pixels, find the local peaks in each image at a very low threshold, keep their positions and times, and look for sets of peaks that line up along a common motion. A faint body produces a statistically coherent sequence of peaks even when no single exposure could establish its existence. Random noise peaks rarely line up.

From bubble chambers to asteroid catalogues

The tool for “find the parameters that many points agree on” is the Hough transform, invented to find particle tracks in bubble-chamber photographs. Each observation votes for every parameter combination that could explain it, and true structures appear as peaks in the vote accumulator. Christoph Dalitz, Tilman Schramke and Manuel Jeltsch’s 2017 iterative Hough transform for lines in three-dimensional point clouds, developed on particle-track data among other examples, added two important steps: refine each coarse peak with orthogonal least squares, then subtract the votes of the points assigned to it before looking for the next (Dalitz, Schramke & Jeltsch 2017).

The statistical foundation for using faint catalogues came from Tamás Budavári, Alexander Szalay and Thomas Loredo in 2016–2017 (Budavári, Szalay & Loredo 2017). They asked whether a survey could publish low-threshold catalogues and recover faint stationary sources by combining them. In an idealised simulation with nine epochs and a one-sigma threshold, the answer was yes, provided both flux and position agreement were used. Two details of their model matter here. The noise hypothesis describes detected peaks, whose heights are biased upward compared with Gaussian pixels. And epochs with no catalogue entry also carry evidence.

Pan-STARRS supplied a practical demonstration. Kai-Jie Lo, Chan-Kao Chang, Hsing-Wen Lin and colleagues reused six nights of high-cadence Pan-STARRS1 observations from October 2016, taken for asteroid rotation periods, with about thirty exposures per field per night (Lo et al. 2020). Overlaying each night’s catalogue positions in a “master frame” made asteroid tracks visible as straight lines, which a two-dimensional Hough transform extracted. They recovered 3,819 of 3,870 known asteroids accessible to the comparison, and reported 3,574 previously unrecognised asteroid tracks, mostly at $21.5\lesssim w_{P1}\lesssim22.5$. The authors caution that some short tracks may be pieces of the same body and that short arcs hinder recovery, so this is a count of reported tracks, not of secure orbits.

Pedro Bernardinelli, Gary Bernstein, Masao Sako and the Dark Energy Survey (DES) team showed the other half of the idea: confirm catalogue candidates with pixels (Bernardinelli et al. 2020). Candidate TNO orbits came from catalogue detections. Each was then tested by measuring flux at the predicted positions in images where the object had not been detected. This sub-threshold significance used information not involved in proposing the orbit. Of 424 candidate linkages, 316 passed, including 245 DES discoveries and 139 first reported in that paper.

Steven Stetzler, Mario Jurić, Pedro Bernardinelli, Dino Bektešević, Colin Orion Chandler, Andrew Connolly and the DECam Ecliptic Exploration Project (DEEP) collaboration put the threads together in An Efficient Shift-and-Stack Algorithm Applied to Detection Catalogs, posted on 30 September 2025 (Stetzler et al. 2025). They took one long stare from DEEP, field A0c on 3 April 2019: 104 exposures of 120 seconds each in the VR band. They injected synthetic asteroids and TNOs, built source catalogues at thresholds from one to five sigma, and searched them. The same data had been searched with GPU KBMOD. The result is the most complete published measurement of what a catalogue search gains and loses against a pixel search on the same images.

How much signal survives a threshold

A source with true single-exposure SNR $\nu_i$ appears in a catalogue with threshold $\tau$ whenever noise lifts its measured SNR above $\tau$. In the Gaussian approximation,

A source need not be above threshold on average to appear sometimes. Combining those intermittent appearances is what lets a catalogue search beat the single-exposure limit. It also biases the measured fluxes upward.

The peak also has to land where the search expects it. For a circular Gaussian PSF with Gaussian scale $a_i$, Stetzler uses the positional variance $\sigma_{x,i}^2=\sigma_{y,i}^2\simeq2a_i^2/\nu_i^2$. If the true projected position sits at the centre of a square bin of width $d$, the chance of being caught is

The number of signal votes $S=\sum_iX_i$, with $X_i$ a Bernoulli variable of probability $p_i$, is Poisson-binomial with mean $\sum_ip_i$. The figure evaluates this for 104 exposures. At a one-sigma threshold, a source with true SNR 0.5 is expected to leave 32 peaks, but only about 5 of them land in a centred bin of width $3a$. At three sigma it leaves under one.

This is an instructive approximation, not an exact likelihood for faint peaks. It assumes the right velocity, centred bins, Gaussian centroid errors and independent selection and position. The divergence of $\sigma_x$ as $\nu\to0$ is a diagnosis: at very low SNR the measured peak need not be a perturbed image of the source at all. Widening the bin catches more true positions but also collects more noise peaks, over an area $d^2$.

The noise side is a count of local maxima of a smoothed random field. With $\lambda_i(\tau)$ the surface density of noise peaks above threshold and $A_i$ the valid area of the projected bin ($d^2$ with full coverage), the number of noise votes in a cell is approximately $\mathrm{Pois}(\Lambda)$ with $\Lambda=\sum_i\lambda_i(\tau)A_i$. The peak density is not the Gaussian tail probability per raw pixel. Smoothing by the PSF reduces the number of independent opportunities, and peak selection shifts the height distribution upward. Budavári and colleagues give the analytic peak density in their Appendix A, which is the place to implement it from. (The specialised form in Stetzler’s v1 preprint differs algebraically from Budavári’s equations, including a missing square root. This is an observation from comparing the two texts, not an established erratum.)

A count threshold $K$ is then chosen so that $P[\mathrm{Pois}(\Lambda)\ge K]\le N_{\rm false}/H$ for $H$ accumulator cells. That gives the expected number of false cells. The number of distinct false objects also depends on clustering and on how candidates are removed. Smoothed-field peaks are correlated, survey artefacts are worse than Poisson, and empirical null runs remain necessary.

The algorithm

Stetzler’s method projects every detection back to a reference time for each trial velocity, bins the projected positions, and treats the stack of velocity-indexed histograms as one accumulator:

C = all catalogue peaks above threshold tau, from every exposure
V = velocity lattice over the search region, spacing about d / T
H = zero counts over (velocity, x bin, y bin)
for v in V:
    for detection in C:
        u = detection.position - v * (detection.time - t0)
        H[v, bin(u, d)] += 1
while max(H) >= K and fewer than the output cap:
    cell   = argmax(H)
    voters = detections that voted for cell
    remove every vote of these voters from H, at all velocities
    track  = robust MCD regression of voter positions on time
    gather all detections within d of the track; optionally refit
    output track

The fit treats time, which is known accurately, as the independent variable and sky position as the response. It uses minimum-covariance-determinant (MCD) estimates of mean and covariance instead of ordinary least squares. MCD looks for the subset of points whose covariance has the smallest determinant, and with a support fraction near one half it can tolerate nearly half the points being outliers. That robustness is what makes coarse, contaminated clusters usable. It is not a detection guarantee: structured outliers or two real objects in one cluster can defeat it.

The search parameters were as follows (Stetzler et al. 2025, Table 1). Bin widths are in units of the PSF scale $a$.

The TNO velocity search covered 0.003–0.04 deg/day (about 0.45–6 arcsec/hour) in all directions; the main-belt search covered 0.1–0.5 deg/day. Sixty detectors were searched. The depth comparisons used a count cut giving 200 results per detector. Main-belt thresholds below three were not attempted because of their cost.

Where the speed comes from

Let $M$ be the number of catalogue detections, $N_v$ the number of trial velocities, $B\simeq\Omega/d^2$ the number of spatial bins per velocity over a field of solid angle $\Omega$, and $R$ the number of extracted candidates. A straightforward implementation costs

plus vote removal and fitting. The terms are source extraction over $N_t$ images of $P$ pixels, voting, and repeated scans of the accumulator for the maximum. A brute-force pixel search costs $O(N_tPN_{v,\rm img})$.

The gain has two origins. The first is obvious: there are far fewer peaks than pixels. Stetzler’s example $2048\times4096$ detector shows peak counts smaller than pixel counts by factors of roughly $10^2$ to $10^6$ across the plotted thresholds. The second is less obvious and just as large. A velocity grid of spacing $d/T$ over a rectangle of velocity spans $W_x\times W_y$ has $N_v\simeq T^2W_xW_y/d^2$ velocities, and the accumulator has

cells. Because the catalogue cluster is refined afterwards by regression, the search can use a much larger $d$ than a pixel likelihood could tolerate, and the cost falls as $d^{-4}$. As an illustration, not a reported run: for a 162 arcmin² field, a four-hour baseline, a 1-arcsec bin and the all-direction TNO annulus above, the velocity lattice has about 1,800 points and the accumulator about $1.05\times10^9$ cells. A five-times-larger bin cuts this by $5^4=625$, to about $1.7\times10^6$.

The sparse input does not guarantee small memory. Storing every detection’s projected bin for every velocity takes $O(MN_v)$, and the accumulator $O(BN_v)$. At the lowest thresholds, this is what limits the method.

The measured trade

Here is Stetzler’s central table. Search time is the paper’s total catalogue CPU-time column. The comparison baseline is a specially written CPU implementation of the image search’s core loop, which took 2.54 core-hours and 17.5 GB, while the comparison depth $m_{50}=25.47\pm0.01$ comes from the GPU KBMOD run. Speedups and memory ratios keep the paper’s inequality signs, because the CPU baseline excludes reporting and filtering results.

The paper summarises the threshold-two search as a speedup of 30.8 times with 88% of the memory, at a cost of 0.25 mag. (Subtracting the printed numbers gives 0.26; the paper’s printed difference is 0.25.) For reference, the measured single-exposure $m_{50}$ was 23.7, and the ideal gain from 104 equal exposures is $1.25\log_{10}104=2.52$ mag, giving an idealised coadd limit of about 26.2. KBMOD itself reached 25.47, well short of that ideal: a blind search over many trajectories needs a higher threshold than a test of one known trajectory, and its candidate filtering has losses of its own.

What does 0.25 mag cost in observing time? A source must be $10^{0.1}=1.26$ times brighter to be detected, and in the background-limited regime recovering that by integrating longer takes $1.26^2=1.58$ times the exposure. So a 30-fold saving in search compute was bought with the equivalent of a 58% increase in telescope time, or with a loss of the faintest objects. Whether that is a good trade depends on the science: it is cheap for a survey whose bottleneck is compute and expensive for one whose bottleneck is photons.

The paper also measured astrometry. Both methods give sub-pixel trajectories over most of the tested range, but catalogue-fit accuracy worsens at about magnitudes 25–26, while KBMOD stays near one DECam pixel (0.263 arcsec). Faint centroids are noisy and selection-biased, so a regression advantage on bright objects does not carry to the discovery limit.

Forced photometry: go back to the pixels

Once a candidate exists, the pixels are cheap to revisit. Measure the flux $f_i$, with variance $s_i^2$ and expected scaling $g_i$, at the predicted position in every usable image, including negative values and epochs with no detection:

At fixed PSF and position this is the same linear statistic as $A/\sqrt B$. Its attraction is that it is evaluated only for $C$ candidates, at cost $O(CN_t)$ reads from precomputed matched-filter maps. If $C$ grew to the whole trajectory bank, it would be the pixel search again.

DES is the worked example. It measured windowed, sky-subtracted fluxes at predicted positions with a one-arcsec Gaussian window and combined them by inverse variance, using only images without a linked detection and skipping exposures within an hour of one. Candidates were accepted with an empirical boundary $\mathrm{STS}>2\sqrt{N_{\rm images}}$, plus blinded visual inspection. That boundary is specific to the survey, not a universal Gaussian prescription. Two cautions apply to any new catalogue search. Measuring high flux after fitting the trajectory to the same noisy data is not independent evidence: hold out exposures, or calibrate the joint proposal-plus-photometry procedure on null data. And forced photometry can rescue weak epochs of a proposed object but cannot recover one the proposal stage never nominated.

Stacked detections are not stacked images

The phrase “stacked detections” has other meanings in Pan-STARRS work, and they should not be confused. Lo’s master frame is an overlay of catalogue positions. After removing noisy detector regions, points near stationary sources and isolated points, a typical nightly field still held about 260,000 detections. The Hough transform sampled angles from 1 to 179 degrees in 0.1-degree steps and needed at least five detections per segment. It took about twenty-five minutes per master frame on a five-node cluster of four-core Intel i7-2600 machines (Lo et al. 2020). Separately, Pan-STARRS image-stack catalogues come from ordinary coadds; positions found in stacks seed forced photometry on the individual warped exposures (Magnier et al. 2020). That preserves weak evidence at fixed positions, not along unknown moving paths, and a stellar-aligned stack smears a mover. The 2025 Pan-STARRS distant-planet search by Matthew Holman, Kevin Napier, Matthew Payne and Jacob Kurlander used ordinary single-exposure catalogues with a stationary-source veto that removed 80–95% of detections, and found 692 Solar System objects including 642 TNOs (Holman et al. 2025). It is a large catalogue search, but not a sub-threshold one; it is treated in Chapter 9 and Chapter 12.

Landmark results

None of the tree-specific sources supplies a survey-wide count of new bodies attributable to the recursion itself.

Limitations and open problems

The digital path can spend the sensitivity budget before any GPU work begins. FaXT’s economy comes from paths that drift from the physical line. At the frame counts where it is fastest, its extra loss exceeds a 2% SNR budget in the project’s matched-filter tests and is visible in the published aperture losses.

Normalisation is the common failure. Interpolating information maps as if they were images, or ignoring the correlation between neighbouring table entries, miscalibrates the null. Near masks, chip gaps and changing PSFs the correction varies from place to place. A miscalibrated early gate loses sources that exact final scoring cannot bring back.

Real cadences and real orbits are not details. Gaps, unequal exposures, parallax, curvature and distance uncertainty add resolved directions that a frame-index recurrence does not represent. Conversely, physical constraints can reduce a generic polynomial bank, so the pessimistic six-coordinate count above is not an impossibility proof.

Intermediate tables can dominate memory. The project’s audit notes that a $256^2$ origin tile times 65,536 motions, storing one FP32 numerator–information pair per root, already takes 34.4 GB, and the two half-length children add 17.2 GB. Root streaming and tiling are required even when the images fit easily.

Pruning changes completeness, and refinement can eat the savings. A low partial score does not imply a low final score. Lowering a seed threshold to keep faint sources admits rapidly more noise candidates, because the Gaussian tail grows steeply as the threshold falls. Matthews and Nguyen’s GEO case shows refinement costing more than the search it follows.

For catalogues, the threshold is irreversible. No association method can recover a signed flux from a record that says only “no peak above threshold.” Censoring models recover part of the lost evidence, not all of it. The data-processing inequality bounds the best catalogue statistic by the best image statistic under the same model. The size of the loss depends on brightness, PSF sampling, threshold, stored metadata and follow-up access; it is not a universal quarter magnitude.

Artefacts break the clean null. Difference residuals recur at stars, shared templates correlate exposures, and diffraction spikes and detector edges make structured patterns of peaks. A constant-density Poisson model understates such tails. Null tests must preserve the instrument’s structure; scrambling times alone may destroy the very artefacts that create real false candidates.

Greedy removal and bin phase can hide objects. A real cluster can straddle a bin edge. A noisier maximum can remove some of its detections first. Overlapping bins or multiple grid offsets help, at a cost in work and in altered null statistics.

Reporting is incomplete. The original FaXT TESS configuration is not accessible; the detailed 2026 FaXT study is simulated, undersampled and run on unstated hardware; Stetzler does not give a processor model or GPU-hours for the catalogue benchmark; and the project’s GPU numbers are models. These are specific missing measurements, not reasons to dismiss demonstrated reuse.

Where the room is

The physical ceiling is the correctly weighted matched filter on the true trajectory, with the right light-curve model. Trees cannot create photons or restore information removed by masking, trailing or over-aggressive background subtraction. A larger parameter space raises both the bank size and the detection threshold. Everything else here is algorithm and engineering: the exposure loop, repeated coordinate evaluation, duplicate short tracks, intermediate storage, memory traffic, repeated accumulator scans and data movement. Those can change by large factors without new photons. The question is always against what.

A decision rule for trees

For $M$ final tests with an average of $\bar E$ useful exposures each, an ideal FP32 direct search reads about $8\bar EM$ bytes. An exact two-channel tree with $cM$ merge states moves about $24cM$ logical bytes. Before caches and other work, the tree beats the direct search on traffic only when

and wins by a factor of ten only when $c<\bar E/30$. If the direct search reads a two-byte encoded pair while the tree keeps FP32 partials, the tenfold condition tightens to $c<\bar E/120$. These are screening inequalities, but they are cheap to evaluate before writing a kernel. At 256 frames the project’s roofline model gave FaXT 31.2 times over an FP32 direct gather, but only 7.8 times over an encoded one; exact nearest paths gave 2.8 and 0.69 times. The comparator matters as much as the algorithm.

The input-reading bound depends on useful coverage, not nominal frame count. A field with thousands of archived exposures but few overlapping any one output trajectory does not suit a dense tree. Charge the direct search only for the useful intersections $K=\sum_\theta N_{\rm valid}(\theta)$, and trim the tree’s empty states in the same way, or the apparent gain will consist largely of work the baseline never needed to do. This is why dense, long stacks of one field (Chapter 19) are the natural home for trees, and sparse multi-year archives (Chapter 20) are not.

Catalogues: engineering limits first

Most of the catalogue method’s weaknesses are representation and engineering, and several have clear fixes.

The threshold-one memory problem comes from materialising projections for every velocity and a dense accumulator. Streaming velocity tiles, computing projections in registers, voting into a tile and extracting its candidates before discarding it, would move memory from $O(MN_v+BN_v)$ toward $O(M+BN_{v,\rm tile})$. The subtlety is that the paper removes a cluster’s votes globally before choosing the next; independent tiles cannot reproduce that order exactly, so an overcomplete set of per-tile maxima has to be merged and conflicts resolved, and injection tests must confirm that the changed order does not cost completeness.

Repeated full scans for the maximum create the $RBN_v$ term. Per-tile maxima, a heap with lazy validation, or GPU compaction followed by a segmented top-$k$ would remove it. Whether dense histograms or sparse sorted keys win depends on occupancy, which should be measured. On a GPU the bottleneck is scattered histogram updates and membership bookkeeping, not floating-point arithmetic; tensor cores do not accelerate integer binning.

The sensitivity loss is the harder problem, and the most interesting one. Stetzler’s score is an integer count: a bright peak and a barely-above-threshold noise peak each get one vote. Storing each peak’s flux, centroid covariance, local noise and coverage, and accumulating a likelihood-weighted score with the censored evidence of Budavári’s model, would recover some of what equal counting discards. The authors list likelihood accumulation as future work. The weights must be robust, or bright subtraction residuals will dominate. And the right endpoint is a hybrid: a cheap sparse proposal stage tuned for recall at a fixed candidate budget, followed by signed forced photometry or a local pixel likelihood for the survivors. Cutting the proposal set tenfold could plausibly cut the follow-up work tenfold, but that is an engineering target, not a demonstrated speedup, and it holds only if measured recall and final false-alarm rate stay fixed.

The experiment that would settle it

For both families the missing measurement is the same. Take identical calibrated exposures and injections with identical spatial and velocity coverage. Run an optimised dense GPU likelihood search, the published tree or catalogue method, an optimised GPU version of it, and a hybrid. Report completeness against flux and rate at a common false-object density after the same final validation. Include extraction, I/O, transfers, refinement and pixel follow-up in the end-to-end time, and report kernel time, memory, device model and code revision separately. Until that table exists, a “30×” or “100×” is a statement about one program against another, not about what the method can do.

Chapter 42 puts these trees and catalogues alongside the other routes to a large equal-sensitivity speedup. The short version from this chapter: reuse is real and old; its price is paid in path fidelity, normalisation and memory traffic; and the most promising designs keep signed sums all the way up, spend exact pixel work only where cheap evidence points, and measure every stage’s survival probability.

References

1. Bernardinelli, P. H., Bernstein, G. M., Sako, M. et al. (2020). Trans-Neptunian objects found in the first four years of the Dark Energy Survey. ApJS 247, 32. arXiv:1909.01478
2. Budavári, T., Szalay, A. S., Loredo, T. J. (2017). Faint Object Detection in Multi-Epoch Observations via Catalog Data Fusion. ApJ 838, 52. arXiv:1611.03171
3. Dalitz, C., Schramke, T., Jeltsch, M. (2017). Iterative Hough Transform for Line Detection in 3D Point Clouds. Image Processing On Line 7, 184–196. doi:10.5201/ipol.2017.208
4. Geringer-Sameth, A., Golovich, N., Iwabuchi, K. (2025). Multi-year Stacking Searches for Solar System Bodies. arXiv preprint. arXiv:2509.25428
5. Gupta, V., Bannister, K. W., Flynn, C., James, C. (2024). Efficient Summation of Arbitrary Masks: ESAM. arXiv preprint. arXiv:2412.10678
6. Holman, M. J., Napier, K. J., Payne, M. J., Kurlander, J. A. (2025). A Pan-STARRS Search for Distant Planets: Part 1. arXiv preprint. arXiv:2506.02144
7. Kumar, P., Zackay, B. (2026). Coherent Signal Detection with Pruning I: Finding Short-Period Binary Pulsars in Circular Orbits. arXiv preprint. arXiv:2607.07700
8. Lo, K.-J., Chang, C.-K., Lin, H.-W. et al. (2020). Asteroid Discovery and Light Curve Extraction Using the Hough Transform: A Rotation Period Study for Sub-Kilometer Main-Belt Asteroids. AJ 159, 25. arXiv:1910.07146
9. Magnier, E. A., Chambers, K. C., Flewelling, H. A. et al. (2020). The Pan-STARRS Data Processing System. ApJS. arXiv:1612.05240
10. Matthews, S., Nguyen, T. (2026). Performance Assessment of Track-Before-Detect Algorithms and Post-Search Refinement Methods. AMOS Technical Conference. amostech.com
11. Nguyen, T., Woods, D. F., Ruprecht, J., Birge, J. (2024). Efficient Search and Detection of Faint Moving Objects in Image Data. AJ 167, 113. doi:10.3847/1538-3881/ad20e0
12. Nir, G., Zackay, B., Ofek, E. O. (2018). Optimal and Efficient Streak Detection in Astronomical Images. AJ 156, 229. arXiv:1806.04204
13. Parker, A. H., Kavelaars, J. J. (2010). Pencil-Beam Surveys for Trans-Neptunian Objects: Novel Methods for Optimization and Characterization. PASP. arXiv:1004.3287
14. Smith, K. M. (2016). New Algorithms for Radio Pulsar Search. arXiv preprint. arXiv:1610.06831
15. Stetzler, S., Jurić, M., Bernardinelli, P. H. et al. (2025). An Efficient Shift-and-Stack Algorithm Applied to Detection Catalogs. AJ 170, 352. arXiv:2509.26279
16. 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. AJ 157, 119. arXiv:1901.02492
17. Zackay, B., Ofek, E. O. (2017). An Accurate and Efficient Algorithm for Detection of Radio Bursts with an Unknown Dispersion Measure, for Single-Dish Telescopes and Interferometers. ApJ 836, 187. arXiv:1411.5373
18. Finding Planets project (2026). Design studies: full-sum trajectory trees; recursive stacking prior-art audit; LOKI transfer study; regime-D table-ratio test. Unpublished internal reports (CPU experiments and GPU performance models; no GPU measurements).