Bayesian Deep Gaussian Processes for Correlated Functional Data: A Case Study in Cosmological Matter Power Spectra

Technical Summary: Bayesian Deep Gaussian Processes for Correlated Functional Data

Problem Statement
The paper addresses the challenge of estimating the underlying "infinite volume" matter power spectrum, $P_\infty(k)$, from cosmological computer simulations. Specifically, the authors focus on the Mira-Titan simulation suite, which provides multiple response curves of varying fidelities for a given set of cosmological parameters (a "cosmology"). These outputs include:

1. Perturbation theory spectra ($P_p$): Analytically derived, accurate only at large spatial scales (low wavenumbers, $k < 0.04$).
2. Low-resolution simulations ($P_{\ell r}$): 16 correlated realizations per cosmology, accurate for intermediate scales ($0.04 \le k \le 0.25$) but subject to initialization randomness.
3. High-resolution simulations ($P_h$): A single realization per cosmology, accurate for a broader range ($0.04 \le k \le 5$) but still subject to finite-volume effects at low $k$.

The core difficulty lies in synthesizing these heterogeneous, correlated, and noisy functional outputs to estimate the smooth, non-random $P_\infty(k)$ while providing rigorous Uncertainty Quantification (UQ). Existing methods often treat these outputs as independent scalar observations or fail to account for the spatial dependence inherent in the functional curves across wavenumbers. Furthermore, the authors aim to build a surrogate model (emulator) to predict spectra for unobserved cosmologies without running expensive new simulations.

Methodology
The authors propose a Bayesian Hierarchical Deep Gaussian Process (DGP) model, termed DGP.FCO, designed specifically for correlated functional outputs.

1. Hierarchical Structure:

   • Observation Layer: The observed weighted average of the simulation outputs ($\bar{Y}$) is modeled as a noisy realization of the underlying spectrum $S$. Unlike previous approaches that assume diagonal noise covariance, this model employs a dense covariance matrix $\Sigma_\epsilon$ (based on a Matérn kernel) to capture the smooth spatial dependence and correlations across wavenumbers.
   • Latent Spectrum Layer: The underlying spectrum $S$ is modeled as a Gaussian Process centered at zero, $S \sim \text{GP}(0, \Sigma_S(W))$.
   • Deep Latent Layer: To handle non-stationarity (specifically baryonic acoustic oscillations in the range $-1.4 \le \log_{10}(k) \le -0.8$), the input space $X$ (wavenumbers) is warped by a latent monotonic Gaussian Process, $W \sim \text{monoGP}$. This warping allows the model to adapt the smoothness of the spectrum realization across different scales.

2. Inference:

   • The authors use Elliptical Slice Sampling (ESS) to infer the latent warping function $W$ and Metropolis-Hastings steps within a Gibbs framework to update kernel hyperparameters ($\theta_S, \theta_W$).
   • The posterior distribution of the spectrum $S$ is derived in closed form given the latent variables, allowing for direct sampling to compute posterior means and 95% credible intervals.

3. Emulation for Unobserved Cosmologies:

   • To predict spectra for new cosmologies, the authors employ a two-stage approach. First, the DGP.FCO model is used to estimate the spectra for the training cosmologies.
   • These estimated spectra are then decomposed using Principal Component Analysis (PCA).
   • A separate Gaussian Process emulator is trained on the principal component weights as a function of the eight cosmological parameters. This surrogate model allows for rapid prediction of spectra for unobserved cosmological configurations.

Key Contributions

• Extension of Bayesian DGPs: The paper extends Bayesian Deep Gaussian Processes from scalar responses to correlated functional outputs. This is a novel application of DGPs in the context of computer experiments.
• Synthesis of Multi-Fidelity Data: The model effectively synthesizes perturbation theory, low-resolution, and high-resolution simulation data, respecting the specific wavenumber ranges where each data source is unbiased.
• Incorporation of Domain Knowledge: The model explicitly incorporates expert knowledge regarding the reliability of different data sources (via precision matrices) and the non-stationary nature of the spectra (via the monotonic latent GP).
• Benchmarking: The method is validated against synthetic data and the CAMB (Code for Anisotropies in the Microwave Background) model, which provides a known "infinite volume" ground truth. It is also benchmarked against the state-of-the-art Cosmic Emu (a deep process convolution approach) on the Mira-Titan dataset.

Results

• Synthetic and CAMB Experiments: In simulation studies where the true underlying curve is known, the DGP.FCO model outperforms competing methods (standard DGP, heteroskedastic GP, and deep process convolution) in terms of log score, a proper scoring rule that evaluates both accuracy and uncertainty calibration. The authors note that accounting for spatial dependence between observations leads to more effective UQ.
• Mira-Titan Emulation: When predicting spectra for held-out cosmologies, the proposed method (DGP.FCO + PCA + GP emulator) performs favorably compared to the existing Cosmic Emu.
• Uncertainty Quantification: The hierarchical model successfully provides credible intervals that reflect the uncertainty arising from finite simulation volumes and initialization randomness, which simple weighted averages fail to capture adequately.

Significance and Claims
The authors claim that their approach provides a robust framework for estimating the infinite volume matter power spectrum from complex, multi-fidelity simulation suites. By moving beyond scalar responses to correlated functional outputs, the method captures the intrinsic smoothness and dependence of cosmological spectra that other emulators miss. The paper positions this work as a critical step toward efficient analysis of large-scale cosmological surveys (e.g., Euclid, LSST), where accurate and fast prediction of matter power spectra with reliable uncertainty quantification is essential for constraining cosmological parameters. The authors provide the data, reproducible code, and an R package to facilitate the adoption of this Bayesian DGP model.