Generative Bayesian Computation as a Scalable Alternative to Gaussian Process Surrogates

Nick Polson; Vadim Sokolov

cs.LGFeb 24, 2026

Generative Bayesian Computation as a Scalable Alternative to Gaussian Process Surrogates

Authors: Nick Polson, Vadim Sokolov

Organizations: Booth School of Business, University of Chicago. · Department of Systems Engineering and Operations Research, George Mason University.

Abstract

Gaussian process (GP) surrogates are the default tool for emulating expensive computer experiments, but cubic cost, stationarity assumptions, and Gaussian predictive distributions limit their reach. We propose Generative Bayesian Computation (GBC) via Implicit Quantile Networks (IQNs) as a surrogate framework that targets all three limitations. GBC learns the full conditional quantile function from input--output pairs; at test time, a single forward pass per quantile level produces draws from the predictive distribution. Across fourteen benchmarks we compare GBC to four GP-based methods. GBC improves CRPS by 11--26% on piecewise jump-process benchmarks, by 14% on a ten-dimensional Friedman function, and scales linearly to 90,000 training points where dense-covariance GPs are infeasible. A boundary-augmented variant matches or outperforms Modular Jump GPs on two-dimensional jump datasets (up to 46% CRPS improvement). In active learning, a randomized-prior IQN ensemble achieves nearly three times lower RMSE than deep GP active learning on Rocket LGBB. Overall, GBC records a favorable point estimate in 12 of 14 comparisons. GPs retain an edge on smooth surfaces where their smoothness prior provides effective regularization.

Figures & tables

Figure 1: GP vs. GBC on a 1D jump function (

n=100

, noiseless). The stationary GP (left) adapts its length scale near the discontinuity but blurs the jump; GBC (right) widens its quantile bands at the boundary. Details on the IQN architecture are in Section 3 .

Panel of Figure 2 — Figure 2: IQN architecture. The input $x$ and quantile level $\tau$ are processed through separate branches ( $f_{x}$ and $f_{\tau}\circ\phi$ ), merged via elementwise product, and decoded into a location estimate $\hat{\mu}$ and a quantile prediction $\hat{q}_{\tau}$ . The cosine embedding $\phi(\tau)$ enables smooth interpolation over quantile levels.

Method	RMSE	CRPS	90% Coverage
hetGP (Matern-5/2)	$\mathbf{23.81\pm 0.49}$	$12.56\pm 0.24$	$0.875$
GBC (IQN, single)	$23.92\pm 0.57$	$12.56\pm 0.29$	$0.803$
GBC (IQN, $K{=}5$ ens.)	$23.88\pm 0.57$	$\mathbf{12.49\pm 0.29}$	$0.830$

Table 1: Motorcycle benchmark ( MASS::mcycle ,

n=133

): mean

\pm

SE over 50 replicates. hetGP is purpose-built for heteroskedastic data; GBC-IQN, with default loss weights and no custom kernel or noise model, achieves marginally lower CRPS.

Figure 3: Motorcycle benchmark: stationary GP (top left), hetGP (top right), single IQN (bottom left), and

K{=}5

IQN ensemble (bottom right). Shaded region: 90% predictive interval. The stationary GP assigns near-constant uncertainty; hetGP adapts via a latent noise process; GBC adapts without a parametric noise model.

	GP (Matérn-5/2)		GBC (IQN)
$d$	RMSE	CRPS	RMSE	CRPS	90% Cov.
2	$2.514\pm 0.053$	$1.379\pm 0.024$	$\mathbf{2.168\pm 0.028}$	$\mathbf{1.224\pm 0.015}$	$0.894$
3	$2.762\pm 0.045$	$1.515\pm 0.022$	$\mathbf{2.197\pm 0.018}$	$\mathbf{1.236\pm 0.011}$	$0.885$
4	$3.122\pm 0.060$	$1.716\pm 0.030$	$\mathbf{2.245\pm 0.017}$	$\mathbf{1.267\pm 0.012}$	$0.845$

Table 2: BGP benchmarks (Park et al. 2026): mean

\pm

SE over 20 replicates. GBC-IQN achieves lower CRPS at all dimensions; the advantage grows with

d

as the GP’s global length-scale estimate becomes less effective.

Figure 4: BGP benchmark (

d=2

): predicted mean surface from stationary GP (left) vs. GBC/IQN (right). Dashed line: true partition boundary

a^{\top}x=0

. GP smooths across the jump; GBC captures the regime change.

Method	RMSE	CRPS	90% Coverage
GP (Matern 5/2)	$0.590\pm 0.016$	$0.378\pm 0.005$	$0.995$
GBC (IQN)	$\mathbf{0.520\pm 0.009}$	$\mathbf{0.323\pm 0.004}$	$0.991$

