Generative Bayesian Computation as a Scalable Alternative to Gaussian Process Surrogates

1. Problem Statement

Gaussian Process (GP) surrogates are the standard tool for emulating expensive computer experiments (e.g., calibration, sensitivity analysis, Bayesian optimization). However, they face three fundamental structural limitations that restrict their applicability in modern, large-scale, or complex settings:

1. Computational Scalability: Standard GPs require $O(n^3)$ operations for training (Cholesky decomposition of the covariance matrix) and $O(n^2)$ for prediction, making them infeasible for training sets larger than $n \approx 10,000$.
2. Stationarity Assumption: GPs assume a stationary covariance structure (correlation depends only on distance). They struggle with non-stationary phenomena such as regime changes, heteroskedasticity (varying noise levels), and jump discontinuities without complex, problem-specific extensions.
3. Gaussian Predictive Distributions: GPs inherently produce Gaussian predictive distributions. This limits their ability to model non-Gaussian posteriors, provide accurate tail probabilities, or estimate conditional quantiles directly, which are crucial for risk assessment and robust decision-making.

Existing solutions (e.g., Deep GPs, Vecchia approximations, Modular Jump GPs) address these issues individually but often introduce new complexities, require specific inference pipelines, or fail to scale simultaneously across all dimensions.

2. Methodology: Generative Bayesian Computation (GBC)

The authors propose Generative Bayesian Computation (GBC) using Implicit Quantile Networks (IQNs) as a unified surrogate framework that addresses all three GP limitations simultaneously.

Theoretical Foundation: Noise Outsourcing

The method relies on the Noise Outsourcing Theorem (Kallenberg, 1997), which states that any conditional distribution $\pi(\theta | y)$ can be represented as a deterministic function of the conditioning variable and a uniform random variable:
$$\theta = G^*(\tau, y), \quad \tau \sim U[0, 1]$$
In the surrogate context, the goal is to learn the generator $G_\phi(\tau, x)$ that maps inputs $x$ and a noise variable $\tau$ to the response $y$. This effectively learns the full conditional quantile function $Q_\tau(Y|x)$.

Architecture: Implicit Quantile Network (IQN)

• Input Embedding: The quantile level $\tau$ is embedded via a cosine transform $\phi(\tau) = [\cos(j\pi\tau)]_{j=0}^{n_h-1}$, allowing the network to represent the entire quantile function continuously.
• Structure: The network processes the input $x$ and the embedded $\tau$ through separate branches, merges them via element-wise multiplication, and outputs:
  1. A location estimate $\hat{\mu}$ (used as an L1 anchor during training).
  2. A quantile estimate $\hat{q}_\tau$ (the predictive output).
• Advantage: Unlike standard quantile regression which trains separate models for discrete $\tau$, a single IQN represents the entire conditional distribution.

Training Objective

The IQN is trained using a three-term loss function to prevent mode collapse, ensure quantile ordering, and minimize quantile error:
$$\mathcal{L} = w_1 \mathbb{E}[|y - \hat{\mu}|] + w_2 \mathbb{E}[m_\tau] + w_3 \mathbb{E}[\rho_\tau(y - \hat{q}_\tau)]$$

• $L_1$ Anchor: Encourages the network to track the conditional median.
• Ordering Surrogate: A soft penalty to prevent quantile crossing (e.g., ensuring $\hat{q}_{0.1} < \hat{q}_{0.9}$).
• Quantile Loss: The standard check loss ($\rho_\tau$) for accurate quantile estimation.
• Adaptability: Weights $(w_1, w_2, w_3)$ can be tuned. For jump discontinuities, a "quantile-dominant" configuration (reducing $w_1$) improves performance by allowing sharper transitions.

Extensions

• GBC-Aug (Boundary Augmentation): For sharp jump processes, the authors propose a three-step pipeline: (1) EM clustering to identify regimes, (2) training a classifier to predict regime probabilities, and (3) appending this probability as an input feature to the IQN. This allows the network to localize quantile estimates without learning the boundary geometry from scratch.
• Active Learning: Uses a Randomized-Prior Ensemble (K=3) where each member has a frozen random-weight network added to its output. This maintains diversity in predictions (crucial for acquisition functions) where standard bootstrap ensembles fail due to convergence to the same empirical risk minimizer.

3. Key Contributions

1. Unified Framework: Formulates GBC with IQNs as a single architecture addressing cubic scaling, non-stationarity, and non-Gaussianity.
2. Theoretical Guarantees: Proves that IQNs are universal approximators for conditional quantile functions and establishes convergence rates for the Continuous Ranked Probability Score (CRPS).
3. Boundary Augmentation: Introduces GBC-Aug, which matches or outperforms specialized Modular Jump GPs (MJGP) on jump datasets while training 26x faster.
4. Robust Active Learning: Proposes a randomized-prior ensemble strategy that provides reliable uncertainty signals for sequential design, outperforming Deep GP active learning on non-stationary benchmarks.
5. Comprehensive Empirical Evaluation: Systematically compares GBC against four GP-based methods (Stationary GP, hetGP, MJGP, DGP+ALC) across 14 benchmarks (dimensions $d=1$ to $10$, sample sizes $n=133$ to $90,000$).

4. Key Results

The paper evaluates GBC against GP baselines across various regimes:

• Scalability: On the Friedman function ($d=10$), GBC scales linearly to $n=90,000$ with $O(n)$ training cost, whereas dense-covariance GPs become infeasible beyond $n=2,000$. GBC achieves lower RMSE and CRPS at all tested sample sizes.
• Non-Stationarity & Jumps:
  • Piecewise Benchmarks (BGP): GBC improves CRPS by 11–26% over stationary GPs as dimensionality increases.
  • Jump Datasets (Phantom/Star): GBC-Aug reduces CRPS by 46% on the Star dataset compared to MJGP and matches MJGP on Phantom, with significantly faster training times.
  • Heteroskedasticity: On the motorcycle dataset, GBC (without a custom noise model) achieves CRPS comparable to the purpose-built hetGP.
• Active Learning:
  • Rocket LGBB (Non-stationary): GBC-AL achieves 2.95x lower RMSE than Deep GP (DGP+ALC).
  • Satellite GRACE (Smooth): DGP outperforms GBC at large sample sizes ($n=500$) due to the GP's smoothness prior, but GBC leads in the early phase ($n < 300$) where data is scarce.
• Overall Performance: GBC records a favorable point estimate (lower error) in 12 of 14 direct comparisons.

5. Significance and Conclusion

This paper establishes Generative Bayesian Computation as a scalable, flexible alternative to Gaussian Processes for computer experiments.

• When to use GBC: It is the superior choice for problems with discontinuities, high dimensionality, large training sets ($n > 5,000$), or non-Gaussian predictive requirements. It eliminates the need for custom kernel engineering or complex multi-stage pipelines.
• When GPs remain relevant: GPs retain an advantage on smooth, low-to-moderate dimensional surfaces with small sample sizes ($n < 5,000$), where their inductive bias (smoothness prior) provides effective regularization that neural networks cannot learn from limited data.
• Impact: GBC offers a "plug-and-play" surrogate that learns the full conditional distribution in $O(n)$ time, making it suitable for modern large-scale simulations where traditional GPs are computationally prohibitive. The authors suggest a hybrid future where practitioners select the tool based on problem structure or combine them (e.g., using GBC for emulation and conformal prediction for calibration).