Bayesian nonparametric estimation of correlated gravitational wave detector network noise using matrix-gamma process priors

Technical Summary: Bayesian Nonparametric Estimation of Correlated Gravitational Wave Detector Network Noise

Problem Statement
Accurate characterization of noise spectral density is critical for next-generation gravitational-wave (GW) detectors like the Laser Interferometer Space Antenna (LISA) and the Einstein Telescope (ET). Unlike current ground-based networks (e.g., LVK) where channel independence is a reasonable assumption due to spatial separation, LISA and ET operate in regimes where cross-channel correlations are physically expected and significant. These correlations arise from shared noise sources (e.g., magnetic, seismic, Newtonian noise) and instrumental coupling (e.g., test-mass and optical-metrology noise).

Ignoring these correlations or mis-specifying the Power Spectral Density (PSD) can lead to biased GW signal parameter estimation, incorrect model selection, and miscalibrated uncertainties. Existing methods face a trade-off:

• Parametric models are computationally efficient but may miss localized spectral features or complex cross-spectral structures.
• Nonparametric methods (e.g., Welch estimates) are flexible but often provide single fixed estimates without propagating uncertainty, or they fail to guarantee that the estimated spectral density matrix remains Hermitian positive definite (a requirement for likelihood inversion).
• Variational approximations scale well but tend to underestimate posterior uncertainty.

The paper addresses the need for a method that is nonparametric over the full spectral density matrix, guarantees positive definiteness by construction, and provides honest posterior uncertainty for propagation into signal inference.

Methodology
The authors propose a Bayesian nonparametric framework, termed VNP (Vector Nonparametric), and a semi-parametric extension, VNP-P (Vector Nonparametrically corrected Parametric).

1. Likelihood Construction:

   • To handle the extremely long time series generated by these detectors, the data is segmented into non-overlapping blocks.
   • The authors employ a blocked Whittle likelihood in complex-Wishart form. By averaging periodogram matrices across blocks, the likelihood is approximated as a product of independent complex Wishart distributions. This reduces computational complexity while retaining the full cross-spectral structure without requiring a low-dimensional parametric form.

2. Prior Specification (VNP):

   • The core innovation is a matrix-gamma process prior placed on the coefficients of a Bernstein polynomial basis expansion.
   • This generalizes the univariate Bernstein-Dirichlet process to the multivariate setting. The spectral density matrix $S(f)$ is modeled as a mixture of Beta densities (Bernstein polynomials) with matrix-valued weights.
   • The matrix-gamma process ensures that the resulting spectral density estimate is Hermitian positive definite at every frequency, a crucial property for inverting the matrix within the likelihood.
   • The model avoids reversible-jump MCMC by representing the matrix-valued basis coefficients as values of a cumulative distribution function, sampled via an adaptive Markov chain Monte Carlo (MCMC) algorithm.

3. Semi-Parametric Extension (VNP-P):

   • To improve efficiency when a rough parametric model is available, the authors introduce a two-step approach.
   • First, a parametric model (specifically a Vector Autoregressive model, VAR) is fitted to capture the broad spectral shape.
   • Second, the full spectral density is modeled as the matrix geometric mean of the parametric estimate ($S_p$) and a nonparametric correction term ($S_{np}$).
   • The matrix-gamma process prior is placed on the correction factor $C = (S_p^{-1/2} S_{np} S_p^{-1/2})^{1/2}$ rather than the full spectrum. This allows the nonparametric component to focus on modeling residual discrepancies (e.g., sharp peaks or misspecified correlations) while the parametric component handles the bulk structure.

Key Results
The paper validates the methods through simulation studies and applications to simulated LISA and ET data:

• Simulation Studies (VAR and VMA processes):

  • For data generated from a VAR(2) process, the VNP-P method (using a VAR(1) working model) achieved lower $L_2$ error than the pure VNP method, demonstrating that a well-specified parametric "working" model can improve nonparametric estimation accuracy.
  • For data generated from a VMA(1) process, where the VAR(1) working model was misspecified, VNP-P offered no accuracy gain over VNP, confirming that the benefits depend on the working model's ability to capture dominant features.
  • The methods demonstrated robust posterior consistency, with credible intervals showing reasonable frequentist coverage, though coverage was slightly lower than nominal for finite block lengths in the VAR(2) case.

• LISA Application:

  • Applied to simulated second-generation TDI channel noise, the VNP method successfully recovered the true PSD and coherence estimates across the frequency band.
  • The method handled the complex transfer functions and frequency-dependent dips characteristic of LISA data, though slight underestimation at low frequencies was observed, attributed to spectral leakage.

• Einstein Telescope (ET) Application:

  • Simulated ET data included correlated Gaussian noise peaks injected into X, Y, and Z channels.
  • The pure VNP method oversmoothed the correlated peaks, failing to resolve them clearly.
  • The VNP-P method, guided by a high-order VAR(303) working model, successfully captured the correlated peaks and provided a better fit (lower $L_2$ error: 0.068 vs. 0.078 for VNP).
  • This highlights the utility of VNP-P in scenarios where the noise contains specific, correlated features that a parametric model can partially approximate.

Significance and Claims
The paper positions this work as a benchmark for gravitational-wave noise estimation. The authors claim:

1. Theoretical Guarantee: The proposed method is the first to provide a nonparametric estimate of the full spectral density matrix that is guaranteed to be Hermitian positive definite by construction, without relying on reversible-jump MCMC or variational approximations that may underestimate uncertainty.
2. Uncertainty Propagation: Unlike fixed PSD estimates (e.g., Welch), this approach provides honest posterior uncertainty, which is essential for propagating noise uncertainty into GW signal parameter inference.
3. Reference Standard: Due to its computational intensity (MCMC-based), the method is not intended to replace faster approximate methods (like stochastic gradient variational Bayes) for routine analysis. Instead, it serves as a reference standard against which the accuracy and uncertainty quantification of faster, approximate methods can be evaluated.
4. Flexibility: The framework is flexible enough to correct misspecified parametric noise models, offering a pathway to improve estimation in complex, correlated noise environments typical of next-generation detectors.

The authors conclude that while the method is computationally demanding, it provides a necessary tool for validating faster approximations and ensuring that correlated noise in LISA and ET does not bias scientific results. Future work is suggested to integrate this into joint signal-plus-noise (on-source) estimation frameworks and to apply it to more realistic noise scenarios involving gaps and glitches.