Table 3: Friedman

d=10

n=2{,}000

: mean

\pm

SE over 10 replicates. GBC-IQN achieves lower RMSE and CRPS. Both methods produce coverage well above the 90% nominal level, indicating conservative (overly wide) intervals.

	GP (Matern 5/2)		GBC (IQN)
$n$	Time (s)	RMSE	Time (s)	RMSE
500	2.9	1.26	9.6	0.91
1,000	9.7	0.76	10.4	0.71
2,000	47.6	0.64	17.7	0.53
5,000	infeasible		41.2	0.40
10,000	infeasible		70.6	0.36

Table 4: Friedman

d=10

scaling study: training time (seconds) and test RMSE at each sample size. GP is infeasible beyond

n=2{,}000

. Each row is a single run; the

n{=}2{,}000

RMSE values differ from Table 3 ’s 10-replicate means because of sampling variability in the training set.

Figure 5: Friedman

d=10

: training time (left) and test RMSE (right) vs. sample size. GP time grows as

\mathcal{O}(n^{3})

; GBC grows approximately linearly. GBC achieves lower RMSE at all

n

and remains feasible as

n

grows beyond GP’s reach.

	MJGP		GBC		GBC-Aug
Dataset ( $n_{\rm train}$ )	RMSE	CRPS	RMSE	CRPS	RMSE	CRPS
Phantom ( $n{=}9\text{k}$ , $d{=}2$ )	$0.053\pm 0.001$	$0.010\pm 0.000$	$0.095\pm 0.001$	$0.031\pm 0.001$ †	$\mathbf{0.044\pm 0.002}$	$\mathbf{0.009\pm 0.000}$ †
Star ( $n{=}9\text{k}$ , $d{=}2$ )	$0.085\pm 0.001$	$0.024\pm 0.001$	$0.124\pm 0.001$	$0.045\pm 0.001$ †	$\mathbf{0.059\pm 0.002}$	$\mathbf{0.013\pm 0.000}$ †
AMHV ( $n{=}8.5\text{k}$ , $d{=}4$ )	$0.027\pm 0.000$	$0.012\pm 0.000$	$\mathbf{0.024\pm 0.000}$	$\mathbf{0.011\pm 0.000}$	—
Michalewicz ( $n{=}90\text{k}$ , $d{=}4$ )	infeasible $(\geq 48$ h)		$\mathbf{0.064\pm 0.001}$	$\mathbf{0.031\pm 0.000}$ †	—

Table 5: Flowers et al. (2026) benchmark datasets: mean

\pm

SE over 20 replicates, full training sets. MJGP uses the authors’ code ( Flowers et al., 2026 ) . GBC-Aug augments the IQN input with a boundary-classifier probability (Section 4.2 ). Bold = best value across methods.

Dataset	$n$	Metric	GBC-AL	DGP+ALC	Ratio
2D Exponential	100	RMSE	$0.087\pm 0.003$	$0.072^{\ast}$	$0.83\times$
Rocket LGBB	300	RMSE	$\mathbf{0.0072\pm 0.0002}$	$0.0213\pm 0.001^{\ast}$	$\mathbf{2.95\times}$
Satellite GRACE	500	RMSPE	$1.19\pm 0.01$	$0.967^{\ast}$	$0.81\times$
SatMini ( $d{=}7$ )	60	RMSPE	$\mathbf{5.30\pm 0.15}$	$7.91\pm 0.28$	$\mathbf{1.49\times}$
7D Proxy	70	RMSE	$\mathbf{0.344\pm 0.012}$	$0.558\pm 0.007$	$\mathbf{1.62\times}$

Table 6: Active learning benchmarks: GBC-AL vs. DGP+ALC. Metric

\pm

SE over independent replicates. Ratio

>1

favors GBC (bold).

Figure 6: Rocket LGBB active learning: RMSE vs. training set size. Lines show the mean over 30 replicates (GBC-AL) and 50 replicates (DGP+ALC); shaded regions are 90% bands (5th to 95th percentile). DGP+ALC replicate count (50) reflects the published experimental design of Sauer et al. (2023) ; GBC-AL uses 30 replicates. The unequal counts do not alter the qualitative conclusion: the bands separate completely by

n\approx 150

, and every GBC replicate achieves lower RMSE than every DGP replicate at

n=300

Figure 7: Satellite GRACE active learning: RMSPE vs. training set size. GBC-AL leads for

n<300

; DGP+ALC overtakes after the crossover (dashed vertical line). For expensive simulators where the acquisition budget is limited, GBC provides better accuracy during the early phase.

