Opening the Black Box: What Neural Networks Learn from Pulsar Timing Array Data

The gist The work uses a simple toy model of a PTA to provide a more pedagogical explanation of how, and under which circumstances, SBI-based detection can improve on classical detection statistics for gravitational-wave background anisotropies.

PTA Model and Signal Description

In the 2D world considered for analytical calculations, the Earth lies at the origin and all pulsars lie on a unit circle centered on the origin (Page 2). A generic metric perturbation is given by a superposition of plane waves (Page 2). A Gaussian and stationary GWB is completely described in terms of the one- and two-point functions of the Fourier modes, h˜(f, ϕ) (Page 2). The sky map P(ϕ) describes the distribution of GWB power on the sky, normalized such that R 2π0dϕ P(ϕ) = 1 (Page 2).

Timing Residuals and Overlap Reduction Function

The timing residuals are given by δta(t), which is related to the Fourier transform of h˜(f, ϕ) Fa(ϕ) (Page 2). The two-point function for the stochastic signal produced by the GWB is 1/2δ(f − f′)H(f)Γab, where Γab is the overlap reduction function (ORF), which in general depends on the distribution of GW power on the sky (Page 3). For an isotropic GWB, this yields a 2D-analog of the Hellings-Downs (HD) correlation function Γiso ab = 1/2cos(γab), where γab = ϕa − ϕb (Page 3). For anisotropic GWBs, the ORF is decomposed as Γab = Γiso ab + Γaniso ab, where Γaniso ab = -πP r2cos(ϕa + ϕb) − P i2sin(ϕa + ϕb) (Page 3).

Likelihood Function and Covariance Matrix

The likelihood function for the timing residuals is given by exp-1/2˜δt†·C−1(θ)·˜δt / det(2πC(θ)) (Page 4). The covariance matrix C(θ) depends on the model parameters θ for either the isotropic or anisotropic signal models (Page 4). Under simplifying assumptions, the covariance matrix for the Fourier series coefficients is approximately diagonal: Caj,bk ≈ Ns/δjk δabσ2n + Γabσ2GWB (Page 4). This leads to a sufficient statistic ρab ≡ XNf j=−Nf˜δtaj ˜δt∗bj (Page 4).

Detection Statistics and Neural Networks

Several detection statistics are compared, including the Bayes factor, the sky-map-based frequentist statistic, and linear and non-linear neural-network classifiers (Page 5). For simple hypotheses, the log of the Bayes factor reduces to a linear function of the correlations ρ or equivalently quadratic in the timing residual data δt (Page 5). The BCE minimization drives a flexible classifier toward the log Bayes factor, and for simple hypotheses, both linear and non-linear classifiers can approximate it exactly (Page 8). For composite hypotheses where anisotropy parameters are unknown, only the MLP can efficiently approximate the log Bayes factor because the linear network performs no better than chance (Page 9).

Conclusion on Neural Network Learning

The training of a classifier is shown to amount to constructing a frequentist detection statistic that, given sufficient network flexibility, reaches the sensitivity of a full Bayesian analysis at a fraction of its cost (Page 10). The linear network learns nothing from the timing residual data because it is only capable of building statistics linear in the input data (Page 11). However, when trained on quadratic data like correlations ρab with non-zero mean priors on anisotropy parameters, the linear network can learn non-trivial weights and obtain a quadratic approximation to the log Bayes factor (Page 11). The non-linear classifier tracks the log Bayes factor up to values of log B10 ∼ 2 before it starts to become under-confident (Page 8).

How it works

The paper introduces several detection statistics for comparison:

1. Bayes factor: defined as the ratio of evidences B10(δt) = p(δtH1) / p(δtH0) (Page 5).

2. Sky-map-based detection statistics: constructed from an estimator of the GWB sky map P, which is maximized using a Gaussian likelihood (Page 5).

3. Neural network detection statistics: trained by minimizing the binary cross-entropy (BCE) functional, which drives the network output toward the log Bayes factor (Page 6).

The neural network output fλ(ρ) is defined as σzλ(ρ), where zλ(ρ) is the pre-activation output of the final neuron and σ is the sigmoid activation function (Page 6). Minimizing BCE leads to z⋆(ρ) = ln B10(ρ), which means training drives the network toward the optimal statistic (Page 6).

The linear network learns nothing from timing residuals because a linear network is only capable of building statistics linear in the input data, while the information regarding the GWB is encoded in quadratic combinations (i.e., the correlations) of the timing residual data (Page 11). Conversely, when trained on correlations ρab with non-zero mean priors on P r2 and P i2, a linear neural network can learn non-trivial weights under training and obtain a quadratic approximation to the log Bayes factor (Page 11). The non-linear MLP is able to match the log Bayes factor across the full range of false alarm probabilities when marginalizing over unknown anisotropy parameters (Page 9).

Results Summary

