Vecchia approximated Bayesian heteroskedastic Gaussian processes

1. Problem Statement

Computer simulations in fields like ecology, epidemiology, and engineering are increasingly stochastic, meaning identical inputs can yield different outputs due to pseudo-random number generators or inherent system variability. A critical characteristic of these simulations is input-dependent noise (heteroskedasticity), where the variance of the output changes based on the input configuration.

Existing methods for modeling such data face two primary limitations:

1. Computational Intractability: Standard Heteroskedastic Gaussian Processes (hetGPs) require decomposing dense $N \times N$ covariance matrices, resulting in $O(N^3)$ computational complexity. This limits applicability to small datasets (typically $N < 1,000$).
2. Inferential Fidelity: To handle larger datasets, existing scalable methods (e.g., Binois et al., 2018) rely on point estimates (e.g., Maximum Likelihood Estimation or Expectation Maximization) for latent variances rather than full Bayesian posterior integration. This undercuts Uncertainty Quantification (UQ), as it ignores the uncertainty inherent in the variance estimates themselves.
3. Scale of Modern Campaigns: Recent applications, such as the NOAA-GLM lake temperature forecasting project, involve massive simulation campaigns with millions of data points ($N \approx 9$ million) and high replication, rendering standard approaches impossible.

2. Methodology

The authors propose bhetGP, a fully Bayesian framework that synthesizes three established techniques with a novel integration strategy to achieve scalability and full uncertainty quantification.

A. Core Components

1. Coupled Gaussian Processes: Following Goldberg et al. (1998), the model uses two coupled GPs: one for the mean function and another for the latent log-variance function.
2. Woodbury Identity (Replication): Leveraging Binois et al. (2018), the method utilizes sufficient statistics derived from replicated simulations. Instead of modeling $N$ total observations, it models $n$ unique input locations. This reduces the computational burden from $O(N^3)$ to $O(n^3)$, where $n \ll N$ in highly replicated campaigns.
3. Vecchia Approximation: To handle cases where the number of unique inputs ($n$) is still large (e.g., $n > 300,000$), the authors employ the Vecchia approximation. This induces sparsity in the covariance matrices by factorizing the joint likelihood into a product of low-dimensional conditional probabilities, reducing complexity to $O(nm^3)$, where $m$ is the number of conditioning neighbors.

B. Novel Integration: Elliptical Slice Sampling (ESS)

The paper's primary methodological innovation is the application of Elliptical Slice Sampling (ESS) for posterior inference of the latent variance field.

• Challenge: Traditional MCMC (Metropolis-Hastings) for high-dimensional latent variables often suffers from poor mixing and "sticky" chains.
• Solution: ESS is designed specifically for posteriors with Gaussian priors. It samples the latent variance vector $\Lambda_n$ by tracing an ellipse between the current state and a prior draw, ensuring high acceptance rates without tuning proposal distributions.
• Integration: The authors successfully integrate ESS with the Woodbury likelihood and Vecchia approximation. This allows for full posterior integration of latent variances rather than point estimation, providing a complete distribution of uncertainty.

C. Algorithmic Workflow

The inference proceeds via a Metropolis-within-Gibbs scheme:

1. Hyperparameters: Sampled via random-walk Metropolis.
2. Latent Variances ($\log \Lambda_n$): Sampled via ESS using the sparse Vecchia-approximated likelihood.
3. Prediction: Aggregates samples from the posterior to produce predictive means and variances, utilizing the law of total variance to account for uncertainty in both the mean and variance processes.

3. Key Contributions

1. Scalable Bayesian hetGP: The first implementation of a fully Bayesian hetGP that can handle millions of data points ($N$) and hundreds of thousands of unique inputs ($n$) by combining Woodbury identities and Vecchia approximations.
2. Full Uncertainty Quantification: Unlike previous scalable methods that rely on point estimates, bhetGP provides a full posterior distribution for latent variances, leading to more robust and calibrated uncertainty intervals.
3. Efficient ESS Implementation: Demonstrates that ESS is highly effective for high-dimensional latent GP fields, overcoming the mixing issues that plagued earlier MCMC approaches in this domain.
4. Open-Source Tool: The method is implemented in the bhetGP R package on CRAN, making these advanced capabilities accessible to practitioners.

4. Results and Evaluation

The authors evaluated bhetGP on two benchmarks:

A. Benchmark: Assemble-to-Order (ATO) Inventory Problem

• Setup: A stochastic simulation with $N=120,000$ observations.
• Comparison: Compared against standard hetGP (MLE-based) and a "full" dataset version.
• Findings: bhetGP achieved lower Root Mean Square Error (RMSE) and higher proper scores (better UQ) than the standard MLE-based hetGP. Crucially, bhetGP was the only method capable of utilizing the full dataset ($n=11,000$ unique inputs) effectively, whereas standard hetGP failed to scale.

B. Motivating Application: NOAA-GLM Lake Temperature

• Setup: A massive campaign with 9 million simulations ($N \approx 9,000,000$) across 300,000 unique input combinations (depths, days of year, forecast horizons).
• Comparison: Compared against Holthuijzen et al. (2024), a state-of-the-art SK-like method using Vecchia but relying on point estimates.
• Findings:
  • Accuracy: bhetGP consistently outperformed the SK-like method in RMSE.
  • Uncertainty: bhetGP provided narrower and more accurate 90% prediction intervals. The SK-like method often produced intervals that were either too wide (conservative) or failed to cover the true values (under-coverage) due to the lack of variance uncertainty integration.
  • Speed: bhetGP was significantly faster than full-N calculations and competitive with point-estimate methods, running in minutes rather than days.

5. Significance

This work represents a major advancement in surrogate modeling for stochastic computer experiments.

• Bridging the Gap: It bridges the gap between the theoretical rigor of fully Bayesian inference and the practical necessity of handling massive datasets.
• Ecological Impact: By enabling accurate, uncertainty-aware forecasting for lake temperatures at a continental scale, the method supports critical resource management and climate adaptation strategies.
• Future Directions: The framework opens the door for applying Bayesian optimization and active learning to massive stochastic simulation campaigns, tasks that were previously computationally prohibitive or statistically unreliable.

In summary, bhetGP solves the "large $N$, large $n$" problem in heteroskedastic modeling by replacing point estimation with efficient, scalable MCMC (ESS) and sparse matrix approximations, delivering superior predictive performance and uncertainty quantification.