Benchmark	GP Baseline	GP Metric	GBC Metric	$\Delta$	Note
Emulation (CRPS, lower is better)
BGP $d{=}2$ ( $n{=}1.6$ k)	Stationary GP	1.379	1.224	$\mathbf{+11\%}$	20/20 reps
BGP $d{=}3$ ( $n{=}1.6$ k)	Stationary GP	1.515	1.236	$\mathbf{+18\%}$	20/20 reps
BGP $d{=}4$ ( $n{=}1.6$ k)	Stationary GP	1.716	1.267	$\mathbf{+26\%}$	20/20 reps
Motorcycle ( $n{=}133$ )	hetGP	12.56	12.49	$\mathbf{+0.5\%}$ §	50 reps
Friedman $d{=}10$ ( $n{=}2$ k)	Stationary GP	0.378	0.323	$\mathbf{+14\%}$	10/10 reps

Table 7: Summary of GBC performance across all benchmarks. “GP baseline” denotes the method-of-choice GP variant for each problem type.

\Delta

metric is the percentage improvement of GBC over the GP baseline (positive = GBC better; CRPS for emulation, RMSE/RMSPE for AL);

k\times

denotes a ratio when the improvement exceeds 100%. Bold indicates GBC is better.

Surrogate Modeling

Surrogate modeling builds faster or more interpretable approximations of expensive simulations and black-box models.

See all 218 papers

Frugal Bayesian Optimization: Scalable Surrogates for Data- and Resource-Limited Discovery
Jul 31, 2026Panagiotis Krokidas, Christoforos Rekatsinas, Vassilis Sioros +3Surrogate Modeling Bayesian Optimization
Institute of Informatics and Telecommunications, National Centre for Scientific Research "Demokritos", Agia Paraskevi, Greece. · Department of Mechanical Engineering and Aeronautics, University of Patras, Patras Greece. · Department of Informatics and Telecommunications, National and Kapodistrian University of Athens, Athens, Greece. +4
Sequential sparse Gaussian process quantile regression
Jun 30, 2026Hugo Nicolas, Olivier Le MaîtreQuantile Regression Uncertainty Quantification
PLATON, CMAP · Inria, Center for Applied Mathematics, ´Ecole polytechnique, Institut Polytechnique de Paris, Palaiseau, France · CNRS, Center for Applied Mathematics, ´Ecole polytechnique, Institut Polytechnique de Paris, Palaiseau, France
Exploiting Gradients in Bayesian Inference of Expensive Simulators
Oct 8, 2026Šimon Soldát, Václav ŠmídlSimulation-Based Inference Bayesian Optimization
Department of Computer Science, Faculty of Electrical Engineering, Czech Technical University in Prague, Czech Republic · Institute of Plasma Physics of the Czech Academy of Sciences, Prague, Czech Republic

Jul 31, 2026J/K move · Enter open · S save

Frugal Bayesian Optimization: Scalable Surrogates for Data- and Resource-Limited Discovery

Panagiotis Krokidas, Christoforos Rekatsinas, Vassilis Sioros, Grigorios M. Chatziathanasiou+2

Institute of Informatics and Telecommunications, National Centre for Scientific Research "Demokritos", Agia Paraskevi, Greece. · Department of Mechanical Engineering and Aeronautics, University of Patras, Patras Greece. · Department of Informatics and Telecommunications, National and Kapodistrian University of Athens, Athens, Greece.+4

Surrogate Modeling Bayesian Optimization

Open paper AskPDF arXiv

Bayesian Optimization (BO) is widely adopted for data-efficient optimization in scientific and engineering applications, yet its computational cost is rarely evaluated alongside optimization performance. Here we present a systematic, compute-aware study of BO that evaluates surrogate models along two axes: optimization quality and computational frugality. Across eight benchmark functions and nine real-world datasets spanning materials science, mechanics, robotics, chemistry, and machine learning, we benchmark four surrogate models: Gaussian Processes, Random Forests, NGBoost, and Bayesian Adaptive Spline Surfaces. We show that Gaussian Process-based BO consistently incurs the highest time and memory overhead without delivering superior optimization or sample efficiency. In contrast, scalable alternatives achieve equal or better performance at a fraction of the computational cost. Motivated by these findings, we introduce a surrogate-recommendation framework that predicts the most suitable BO surrogate from inexpensive dataset characteristics. Together, these results establish FruBO as a reproducible, compute-aware baseline for Bayesian Optimization and provide practical guidance for surrogate selection under limited computational and experimental budgets.

Abstract

Surrogate Modeling

Explore similar work

Frugal Bayesian Optimization: Scalable Surrogates for Data- and Resource-Limited Discovery

Sequential sparse Gaussian process quantile regression

Exploiting Gradients in Bayesian Inference of Expensive Simulators