For simple hypotheses, both linear and non-linear classifiers perform identically, reaching the Bayes factor exactly because it is a linear function of correlations (Page 8). In the weak-signal limit for simple hypotheses, the learned weights of the linear network trace the expected cos(ϕa + ϕb + π/4) profile (Page 8). For composite hypotheses, only the MLP can efficiently approximate the log Bayes factor because averaging over zero-mean priors on anisotropy parameters makes any statistic linear in ρab indistinguishable (Page 9). The map-based statistic performs better than the linear classifier but still falls short of the optimum due to treating cross-correlations as Gaussian variables (Page 9). The non-linear classifier tracks the log Bayes factor up to values of log B10 ∼ 2 before it starts to become under-confident (Page 8).

Conclusion

Training such a classifier amounts to constructing a frequentist detection statistic that, given sufficient network flexibility, reaches the sensitivity of a full Bayesian analysis at a fraction of its cost (Page 10). The work demonstrates that training such a classifier amounts to constructing a frequentist detection statistic that, given sufficient network flexibility, reaches the sensitivity of a full Bayesian analysis at a fraction of its cost (Page 10). This pedagogical work has used a simple PTA toy model as a case study to understand what neural classifiers learn when used as a detection statistic for GWB anisotropies (Page 10). The paper found that training such a classifier amounts to constructing a frequentist detection statistic that, given sufficient network flexibility, reaches the sensitivity of a full Bayesian analysis at a fraction of its cost (Page 10). The work establishes that both classifiers recover the optimal statistic for simple hypotheses and in the linear case they do so by learning the Neyman-Pearson weights themselves (Page 8). The overall conclusion is that training such a classifier amounts to constructing a frequentist detection statistic that, given sufficient network flexibility, reaches the sensitivity of a full Bayesian analysis at a fraction of its cost (Page 10). The work provides intuition for what SBI is doing and how it does it (Page 10).

--- Page 1 ---

Imperial–TP–2026–AM–04 James Alvey, Andrea Mitridate, Mauro Pieroni and Joseph D. Romano (Page 1) In recent years, simulation-based inference (SBI) methods have been proposed to address several data-analysis challenges faced by existing and planned gravitational-wave experiments (Page 1). Over the past few years, pulsar timing array (PTA) collaborations have reported mounting evidence for a nanohertz gravitational-wave background (GWB), whose characteristic signature is the Hellings-Downs pattern of correlations between the timing residuals of pairs of pulsars [7–13] (Page 1). In this work, we use a simple toy model of a PTA to provide a more pedagogical explanation of how, and under which circumstances, SBI-based detection can improve on classical detection statistics for gravitational-wave background anisotropies (Page 1).

--- Page 2 ---

2 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, United Kingdom (Page 2) A. Metric perturbations and two-point functions (Page 2) For this 2D example, a generic metric perturbation is given by a vector⃗h(t,⃗x), which can be decomposed as a superposition of plane waves (Page 2). A Gaussian and stationary GWB is completely described in terms of the one- and two-point functions of the Fourier modes, h˜(f, ϕ) (Page 2). We will assume that the one-point function vanishes, i.e., ⟨h˜(f, ϕ)⟩ = 0, where ⟨·⟩ denotes an average over the Gaussian ensemble (Page 2).

1. Bold header: Learn Optimal Statistics via Neural Networks

Neural networks can be trained to minimize a cost function derived from binary cross-entropy, which drives the network output toward the log Bayes factor. This allows them to recover the optimal statistic without evaluating a single evidence integral, effectively matching the sensitivity of a full Bayesian analysis at a fraction of its computational cost.

2. Bold header: Recover Anisotropy Information in Composite Hypotheses

The paper demonstrates that while linear networks underperform the Bayes factor when marginalizing over anisotropy parameters, the non-linear MLP should be able to match it. This capability allows the AI system to correctly distinguish between null and signal hypotheses even when model parameters like anisotropy are unknown.

3. Bold header: Learn Optimal Weights for Simple Hypotheses

For simple hypotheses, the linear network can learn the Neyman-Pearson weights themselves, as its learned weights trace the expected cos(ϕa + ϕb + π/4) profile. This provides a direct, analytically verifiable comparison between the learned statistic and the optimal analytic prediction.

4. Bold header: Estimate Sky Map Parameters from Cross-Correlations

The map-based detection statistic can be used to find estimators for anisotropy parameters like Pˆm = M−1 m Xm, which are then used to construct the detection statistic Λmap = Pˆ22 + Pˆ−22. This allows the AI system to estimate sky map components even when they are not known a priori.

In recent years, simulation-based inference (SBI) methods have been proposed to address several data-analysis challenges faced by existing and planned gravitational-wave experiments. For example, SBI classification has recently been shown to significantly improve the prospects for detecting anisotropies in pulsar timing array (PTA) data. In this work, we use a simple toy model of a PTA to provide a more pedagogical explanation of how, and under which circumstances, SBI-based detection can improve on classical detection statistics for gravitational-wave background anisotropies.