Active Learning for Deep Gaussian Process Surrogates

1 Introduction

• Computer experiments can be costly, motivating carefully designed, limited simulation campaigns and fast surrogates with quantified uncertainty.Stationary GPs can struggle with regime changes, while existing non-stationary approaches often target specific cases or low-dimensional inputs.Source: §1 · pp. 1–2
• Hidden Gaussian layers warp input distances, enabling deep GPs to model non-stationarity, but their use as surrogates for simulation experiments remains underdeveloped.The setup combines warping, latent inputs, and composition within otherwise conventional kernel structures.Source: §1 · p. 2
• A hybrid Gibbs–Metropolis and elliptical slice sampling scheme enables DGP posterior sampling, retaining full uncertainty quantification while remaining tractable for typical surrogate-modeling scenarios.The authors favor this more expensive approach over approximate variational inference, which sacrifices full UQ.Source: §1 · p. 2
• The proposed DGP acquisition strategies use integrated mean-squared error and active learning Cohn, which the authors describe as better-performing alternatives to maximum-variance active learning.They defer DGP-based Bayesian optimization using expected improvement to future work.Source: §1 · p. 3
• The authors advocate a simple, limited DGP specification that outperforms ordinary GPs and deeper, more complex DGPs on benchmark and real-data surrogate modeling and design exercises.They provide the implementation in the deepgp R package on CRAN.Source: §1 · p. 3

2.1 Gaussian process surrogates

• Full posterior sampling lets GP surrogates average over hyperparameter uncertainty and avoid the overfitting seen with maximum-likelihood estimates, but stationary covariance can miss distinct regimes and GP computations scale cubically in training size.The review identifies dense n × n matrices as requiring quadratic storage and cubic decomposition.Source: §2.1 · p. 5

2.2 Active learning

• Sequential active learning adds one input at a time to reduce prediction error, while stationary GP variance criteria may favor space-filling locations instead of regions where the response is wiggly.The review describes active learning Cohn (ALC) as approximating integrated mean-squared prediction error using a reference set.Source: §2.2 · pp. 6–7

2.3 Deep Gaussian processes

• DGP layers warp inputs to enable region-dependent uncertainty and steer acquisitions toward interesting regions; up to three layers and no more latent nodes than input dimensions are sufficient for surrogate modeling and active learning.Two- and three-layer models compose intermediate GPs, redistributing inputs and allowing response values to be nonstationary.Source: §2.3 · p. 8
• Bayesian inference for two- and three-layer DGPs cannot analytically integrate out hidden layers; variational methods are cheaper but may oversimplify uncertainty, whereas MCMC represents uncertainty more thoroughly at greater computational cost.Efficient Markov-chain mixing is essential to mitigate MCMC costs.Source: §2.3 · p. 8

3 Modeling and inference

• The proposed inference scheme uses elliptical slice sampling (ESS) for posterior inference with DGPs.Source: §3 · p. 8

3.1 Model specification

• Noiseless hidden layers with unit scale, while placing scale τ 2 and nugget g only on the outer layer, supports stable posterior inference.Source: §3.1 · p. 8
• For latent W, the default is dim(W) = dim(Xn) = d, though smaller latent dimensions may suit high-dimensional inputs.Source: §3.1 · p. 8
• The model simplifies latent and output covariance by making W nodes conditionally independent, isotropic in inputs, and Y isotropic across W nodes.Source: §3.1 · p. 9
• These isotropy restrictions rein in complexity while allowing latent W to imitate separable or anisotropic kernels through its learned arrangement.Source: §3.1 · p. 9
• In the three-layer model, Z duplicates W’s conditional-independence and isotropy assumptions, has p ≤ d nodes, and connects each Z input to every W output.Source: §3.1 · p. 9, Figure 4 · p. 10

Hyperparameter priors

• The scale τ^2 uses the improper reference prior π(τ^2) ∝ 1/τ^2, whose posterior is proper for zero-mean models when n ≥ 1.Source: Hyperparameter priors · p. 10
• Gamma(3/2,b[·]) priors penalize extreme lengthscales and nugget values, with b[θy] > b[θw] > b[θz] favoring progressively less-wiggly deeper layers.For deterministic simulations, the outer-layer nugget is fixed at a small nonzero ε to interpolate the runs.Source: Hyperparameter priors · p. 10

