---
title: 'PyStoch: Efficient Pipeline for SGWB Mapmaking'
url: https://www.emergentmind.com/topics/pystoch
type: topic
---

# PyStoch: Efficient Pipeline for SGWB Mapmaking

Searching arXiv for PyStoch-related papers to ground the article in the cited literature.
PyStoch is an open-source, Python-based pipeline for gravitational-wave radiometry and anisotropic stochastic gravitational-wave background (SGWB) mapmaking that was introduced to exploit sidereal-day data folding, matrix-based computation, and HEALPix-based sky representations in order to make directional cross-correlation analyses computationally practical at scale [1803.08285]. In subsequent work, its role broadened from fast pixel-basis radiometer mapping to a unified framework that also supports spherical-harmonic (SpH) analyses, all-sky all-frequency directional searches, joint multicomponent anisotropy inference, and directed continuous-wave (CW) searches from known sky locations [2011.05969].

## 1. Historical development and problem setting

PyStoch was developed in the context of SGWB searches, where the observable is a weak excess correlation between spatially separated detectors with approximately uncorrelated instrumental noise. Conventional anisotropic SGWB analyses had already used cross-correlation techniques, but they were computationally expensive and, in practice, implemented with distinct numerical pipelines for different sky bases. The 2018 presentation of PyStoch emphasized three design decisions: folding data into one sidereal day, replacing explicit loops with matrix multiplications, and incorporating HEALPix pixelization tools; together, these made the radiometer analysis “a few thousand times” faster and allowed all-sky maps to be produced “in just a few minutes on an ordinary laptop” [1803.08285].

The pipeline’s scope then expanded along three directions. First, it became part of the anisotropic SGWB workflow used with Advanced LIGO and Advanced Virgo folded data for broadband radiometer, narrow-band radiometer, and spherical-harmonic decomposition searches over the first three observing runs [2103.08520]. Second, the algebra of pixel-basis and SpH mapmaking was unified within PyStoch, eliminating the need for manifestly redundant pipelines and making per-frequency SpH searches computationally feasible [2011.05969]. Third, the same folded cross-correlation machinery was repurposed for directed CW searches, where PyStoch serves as a low-cost first-pass filter for nearly monochromatic signals from known sky positions, with more sensitive CW pipelines reserved for follow-up [2507.15027].

A common misconception is that PyStoch is only a broadband stochastic-background mapmaker. The published literature instead describes it as a general directional cross-correlation engine that produces maps at every frequency bin, supports both pixel and SpH representations, and has been used in both stochastic-background and directed CW searches [2110.09834].

## 2. Statistical estimator and sky representations

