BayesWave Initialization for GW Analysis

• BayesWave Initialization Algorithm is a suite of deterministic and stochastic procedures that sets high-quality starting conditions for gravitational-wave burst and glitch analyses.
• It employs fast PSD estimation and iterative Morlet–Gabor wavelet decomposition to distinguish coherent astrophysical signals from transient glitches in non-stationary noise.
• The algorithm reduces RJMCMC burn-in by seeding high-likelihood regions, leading to improved parameter estimation and real-time signal recovery.

The BayesWave Initialization Algorithm is a suite of deterministic and stochastic procedures, central to the Bayesian inference framework for gravitational-wave burst and glitch analysis, that accelerates convergence by providing a high-quality starting point for signal and noise model parameter estimation. Designed to address non-stationary, non-Gaussian noise environments inherent to interferometric gravitational-wave detectors (such as LIGO and Virgo), the algorithm adaptively decomposes detector data into a parametrized sum of Morlet–Gabor (sine–Gaussian) wavelets for both coherent astrophysical signals and incoherent instrumental glitches. Central to its initialization philosophy is the idea of “modeling everything,” integrating robust power spectral density (PSD) estimation (BayesLine or more recent techniques such as FastSpec) and iterative denoising to create favorable initial conditions for the trans-dimensional Reversible Jump Markov Chain Monte Carlo (RJMCMC) sampler.

1. Conceptual Structure and Motivation

BayesWave is structured around the additive model

$s = R \circledast h + g + n_G$

where $s$ is the detector data, $h$ is the gravitational wave signal, $g$ the glitch (non-Gaussian transient) component, $n_G$ stationary colored Gaussian noise, and $R$ the instrument response projection. Initialization aims to rapidly locate regions of high likelihood in the immense, variable-dimension parameter space by leveraging a deterministic approximation to the maximum likelihood for noise and glitch models, followed by rapid identification of excess signal power in the time–frequency domain.

The motivation is to minimize the computationally expensive "burn-in" period of the RJMCMC, efficiently “seeding” the sampler in a part of the parameter space where the posterior density is significant. Such initialization is crucial to accurately distinguish coherent astrophysical signals from incoherent glitches and stationary background noise, and to ensure robust model selection and parameter estimation, especially in non-Gaussian, transient-rich environments (Cornish et al., 2014, Cornish et al., 2020).

2. Fast PSD and Noise Model Initialization

The initialization begins with a deterministic PSD estimation, essential because a mischaracterized noise model can significantly bias subsequent Bayesian inference. The typical sequence is:

• Apply a window function (e.g., Tukey window) to data, then Fourier transform to obtain the raw power spectrum.
• Smooth the periodogram using a running median filter; the width is chosen to balance tracking the broadband PSD and suppressing narrow line artifacts.
• Identify spectral outliers (e.g., bins exceeding the running median by a factor—usually around 10)—to be modeled as Lorentzian lines.
• Fit the underlying broadband noise with a spline, originally a cubic spline but now typically an Akima spline for local smoothness and flexible knot placement (Gupta et al., 2023).
• Use a product combination for the spectral model: $S(f) = S_S(f)[1+S_L(f)]$, where $S_S(f)$ is the broadband spline component and $S_L(f)$ the normalized sum of Lorentzians, improving prior assignment and robustness.

This approach is further optimized in the FastSpec algorithm, which serves as a rapid, fixed-dimension Bayesian estimator for on-source data segments, dramatically reducing runtime compared to full RJMCMC-based PSD exploration (Gupta et al., 2023).

3. Iterative Glitch Model Seeding via Wavelet Decomposition

After whitening data using the PSD estimate, the algorithm proceeds to initialize the glitch (and potentially signal) model:

• Perform a wavelet transform using an overcomplete Morlet–Gabor basis, scanning the data in the time–frequency domain to identify significant excess power (regions where SNR for a single wavelet surpasses a threshold, typically SNR$>5$).
• For each substantial outlier, fit and subtract a wavelet from the data. Because Morlet–Gabor wavelets overlap, this step is iterative, recalculating the time–frequency map after each subtraction.
• Optionally, refine the fit in a maximum-likelihood fashion (e.g., via FFT-based maximization over amplitude, time, and phase) to further approach the global maximum for the glitch model parameters.
• Repeat the subtraction/fit process until no significant outliers remain.

Mathematically, for each wavelet the SNR estimate is

$$\rho^2 \approx \frac{A^2 Q}{2\sqrt{2\pi} f_0 S_n(f_0)},$$

where $A$ is amplitude, $Q$ quality factor, $f_0$ central frequency, and $S_n(f_0)$ the PSD at $f_0$ (Cornish et al., 2020).