3.2 MCMC posterior sampling

• The hybrid sampler uses random-walk Metropolis–Hastings for hyperparameters and ESS for latent layers; ESS improves mixing without tuning parameters, though it may require more likelihood evaluations.Source: §3.2 · pp. 10–12

3.3 Illustration

• ESS latent-layer samples bend inputs to separate flat interiors from wiggly extremes, while MH is initialization-sensitive and yields poorer out-of-sample RMSE under similar computational budgets.The MH chains also required step-size tuning and re-initialization, making burn-in difficult to assess.Source: Figure 5 · p. 13, §3.3 · p. 14

4.2 Acquisition

• The method averages each acquisition criterion across MCMC samples, then searches a candidate set; IMSE integrates over the sampled latent-space range, and ALC uses a reference set.The authors prefer the candidate-set search for pre-calculation, averaging across iterations, and parallel implementation.Source: §4 · p. 14, §4.2 · pp. 16–17

4.1 Prediction

• For prediction, latent-layer samples are propagated through the outer GP; substituting the latent posterior mean slightly underestimates uncertainty but has no detectable effect on acquisitions.The shortcut avoids sampling latent values at candidates, whose covariance calculation remains cubic in the number of candidates.Source: §4.1 · pp. 14–15

4.2 Acquisition

• In the Figure 2 example, both IMSE and ALC favor acquisitions in the left, wigglier region, while warping the ALC reference set produces slight criterion differences.The discrepancy can become more pronounced in three-layer models and higher.Source: §4.2 · p. 17

5.1 Implementation details

• The implementation uses layer-specific lengthscale priors to encourage near-identity latent mappings in inner layers while allowing more nonlinearity in outer layers.Hidden-layer dimension defaults to p = d; smaller p can help higher-dimensional problems uncover lower-dimensional latent structure, but may hurt low-dimensional problems unless most inputs are irrelevant.Source: §5.1 · p. 18
• After serial MCMC, predictions and acquisition criteria can be parallelized across iterations and candidate locations using provided S3 methods.When only pointwise means and variances are needed, predictions can instead be parallelized across test locations.Source: §5.1 · p. 18

5.2 Simulated data

• In the 1D simulation, the two-layer DGP outperformed the one-layer GP at small data sizes, and both deep models concentrated more runs in the wiggliest region.At larger n, the one-layer model caught up and edged the two-layer DGP on RMSE, while the two-layer model led on score; the three-layer model had wide quantiles from over-fitting.Source: §5.2 · pp. 19–20, Figure 6 · p. 19, Figure 7 · p. 20
• In the 2D simulation, the three-layer DGP led both RMSE and score on average but had the widest quantiles; both deep models focused acquisitions in the region of interest.The two-layer DGP achieved consistently higher scores than the one-layer GP after n > 35 acquisitions while having similar RMSE.Source: Figure 8 · p. 20, Figure 9 · p. 21, §5.2 · p. 21

5.3 Langley Glide-Back Booster

• Random 300-point designs substantially worsened RMSE for two- and three-layer DGPs, whereas their active-learning designs performed better.One-layer results were comparable between static and active designs; the passage attributes TGP performance to possible serendipitous accuracy rather than active-learning reinforcement.Source: §5.3 · p. 22, Figure 11 · p. 23
• Both deep models outperformed the one-layer GP but not TGP on LGBB RMSE, while allocating more than half their inputs to mach < 2.The one-layer design was similar to space filling, and TGP also favored mach < 2 but placed most points in the flatter region.Source: §5.3 · p. 22, Figure 11 · p. 23

5.4 Satellite drag

• DGPs with p = 3 achieved the 1% RMSPE benchmark with high probability in fewer than 500 runs, while one-layer GP and TGP did not.Increasing latent dimension improved results but cost more computation; p = 3 ran twice as fast as p = 7 on an eight-core machine.Source: Figure 12 · p. 24, §5.4 · p. 24