PyStoch implements the standard radiometer maximum-likelihood mapmaker. Let \(s_1(t;f)\) and \(s_2(t;f)\) denote Fourier components of the strain in two detectors for baseline \(I\). In the pixel basis, the “dirty” map \(X_p\) and Fisher matrix \(\Gamma_{pp'}\) are
\[
X_p  =  \frac{4}{\tau}  \sum_{I,f,t} \frac{ H(f)\,\gamma^*_{I f t, p} }{P_{I_1}(t;f) P_{I_2}(t;f)}\,s^*_{I_1}(t;f)\,s_{I_2}(t;f),
\]
\[
\Gamma_{pp'}  =  4 \sum_{I,f,t} \frac{H^2(f)}{P_{I_1} P_{I_2}}\,\gamma^*_{I f t, p}\,\gamma_{I f t, p'}.
\]
Here \(\tau\) is the total observation time, \(H(f)\) is the assumed SGWB frequency spectrum, \(P_{I_1,2}(t;f)\) are the one-sided noise PSDs, and \(\gamma_{Ift,p}\) is the overlap reduction function (ORF) for baseline \(I\) and sky pixel \(p\) [2011.05969]. The maximum-likelihood estimate in pixel space is
\[
\hat d_p = \sum_{p'} (\Gamma^{-1})_{pp'} X_{p'}.
\]

In the SpH basis, the anisotropy is expanded as
\[
\mathcal{P}(\Omega) = \sum_{\ell m} a_{\ell m} Y_{\ell m}(\Omega),
\]
with dirty coefficients and Fisher matrix
\[
X_{\ell m}  =  \int_{S^2} d\Omega\,Y^*_{\ell m}(\Omega)\,X(\Omega),
\qquad
\Gamma_{\ell m,\ell'm'} = 4 \sum_{I,f,t} \frac{H^2}{P P}\,\gamma^*_{I f t,\ell m}\,\gamma_{I f t,\ell'm'}.
\]
The corresponding maximum-likelihood estimator is
\[
\hat a_{\ell m} = \sum_{\ell'm'} (\Gamma^{-1})_{\ell m,\ell'm'} X_{\ell'm'}.
\]
This formulation supports both localized and diffuse searches: pixel maps are used to probe pointlike hot spots, while SpH moments constrain diffuse anisotropies [2011.05969].

The same formalism also underlies isotropic searches. In the unified PyStoch treatment, the isotropic limit is recovered by setting \(\ell_{\max}=0\), which yields the standard isotropic SGWB estimator [2011.05969].

## 3. Folded-data architecture and computational organization

The key algorithmic acceleration in PyStoch is sidereal folding. Because antenna patterns, time-delay phase factors, and related baseline-dependent quantities repeat every sidereal day, months of data can be co-added into a single “folded day.” In the all-sky, all-frequency (ASAF) implementation, ordinary analyses step through segments with \(\Delta t=192\,\mathrm{s}\), whereas folding reduces the time axis to \(N_{\mathrm{fold}}=T_s/\Delta t \approx 450\) sidereal-time slices; subsequent ORF evaluation, noise weighting, and mapmaking are then performed only on the folded data [2110.09834].

PyStoch expresses the mapmaking operations as linear algebra. In the original folded-data implementation, the baseline-dependent sky response is compressed into two “seed” matrices: an antenna-pattern product matrix and a time-delay seed matrix. For each frequency, the full ORF is generated on the fly by multiplying the antenna-pattern seed by the phase factor \(\exp[2\pi i f\,(\Delta x\!\cdot\!\Omega/c)]\), so the full time-frequency-pixel ORF table never needs to be stored [1803.08285]. This compression is central to the claim that PyStoch can stream over frequency while keeping memory requirements modest.

HEALPix is integrated throughout the pipeline. PyStoch stores sky-dependent quantities as HEALPix maps, uses standard angular-to-pixel conversions, and exploits HEALPix spherical-harmonic transforms. In the unified anisotropy pipeline, the antenna-pattern product \(F_{I_1}F_{I_2}(\Omega,t)\) and the time delay \(\Delta x\!\cdot\!\Omega/c\) are stored as frequency-independent “seed” maps for each folded time segment, while conversion to the SpH basis uses HEALPix’s fast `map2alm` routines; because the ORFs are complex, a thin wrapper around `map2alm` was added [2011.05969]. In the ASAF search, PyStoch further factors the ORF into a “time-dependent modulation” tensor and a “frequency-dependent phase” tensor, avoiding a naïve nested loop over \((\tau,p,f)\) [2110.09834].

The sky resolution depends on the analysis. The ASAF search used HEALPix \(N_{\rm side}=16\), corresponding to \(N_{\rm pix}=3072\) equal-area pixels, while the O1–O3 broadband radiometer analysis used \(n_{\rm side}=32\) [2110.09834; 2103.08520].

## 4. Unified mapmaking, basis transformations, and multicomponent inference

A central technical result associated with PyStoch is that pixel-basis and SpH estimators are algebraically transformable. The map coefficients satisfy
\[
a_{\ell m} = \sum_p Y^*_{\ell m}(\Omega_p)\,d_p,
\qquad
d_p = \sum_{\ell m} Y_{\ell m}(\Omega_p)\,a_{\ell m},
\]
and the Fisher matrices transform analogously,
\[
\Gamma_{\ell m,\ell'm'} = \sum_{p,p'} Y^*_{\ell m}(\Omega_p)\,Y_{\ell'm'}(\Omega_{p'})\,\Gamma_{pp'}.
\]
This relation implies that a pixel-basis analysis can be projected into the SpH basis after the fact, or the SpH Fisher matrix can be built directly from \(\gamma_{Ift,\ell m}\) [2011.05969].

Earth’s rotation introduces an especially useful simplification:
\[
\gamma_{I f t,\ell m} = \gamma_{I f 0,\ell m}\,e^{i m 2\pi t/T}.
\]
As a result, only the \(t=0\) transform needs to be precomputed, with the remaining time dependence obtained by a phase factor. In the reported implementation, this reduces the cost enough that per-frequency ORFs in the SpH basis are generated in \(\lesssim 10\,\mathrm{ms}\), and pixel and SpH maps with full Fisher matrices for one observing run can be built in \(\lesssim 2\,\mathrm{hr}\) on a modern multi-core CPU [2011.05969].

The same output objects support multicomponent anisotropy inference. For multiple assumed spectral shapes \(H_\alpha(f)\), PyStoch constructs dirty maps \(X_u^\alpha\) and Fisher blocks \(\Gamma^{\alpha\beta}_{uu'}\). The formal joint clean-map solution is
\[
\widehat{\mathcal P}^{\alpha}_{u}
=
\sum_{\beta,u'} [C^{-1}]^{\alpha\beta}_{uu'}\,X^\beta_{u'},
\]
with covariance \([C^{-1}]^{\alpha\beta}_{uu'}\), where \(C^{\alpha\beta}_{uu'}\equiv\Gamma^{\alpha\beta}_{uu'}\) [2106.09593]. At present detector sensitivities, the off-diagonal pixel-pixel terms are reported to be very small compared to the diagonal, and the standard analyses therefore drop them and invert only the much smaller \(N\times N\) matrices at each pixel. The published argument is that this changes neither upper limits nor point-spread recovery at any level significant for O1–O3, while making joint-index anisotropic mapmaking tractable [2106.09593].

This unified architecture also underlies later downstream products. The literature explicitly notes clean maps via SVD or norm regularization, decorrelated \(C_\ell\) spectra, polarization studies, and lensing as analyses that can be built on the common PyStoch output [2011.05969].

## 5. Role in anisotropic SGWB searches

PyStoch has been used in the major directional SGWB searches with Advanced LIGO and Advanced Virgo data. In the O1–O3 anisotropic search, the pipeline was run on data folded over one sidereal day and used for broadband radiometer, narrow-band radiometer, and spherical-harmonic decomposition analyses. No evidence for gravitational-wave signals was found. The reported \(95\%\) confidence-level upper limit sky maps on the gravitational-wave energy flux from broadband point sources ranged from \(F_{\alpha,\Theta} < (0.013 - 7.6)\times10^{-8}\,\mathrm{erg\,cm^{-2}\,s^{-1}\,Hz^{-1}}\), while upper limits on the normalized gravitational-wave energy density spectrum from extended sources ranged from \(\Omega_{\alpha,\Theta} < (0.57 - 9.3)\times10^{-9}\,\mathrm{sr^{-1}}\), depending on direction and spectral index; these limits improved on previous limits by factors of \(2.9 - 3.5\) [2103.08520].

In that analysis, the radiometer output was translated into physically interpretable quantities through
\[
\hat{\mathcal F}_{\alpha,\Theta} = \frac{c^3\pi}{4G}\,f_{\rm ref}^2\,\hat{\mathcal P}_{\alpha,\Theta},
\qquad
\hat{\Omega}_{\alpha,\Theta} = \frac{2\pi^2}{3H_0^2}\,f_{\rm ref}^3\,\hat{\mathcal P}_{\alpha,\Theta},
\]
with the first expression used for point-source broadband radiometer maps and the second for extended-source SHD maps [2103.08520].

PyStoch was also the enabling pipeline for the first all-sky all-frequency anisotropic SGWB search. That ASAF search analyzed every \(1/32\,\mathrm{Hz}\) frequency bin from \(20\) to \(1726\,\mathrm{Hz}\), except bins likely to contain instrumental artefacts and hence notched, over \(3072\) HEALPix sky pixels. No statistically significant evidence for narrowband gravitational-wave signals was found, and \(95\%\) confidence upper limits on the gravitational-wave strain were set for each pixel-frequency pair, in the range \((0.030 - 9.6)\times10^{-24}\) [2110.09834]. The same work stressed that a broadband analysis can miss narrowband signals because the signal-to-noise ratio of a narrowband signal can be significantly reduced when combined with detector output from other frequencies; this is the main methodological motivation for per-frequency mapmaking [2110.09834].

Candidate identification in ASAF used a noise-background estimation based on unphysical time shifts and a per-frequency maximum-SNR scan over contiguous \(10\,\mathrm{Hz}\) sub-bands. Hardware injections in O2 were recovered as outliers by this procedure, providing an internal validation of the follow-up logic [2110.09834].

## 6. Directed continuous-wave searches

PyStoch has also been repurposed for directed CW searches from known sky locations, where it functions as a radiometer-style coherent cross-correlation statistic over folded data. In the O3 study targeting four supernova remnants—Vela Jr., G347.3-0.5, Cassiopeia A, and the neutron star associated with the 1987A supernova remnant—the pipeline scanned the \(20\)–\(1726\,\mathrm{Hz}\) band with default settings \(T_{\rm coh}=192\,\mathrm{s}\) and \(\delta f_{\rm def}=1/32\,\mathrm{Hz}\). No interesting candidates were identified in the real O3 analysis, and \(95\%\) confidence-level upper limits on the CW strain amplitude \(h_0\) were set. The most stringent limit was \(h_0 = 1.13\times10^{-25}\) at \(201.57\,\mathrm{Hz}\) for Cassiopeia A; the best limits for the other targets were \(1.20\times10^{-25}\) at \(202.16\,\mathrm{Hz}\) for G347.3-0.5, \(1.20\times10^{-25}\) at \(217.81\,\mathrm{Hz}\) for Vela Jr., and \(1.47\times10^{-25}\) at \(186.41\,\mathrm{Hz}\) for the 1987A target [2507.15027].

The O3 CW implementation used a “bincomb” strategy to cope with spin-down broadening: adjacent \(1/32\,\mathrm{Hz}\) bins were combined in a sliding window of width \(\delta f_{\rm comb}=(2N+1)\delta f_{\rm def}\), with candidates flagged when the empirical \(p\)-value was \(p\le0.10\), roughly \(\mathrm{SNR}>4.5\). Follow-up pipelines were 5-vector Resampling and Band-Sampled Data Frequency-Hough [2507.15027].

In the O4a supernova-remnant search, PyStoch was one of five directed pipelines and was described as a cross-correlation-based method that uses \(T_{\rm coh, sid}=192\,\mathrm{s}\) and \(\delta f=1/32\,\mathrm{Hz}\) folded Stochastic Intermediate Data. The O4a workflow included target-specific sky demodulation, adaptive bin grouping over window sizes from 1 bin up to an empirical maximum of 48 bins (\(1.5\,\mathrm{Hz}\)), and empirical \(p\)-value thresholding with \(p\le0.10\), equivalently \(\rho\gtrsim4.5\). In that study, none of the PyStoch candidates survived coherence follow-up. The tightest O4 Pan-PyStoch constraint was \(h_0^{95\%}\simeq 9\times10^{-26}\) near \(320\,\mathrm{Hz}\) for Vela Jr., corresponding to a sensitivity depth \(D\equiv\sqrt{S_h}/h_0\approx36\,\mathrm{Hz^{-1/2}}\) [2603.25808].

These CW applications clarify PyStoch’s niche. The O3 supernova-remnant study reported that, relative to standard CW directed searches, PyStoch is approximately \(2\)–\(3\) times less sensitive but approximately \(10\)–\(20\) times faster, making it suitable as a first-pass filter for targets with poorly constrained spin parameters [2507.15027].

## 7. Validation, scaling behavior, and limitations

Validation results reported for PyStoch are unusually explicit. In the unified pixel/SpH mapmaking study, SpH moments obtained by direct cross-spectral-density projection and by pixel-to-SpH transformation matched to \(\sim0.1\%\) rms for \(\ell\le30\), while Fisher matrices agreed to \(\lesssim0.01\%\) between the conventional and PyStoch implementations. On that basis, the study concluded that a single skymap is sufficient to describe anisotropies in a stochastic background [2011.05969].

Performance claims are analysis-dependent but consistently large. In the original folded-data mapmaking paper, a representative O1-like benchmark was reduced from \(\sim10\) CPU-years for the conventional unfolded pipeline to \(\sim10\) CPU-days with folding alone and to \(\sim40\) CPU-minutes with folding plus PyStoch; intermediate storage also shrank from \(\sim450\,\mathrm{GB}\) to a \(\sim1.5\,\mathrm{GB}\) folded dataset [1803.08285]. In the ASAF search, PyStoch plus folding was reported to run \(\sim400\)–\(600\times\) faster than the original Matlab implementation without folding; each baseline-plus-run required \(\lesssim1\,\mathrm{hour}\) on a single modern 2-core laptop, and the full O1+O2+O3 HLV combination took \(O(3\)–\(4\,\mathrm{h})\) wall time on one core [2110.09834]. In the O3 CW supernova-remnant study, the PyStoch stage ran in \(\sim0.5\,\mathrm{hr}\) on a 4-core CPU [2507.15027].

The literature also records concrete limitations. In the 2025 CW implementation, PyStoch had no built-in demodulation for Doppler or spin-down, relied on the bincomb approximation, used fixed \(T_{\rm coh}=192\,\mathrm{s}\), and folded data precluded segment-by-segment vetoes of non-Gaussian artifacts [2507.15027]. In anisotropic multicomponent inference, the approximation of neglecting off-diagonal pixel-pixel covariance is justified only because present detector sensitivities render those terms very small relative to the diagonal [2106.09593]. These are not changes to the basic cross-correlation estimator, but they define the practical envelope within which PyStoch has been deployed.

Taken together, the published record places PyStoch at the intersection of radiometry, folded-data compression, and fast spherical analysis on the two-sphere. Its defining feature is not a new detection statistic, but an implementation strategy that makes frequency-resolved directional inference routine enough to support broadband SGWB searches, all-sky all-frequency scans, multicomponent anisotropy studies, and computationally economical CW candidate generation within a single family of workflows [1803.08285].

Source: https://www.emergentmind.com/topics/pystoch