This sequence yields an initial glitch model near the maximum likelihood, substantially accelerating subsequent convergence of the RJMCMC and reducing the need for numerous trans-dimensional (“birth”/“death”) moves in the early sampling phase.

4. Trans-Dimensional RJMCMC Initialization and Model Selection

Upon completion of deterministic seeding, the RJMCMC is initialized with this near-optimal configuration of the noise and glitch model. The RJMCMC operates over models where the number of wavelets (for both signal and glitch) is itself a parameter, enabling dimensionality to adapt during sampling (Cornish et al., 2014). Key elements include:

• The use of informed proposal distributions for wavelet parameters, often leveraging the local Fisher information matrix to sample efficiently in directions of steepest likelihood increase.
• Trans-dimensional “birth” and “death” moves allowing variable numbers of wavelets, with proposals informed by time–frequency maps and proximity or clustering heuristics to exploit signal coherence.
• Parallel tempering: Running multiple chains at different effective likelihood “temperatures” to ensure robust sampling of posterior regions and transition between modes.
• Careful selection of prior distributions for intrinsic parameters ($A$, $t_0$, $f_0$, $Q$, $\phi_0$), with SNR-based priors (e.g., $$p(\mathrm{SNR}) \propto \mathrm{SNR} \exp[-\mathrm{SNR}/\mathrm{SNR}_*]$$) enhancing efficiency and parsimony.

During burn-in, samples generated are not used for final inference but serve to push the chain rapidly toward global posterior maxima.

5. Signal Model Initialization: Coherence Versus Glitch Null Hypothesis

A defining feature of BayesWave's initialization is the simultaneous setup of models for both coherent signals and incoherent glitches:

• The coherent signal model imposes common intrinsic wavelet parameters across all detectors, modulo extrinsic parameters (source sky location, polarization, etc.), while the glitch model treats each detector’s decomposition independently.
• Initialization thus frames the competition between the hypotheses directly, setting up for subsequent computation of the signal-to-glitch Bayes factor

$$\mathcal{B}_{\mathcal{S},\mathcal{G}} = \frac{p(s|\mathcal{S})}{p(s|\mathcal{G})},$$

with analytic scaling derived as

$$\ln \mathcal{B}_{\mathcal{S},\mathcal{G}} \propto \mathcal{I} N \ln\,\mathrm{SNR_{net}},$$

where $\mathcal{I}$ is the number of detectors, $N$ the number of wavelets, and $\mathrm{SNR_{net}}$ the network SNR (Lee et al., 2021).

6. Algorithmic Enhancements: Modified Wavelet Basis, TF$\tau$ Transform, Heterodyning

Recent improvements focus on reducing computational bottlenecks in inner product calculations and overlapping wavelet subtractions:

• Modified wavelet bases are constructed such that redundant inner products $\langle \psi_i, \psi_j \rangle$ are precomputed or made analytic, accelerating model evaluation (Mathur et al., 18 Aug 2025).
• Transition from traditional Time–Frequency–Quality factor (TFQ) transforms to Time–Frequency–Time extent (TF$\tau$) transforms, reparameterizing the wavelet transform to optimize efficiency and minimize interpolation errors during subtraction.
• Downsampled heterodyned transforms shift the analysis band to baseband and reduce the data sampling rate, enabling rapid, low-latency reconstruction (e.g., in MaxWave), especially important for long-duration signals and machine learning pre-processing (Mathur et al., 18 Aug 2025).

7. Impact on Efficiency, Parameter Estimation, and Applications

The BayesWave initialization paradigm yields profound improvements in the practical application of Bayesian inference for gravitational-wave burst detection and characterization:

• Deterministic PSD and glitch model initialization shortens burn-in by orders of magnitude compared to random prior seeding, allowing rapid convergence of the full trans-dimensional sampler (Cornish et al., 2020, Gupta et al., 2023).
• Posterior exploration becomes highly adaptive, allocating model complexity where data warrant it—quiet data stretches use few wavelets, glitch-dominated or signal-rich stretches trigger additional wavelet components.
• The initialization procedures enhance overlap between injected and reconstructed waveforms, improve sky localization precision, and accelerate discrimination between coherent signals and non-Gaussian artifacts.
• These developments are foundational for real-time and offline analyses, including de-noising for stochastic background and continuous-wave searches, as well as for providing clean time–frequency representations to machine learning classifiers in the identification and understanding of glitches (Mathur et al., 18 Aug 2025).

In sum, the BayesWave Initialization Algorithm is both a deterministic and adaptive prelude to Bayesian inference for gravitational-wave data, ensuring that sophisticated noise, glitch, and signal models interact efficiently and accurately in the high-dimensional, non-linear search space encountered in realistic detector environments.