6 Conclusion

• DGP models outperformed typical GP regression in synthetic and real-simulator active-learning examples up to seven dimensions.The authors present these empirical outcomes as evidence across a relatively broad range of examples.Source: §6 · pp. 24–25
• Common latent-layer dimensionality, depth limited to three layers, and node counts no greater than input dimension regularized models, especially early in active learning.With small training designs, strong regularization encourages stability; later, likelihood information can dominate the prior.Source: §6 · pp. 24–25
• Two-layer DGPs with nodes matching input dimension compared well with more complex alternatives, but real-time selection of depth and node size remains of practical interest.The authors found no common simulation benchmarks demanding deeper models, while noting that deeper GPs can work with abrupt regime shifts and that two-layer models with fewer nodes than inputs were intriguing at higher dimensions. Input-connected networks may increase the utility of deeper models.Source: §6 · p. 25
• For calibration, multi-fidelity modeling, and black-box optimization, approximation via maximization can underestimate or mischaracterize uncertainty, while the proposed MCMC scheme may enable tractable inference and improve UQ.The cited applications are calibration (Marmin and Filippone, 2018), multi-fidelity modeling (Cutajar et al., 2019), and black-box optimization.Source: §6 · p. 25
• DGP fitting requires additional computation, and the proposed MCMC method is unsuitable for abundant data; inducing points, covariance formulations, or sequential model updates may help.Potential remedies include inducing points, covariance function or operator formulations, and sequential model updates.Source: §6 · p. 25

• The supplement illustrates the ALC surface for the 2d example.Source: SUPPLEMENTARY MATERIAL · p. 25
• The supplement provides a derivation of IMSE (18).Source: SUPPLEMENTARY MATERIAL · p. 25
• The supplement provides a derivation of ALC (19).Source: SUPPLEMENTARY MATERIAL · p. 25

A.1 2d ALC Surface

• For the stationary GP, the lower-left corner is highest-interest, but its ALC is essentially a function of proximity to Xn.The surface was computed using a dense Xref grid after an LHS of size n = 30.Source: §A.1 · p. 31
• The two-layer DGP allocates uncertainty to the interesting region and produces non-space-filling acquisitions, unlike the typical GP’s proximity-driven ALC.Both surrogates use the same training data; the ALC surfaces are shown for f(x1, x2) in Eq. (21).Source: Figure 13 · p. 31, §A.1 · p. 31

B.1 Partitioned Matrix Inverse

• The partitioned-matrix derivation uses k_n(w_{n+1}) = K_{θy}(W_n, w_{n+1}), following matrix properties such as Barnett (1979).The supplied passage gives the relation as part of the derivation.Source: §B.1 · p. 31
• It defines v = 1 + g − k_n(w_{n+1})^T C_n^-1 k_n(w_{n+1}) and h = −v^-1 C_n^-1 k_n(w_{n+1}).These expressions specify v and h in terms of the kernel vector and C_n^-1.Source: §B.1 · p. 31

B.2 IMSE Derivation

• The closed-form IMSE derivation of Binois et al. (2019) is extended to uniform measure over domains with coordinate bounds [a_i, b_i].The domain is W_i ∈ [a_i, b_i] for i = 1, . . . , p.Source: §B.2 · p. 31
• The model uses an isotropic Gaussian kernel and sets Σ(w) = 1 rather than 1 + g to target the variance of the mean.The stated model is Y ∼ N_n(0, K_θy(W) + gI).Source: §B.2 · p. 32
• Using W_n+1 and C_n+1, the IMSE expectation reduces to a trace of matrix products.The reduction follows Lemma 3.1 of Binois et al. (2019).Source: §B.2 · p. 32
• The notation defines w_j,i as the jth element of the ith node of W and Φ as the standard Gaussian CDF.The section also defines the elements of H for the derivation.Source: §B.2 · p. 32

B.3 ALC Derivation

• ALC is re-expressed from Eq. (19) using the new latent element w_{n+1}, with the (t) indexing omitted.The quantities k_n(w), v, and h are defined in Appendix B.1.Source: §B.3 · p. 32