#  Contour Location for Reliability in Airfoil Simulation Experiments
using Deep Gaussian Processes

Annie S. Booth Corresponding author: Department of Statistics, NC State
University<EMAIL_ADDRESS>S. Ashwin Renganathan Department of Aerospace
Engineering, Penn State Robert B. Gramacy Department of Statistics, Virginia
Tech

###### Abstract

Bayesian deep Gaussian processes (DGPs) outperform ordinary GPs as surrogate
models of complex computer experiments when response surface dynamics are non-
stationary, which is especially prevalent in aerospace simulations. Yet DGP
surrogates have not been deployed for the canonical downstream task in that
setting: reliability analysis through contour location (CL). Level sets
separating passable vs. failable operating conditions are best learned through
strategic sequential design. There are two limitations to modern CL
methodology which hinder DGP integration in this setting. First, derivative-
based optimization underlying acquisition functions is thwarted by sampling-
based Bayesian (i.e., MCMC) inference, which is essential for DGP posterior
integration. Second, canonical acquisition criteria, such as entropy, are
famously myopic to the extent that optimization may even be undesirable. Here
we tackle both of these limitations at once, proposing a hybrid criteria that
explores along the Pareto front of entropy and (predictive) uncertainty,
requiring evaluation only at strategically located “triangulation” candidates.
We showcase DGP CL performance in several synthetic benchmark exercises and on
a real-world RAE-2822 transonic airfoil simulation.

Keywords: entropy, tricands, Pareto front, uncertainty quantification,
sequential design, active learning, surrogate, emulator

##  1 Introduction

Computer simulations are increasingly utilized for experimentation when direct
manipulation is impossible or infeasible. Although diverse applications abound
(e.g., Tietze,, 2015; Bremer,, 2019; Lippe et al.,, 2019), here we are
motivated by one in particular: physics-based simulations in aeronautics.
These involve computational fluid dynamics (CFD) simulation requiring the
numerical solution of a set of stiffly coupled nonlinear partial differential
equations satisfying conservation laws (e.g., Pamadi et al.,, 2004; Mehta et
al.,, 2014; Vassberg et al.,, 2008). A specific example is the aerodynamic
modeling of the interaction between an aircraft, aircraft wing, or wing
section (a.k.a.,“airfoils”) and the surrounding airflow at transonic speeds
(i.e., flight at $\sim 70-80\%$ of the speed of sound). Transonic flow
straddles the line between subsonic (flying slower than the speed of sound)
and supersonic (faster than sound) flow that makes it particularly complex to
model. When the computational costs of obtaining runs severely limits how many
may be collected, a “surrogate model” or “emulator”, which is a fitted
statistical model trained on data from a computer experiment, may be essential
for downstream tasks requiring model evaluation at un-tried inputs.

###  1.1 Surrogate modeling

The canonical surrogate is a Gaussian process (GP; Santner et al.,, 2018;
Gramacy,, 2020), but in their traditional form, GPs suffer from the assumption
of stationarity. They must impart similar dynamics across the entire input
space which limits their ability to cope with regime shifts, as are common in
aerospace simulations and CFD more widely. Aircraft experience “shocks” when
the sound barrier is broken. Several adaptations have been considered to allow
for non-stationary GP flexibility (see Sauer et al.,, 2023a, for a review),
but they all have drawbacks. Non-stationary kernels suit spatial applications
with low input dimension (Higdon et al.,, 1999; Paciorek and Schervish,,
2003). Partition- or treed-GPs divy up the input space (Gramacy and Lee,,
2008; Bitzer et al.,, 2023), sacrificing global scope. Locally approximate GPs
are data hungry (Gramacy and Apley,, 2015).

The deep Gaussian process (DGP; Damianou and Lawrence,, 2013) shares many of
these same ingredients, but has fewer drawbacks. DGPs use stationary GPs to
spatially warp the original inputs. Warped inputs were first proposed in the
spatial statistics community (Schmidt and O’Hagan,, 2003; Sampson and
Guttorp,, 1992), but have recently been popularized by machine learners who
showcased DGP prowess on large-scale classification and regression tasks
(Damianou and Lawrence,, 2013; Dunlop et al.,, 2018). The challenge in DGP
inference lies in learning high-dimensional latent warping variables. With
large training data sizes and a penchant for computational thrift, many in
machine learning embraced approximate variational inference (VI; Bui et al.,,
2016; Salimbeni and Deisenroth,, 2017).

Yet the supposed thriftiness of approximate VI often disappoints because the
algorithms require careful tuning. Moreover, optimization in lieu of full
posterior integration sacrifices uncertainty quantification (UQ), which is
essential to safety/reliability tasks, especially in data-poor settings. For
our aeronautic work we instead embrace the fully-Bayesian inferential scheme
of Sauer et al., (2023b), leveraging modern advances in Markov Chain Monte
Carlo (MCMC). This fully-Bayesian DGP has been shown to outperform VI-based
approximate DGP competitors (Sauer et al.,, 2022, 2023a) and to excel in
surrogate modeling settings where training data is limited (Sauer et al.,,
2023b).

###  1.2 An illustration

To make this concrete, foreshadow somewhat, and to help frame our novel
methodological contribution, consider a 7d model of aerodynamic flow around an
RAE-2822 transonic airfoil. The model, utilizing Reynolds Averaged Navier-
Stokes (RANS) equations, is solved via SU2 (Economon et al.,, 2016), a public
simulation suite. Here we fit surrogates to simulations obtained in a 500-run
space-filling design in the 7d input space, specifically: the Bayesian DGP
(DGP ESS) of Sauer et al., (2023b); a traditional stationary GP via plug-in
estimates (GP MLE); and the approximate VI DGP (DGP VI) of Salimbeni and
Deisenroth, (2017). Further details about these comparators and acronyms will
be revealed as the paper progresses. Predictive performance is assessed on a
hold-out 4,500-point test set via root mean squared error (RMSE, measuring
accuracy) and continuous ranked probability score (CRPS, measuring UQ;
Gneiting and Raftery,, 2007). See the left panels of Figure 1. Lower is better
for both (App. A). Observe that the Bayesian DGP is both more accurate and
offers better UQ.

Figure 1: Out-of-sample prediction metrics for surrogate fits to static LHS
designs of size $n=500$ from the RAE-2822 simulation. Failures defined at
$L/D<3$. Arrows along $y$-axes indicate which direction is preferred. Boxplots
show 10 MC repetitions.

Improved prediction is just the tip of an iceberg. Surrogate models are most
often a means to another end; they are used for downstream tasks such as
Bayesian optimization (BO; e.g., Jones et al.,, 1998), active learning (AL;
e.g., Cohn,, 1994), and calibration (e.g., Kennedy and O’Hagan,, 2001). Here
we are interested in contour location (CL), a common task in aeronautic
reliability analysis (e.g., Renganathan et al.,, 2022, 2023). In CL, the
surrogate’s task is to identify a specific level set in the response surface,
one that separates system “passes” from “failures”. In the case of our 7d
airfoil simulator, we define a failure contour that separates efficient v.
inefficient flight conditions and wing design. Another common example involves
separating safe v. unsafe conditions based on the amount of vibration, or
“flutter”, experienced by the aircraft (e.g., Stanford et al.,, 2022). A
surrogate must accurately locate failure regions in order to provide precise
quantification of failure probabilities.

The right panels of Figure 1 provide out-of-sample classification metrics
through sensitivity and specificity (App. A). Observe that the Bayesian DGP
also performs well on these metrics (higher is better for both). Further
discussion is reserved for Section 5. Yet there is further scope for
improvement; space-filling designs are, of course, inefficient for finding
lower-dimensional manifolds like level sets.

###  1.3 Sequential design for contour location

We are interested in finding inputs that separate “pass” and “fail” regimes
under a limited simulation budget. In this context, it makes sense to build up
the design sequentially to strategically hone in on this target, i.e., CL – a
form of AL. There are two keys to effective AL: (a) getting a good surrogate
fit to data you already have; and (b) using the surrogate to solve an
acquisition criteria to determine the next run of the simulator which will be
used to augment the data. These steps are repeated until the simulation budget
is exhausted or until some performance metric is satisfied. The (Bayesian) DGP
provides (a), but it is less amenable to (b) compared to the ordinary GP.

For CL, acquisition criteria target the contour through uncertainty in the
pass/fail prediction, which is quantifiable through surrogate posterior
predictions. The most popular criteria, based on entropy (Gray,, 2011),
converts classification probabilities into a metric which may be maximized for
50/50 probabilities, i.e., where classification uncertainty is highest. But
entropy and its variants have two downsides. One is that finding the input
location with highest entropy requires a multi-start numerical optimizer in
order to circumnavigate disparate regions (local optima) of high entropy,
separated by wide gulfs of numerically zero entropy. This difficulty is
compounded when there are disconnected level sets separating regimes. MCMC-
based inferential schemes make things more challenging still: the surrogate
(i.e., a DGP) doesn’t provide a single predictive surface but thousands,
precluding library-based optimization. Second, and perhaps more importantly,
the level set is a continuum, not a singleton as in BO and other AL variants.
Although entropy does indeed target the contour, it is famously myopic (e.g.,
Cole et al.,, 2022) as an acquisition function: it targets the most probable
boundary crossing point, not exploration of the entire manifold. The result is
clumped acquisitions.

These drawbacks have inspired a cottage industry of workarounds. One simple
solution involves evaluating entropy on a discrete set of space-filling
candidates, like Latin hypercube samples (LHS; McKay et al.,, 2000). This is
straightforward, and has some advantages over continuous numerics, including
parallel evaluation and simple modular implementation. (Evaluating entropy on
candidates is a post-processing step, whereas numerical optimization must be
embedded within the surrogate.) Yet after an initial handful of acquisitions,
most candidates reside in zero-entropy regions, resulting in essentially
random selection and a coarse view of the target contour. Alternatively, some
have adapted expected improvement ideas from BO to fit the CL setting,
attempting to balance both exploration and exploitation (Ranjan et al.,, 2008;
Bichon et al.,, 2008; Picheny et al.,, 2010). Others have found success within
stepwise uncertainty reduction (SUR) BO frameworks, in which entropy (and
other criteria) are numerically integrated over the input domain (Marques et
al.,, 2018; Chevalier et al.,, 2014). However, the requisite numerical
quadrature can be hard to do, and thus SUR-based strategies can be
computationally daunting in moderate input dimension. Cole et al., (2022)
suggested a hybrid-candidate-local numerical optimization to favor local
optima in the entropy surface, a sensible but ultimately ad hoc enterprise:
deliberately doing worse at one thing (solving for high entropy) to do better
at another (avoiding myopia). Nevertheless, it out-performs the methods cited
above in many exercises. But it requires derivatives, making it ill-suited to
MCMC.

Here we explore an idea from the BO literature that strategically allocates
candidates between existing design points through the use of Delaunay
triangulation (Gramacy et al.,, 2022). So-called “tricands” are space-filling
in an adaptive way, guaranteeing that there are candidates in promising
regions of the input space for acquisition. Optimization on surrogate
predictive quantities, cumbersome with MCMC-based surrogates, is replaced with
a geometric solver (Barber et al.,, 1996; Habel et al.,, 2019) that doesn’t
consult the surrogate at all. As candidates, tricands enjoy parallel
evaluation downstream with any surrogate, even MCMC-based ones. For BO, they
are superior to LHS candidates and numerically optimized acquisitions. For CL,
well-chosen candidates provide an opportunity to address multiple criteria
simultaneously, for example to correct the myopia in pure entropy-based
acquisition. Toward that end, we propose acquisitions on the “Pareto front” of
tricands evaluations of CL entropy and ordinary predictive uncertainty. This
is both effective, and enjoys a simple implementation as a post-processing
step on surrogate prediction. While similar to ideas of Bryan et al., (2005),
who utilize the product of posterior variance and entropy, our Pareto
acquisitions are agnostic to the relative scale of predictive variance.

###  1.4 Roadmap

Each ingredient in this scheme (DGP, tricands, Pareto front), separately
offers a performance bump over its state-of-the-art analogue (GP, derivative-
based acquisition, entropy-only). Over the course of this paper, we will offer
intuitive explanations of these improvements and demonstrate superiority
empirically with realistic benchmarking exercises. We acknowledge that these
ingredients, taken separately, may not comprise a substantial novel
contribution. But their application in the CL context is novel, both
separately and together, and they have never been used in an airfoil
reliability context. DGPs, whether VI- or MCMC-based, have never been
entertained for CL, and both tricands and Pareto acquisitions are new to CL.
To ease use downstream and support reproducibility, we provide a fully open
source implementation as a supplement to this document, and in public
repositories.111Code to reproduce all examples and figures may be found at
https://bitbucket.org/gramacylab/deepgp-ex/. DGP surrogate entropy
calculations are provided as an update to the deepgp package on CRAN (Booth,,
2023).

The remainder of this paper is laid out as follows. Section 2 provides a
review of the surrogate-CL state-of-the-art: GPs and entropy-based
acquisition. Section 3 details the main, modern ingredients – DGP surrogates,
triangulation candidates, and Pareto front acquisitions – and how they weave
together to form our contribution. Section 4 offers implementation details and
empirical results on synthetic examples. Section 5 returns to our motivating
airfoil simulation. We offer a brief discussion in Section 6.

##  2 Stationary GP Contour Location

Here we review the state-of-the-art in GP-based entropy contour location. The
more modern ingredients of our approach (DGPs, tricands, etc.) will be
reviewed as they are adapted to suit CL in Section 3.

###  2.1 Gaussian process surrogates

Let $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ denote a blackbox computer
simulator; $X_{n}$ be an $n\times d$ matrix of simulation inputs, with rows
$x_{i}^{\top}$; and let $y_{n}=f(X_{n})$ represent an $n$-vector of outputs. A
traditional GP prior assumes

$y_{n}\sim\mathcal{N}_{n}\left(\mu,\Sigma(X_{n})\right)\quad\textrm{where}\quad\Sigma(X_{n})^{ij}=\Sigma(x_{i},x_{j})=\tau^{2}\left(k\left(\sum_{h=1}^{d}\frac{(x_{ih}-x_{jh})^{2}}{\theta_{h}}\right)+\eta\mathbb{I}_{i=j}\right).$
(1)

While $\mu$ may be linear in columns of $X_{n}$, we specify $\mu=0$ without
loss of generality. Popular choices of kernel $k(\cdot)$ are the squared
exponential or Matèrn (Stein,, 1999), but any positive definite function will
suffice. Hyperparameters $\tau^{2}$, $\theta=[\theta_{1},\dots,\theta_{d}]$,
and $\eta$ govern the scale, lengthscale, and nugget/noise respectively. We
focus here on deterministic computer simulations, fixing $\eta=1\times
10^{-6}$ throughout, to interpolate observations. Vectorized $\theta$ allows
for coordinate-wise anisotropic modeling. Fixing
$\theta=\theta_{1}=\dots=\theta_{d}$ enforces isotropy. These unknown
hyperparameters may be estimated through the likelihood, either by maximizing
(MLE) or with MCMC. See Santner et al., (2018); Gramacy, (2020) for review.

Given data $\\{X_{n},y_{n}\\}$, posterior predictions at an $n_{p}\times d$
matrix of untried input locations $\mathcal{X}$ follow

$\displaystyle\mathcal{Y}\mid
X_{n},y_{n}\sim\mathcal{N}_{n_{p}}\left(\mu_{Y},\Sigma_{Y}\right)\quad\textrm{where}\quad\mu_{Y}$
$\displaystyle=\Sigma(\mathcal{X},X_{n})\Sigma(X_{n})^{-1}y_{n}$ (2)
$\displaystyle\Sigma_{Y}$
$\displaystyle=\Sigma(\mathcal{X})-\Sigma(\mathcal{X},X_{n})\Sigma(X_{n})^{-1}\Sigma(X_{n},\mathcal{X}).$

As a visual, consider the 2d “plateau” function displayed in the left panel of
Figure 2 and featured later in Section 4.2.

Figure 2: Left: plateau function in 2d. Middle and right: predictive
distribution along a slice in $x_{2}$ along with entropy calculations for an
ordinary (left) and deep (right) GP.

We trained a GP surrogate to data obtained from a random design of size
$n=15$; a 1d slice of the predictive surface (with fixed $x_{2}=0.2$) is
displayed in the center panel. For now, focus on just the blue lines/shading.
The GP offers an adequate nonlinear fit to the surface, considering the
sparsity of the training data in 2d, but it struggles to balance the disparity
between the flat regions and the steep drop. This GP is restricted by the
assumption of stationarity, a byproduct of $\Sigma(\cdot)$ being solely a
function of relative distances. A stationary GP will struggle to balance steep
drops with flat regions; it must compromise between oversmoothing the drop and
undersmoothing the flats. A DGP, previewed in the right panel, is more
flexible in handling this transition. Surrogate predicted mean and variance
estimates (blue lines) are essential components of AL enterprises like CL.

###  2.2 Entropy contour location

The most popular acquisition criterion for CL is entropy (e.g., Oakley,, 2004;
Marques et al.,, 2018; Cole et al.,, 2022). Let $g$ represent a threshold
level such that $f(x)>g$ indicates failure. Denote
$p_{x}=\mathbb{P}\left(f(x)>g\right)$ as the probability of observing a
failure at location $x$. The entropy criteria

$H(x)=-p_{x}\log(p_{x})-(1-p_{x})\log(1-p_{x}),$ (3)

will be high in regions of pass/fail uncertainty (i.e., $p_{x}\approx 0.5$).
Under a GP surrogate, failure probabilities may be evaluated through a
Gaussian CDF $\Phi$:

$p_{x}=\mathbb{P}\left(f(x)>g\right)=1-\Phi\left(\frac{g-\mu_{Y}(x)}{\sigma_{Y}(x)}\right),$
(4)

where $\mu_{Y}(x)$ and $\sigma_{Y}(x)=\sqrt{\Sigma_{Y}(x)}$ are provided in
Eq. (2). Consider again the center panel of Figure 2. Define $g=0$ such that
$y>0$ represents failure. The true location of the failure contour (where the
true line crosses $g=0$) is marked by the purple star and vertical line. The
red and black lines along the $x$-axis show entropy $H(x)$ evaluated along a
dense grid. Entropy is highest where the surrogate mean prediction crosses
$g=0$, but it is high in all areas where the 95% interval captures $g=0$ as
well. This is entropy’s desirable property – it hones in on the failure
contour. Yet as we acquire more points (i.e., beyond $n=15$ in 2d), entropy
will become peakier and will still favor acquisitions around the vertical
purple bar.

Notice, there is a local optima in the entropy surface at the right-most
boundary where $x_{1}=1$. This arises primarily because the GP predictive
equations indicate higher uncertainty at the edges of the input space. Many AL
criteria are coupled with theory indicating that eventually such disconnected
“high uncertainty” regions will recruit acquisitions because the criteria will
eventually be maximized there, after targeted sampling elsewhere saturates.
This is not true with entropy; it will never “in-fill”. A consequence of this
is that it will not explore the entire level set. Once it finds a highly
probable crossing, like the one along the slice in the figure, it will only
explore near its current best estimate. In the parlance of AL/BO, it is too
exploitative, too myopic (Renganathan et al.,, 2021). There are several
suggested solutions in the literature. For example, Marques et al., (2018)
integrate entropy over the input space to target the largest reduction in
global entropy rather than simply the largest value. This is a good idea in
principle, but unfortunately the integral is not available in closed form, and
the quadrature required is numerically fraught beyond 2d. This approach does
not compare favorably to simpler, more ad hoc methods. For example, Cole et
al., (2022) suggest deliberately coarsening the entropy search by only finding
local rather than global optima of the acquisition surface to encourage
exploration. This works well, and represents the state-of-the-art as far as we
are aware. However, it has the downside of requiring a tuning parameter
(trading off local versus global search) and derivative-based optimization
which is not amenable to MCMC-based surrogates.

##  3 Non-stationary DGP Contour Location

Here we describe our contribution: contour location (CL) with deep Gaussian
process (DGP) surrogates, triangulation candidates, and Pareto front
acquisitions.

###  3.1 Deep Gaussian process surrogates

DGPs are formed by layering stationary GPs. While this may be accomplished
with kernel convolutions (Dunlop et al.,, 2018), we prefer the lens of
functional composition. A two-layer DGP prior is defined as

$\begin{aligned} y_{n}\mid
W&\sim\mathcal{N}_{n}\left(\mu_{y},\Sigma(W)\right)\\\
w_{i}&\stackrel{{\scriptstyle\mathrm{ind}}}{{\sim}}\mathcal{N}_{n}\left(\mu_{w},\Sigma(X_{n})\right)\quad\textrm{for}\;\;i=1,\dots,p,\end{aligned}\quad\quad\textrm{where}\quad
W=\begin{bmatrix}w_{1}&w_{2}&\dots&w_{p}\end{bmatrix},$ (5)

and $\Sigma(\cdot)$ is provided in Eq. (1) with potentially unique
hyperparameters. Each $w_{i}$ is a conditionally independent stationary GP
over the inputs $X_{n}$. Together these “nodes” form the latent $W$ layer,
which serves as a warped version of the inputs and feeds into an outer
stationary GP connecting to the response $y_{n}$. The dimension of the latent
layer need not match the dimension of the input space, but $p=d$ has been
shown to work well (Sauer et al.,, 2023b). Prior means are adjustable, but
common choices are $\mu_{y}=0$ and either $\mu_{w}=0$ or $\mu_{w}=X_{n}$ (the
latter only being applicable when $p=d$). We have found settings where both
choices are advantageous, but here we opt for the simplest option
$\mu_{y}=\mu_{w}=0$.

Latent $W$ drives non-stationary flexibility. Its nodes may “stretch” inputs
in regions of high signal and “squish” inputs in regions of low complexity.
Yet the functional and multi-dimensional nature of this unknown quantity poses
an inferential challenge. One solution is to approximate the intractable DGP
posterior with a “close” one (in terms of Kullback-Leibler divergence) from a
known target family; this is known as approximate variational inference (VI;
Blei et al.,, 2017) and is popular in machine learning applications where
signal-to-noise ratios are low and data are abundant (Damianou and Lawrence,,
2013; Salimbeni and Deisenroth,, 2017). However, in our deterministic, data-
sparse surrogate modeling setting we opt for full posterior integration to
prioritize UQ. We accomplish this using the Bayesian DGP set-up of Sauer et
al., (2023b). Latent $W$ are sampled through elliptical slice sampling (ESS;
Murray et al.,, 2010). Kernel hyperparameters are also sampled through MCMC,
although these are more of a fine-tuning and their sampling may be replaced
with MLE-based alternatives (e.g., Ming et al.,, 2022). See Sauer, (2023) for
a thorough review of Bayesian DGP surrogates.

After burn-in, ESS provides posterior samples $w_{i}^{(t)}$ for $i=1,\dots,p$
and $t\in\mathcal{T}$ (we shall drop reliance on kernel hyperparameters here,
which may also be indexed by $t$). Samples in hand, prediction at new inputs
$\mathcal{X}$ requires first “warping” the predictive locations, $\mathcal{X}$
to $\mathcal{W}_{i}^{(t)}$, for each node ($i$) and each MCMC sample ($t$).
Each component layer is conditionally Gaussian (5), so this means applying Eq.
(2) with sampled $w_{i}^{(t)}$ in place of $y_{n}$. We follow Sauer et al.,
(2023b) in mapping through the posterior predictive mean,

$\mathcal{W}_{i}^{(t)}=\Sigma(\mathcal{X},X_{n})\Sigma(X_{n})^{-1}w_{i}^{(t)}$
(6)

instead of drawing from the full Gaussian posterior. Second, we column bind
$\mathcal{W}_{i}^{(t)}$ into the $n_{p}\times p$ matrix $\mathcal{W}^{(t)}$
and again apply Eq. (2) to obtain posterior moments

$\displaystyle\mu_{Y}^{(t)}$
$\displaystyle=\Sigma(\mathcal{W}^{(t)},W^{(t)})\Sigma(W^{(t)})^{-1}y_{n}$ (7)
$\displaystyle\Sigma_{Y}^{(t)}$
$\displaystyle=\Sigma(\mathcal{W}^{(t)})-\Sigma(\mathcal{W}^{(t)},W^{(t)})\Sigma(W^{(t)})^{-1}\Sigma(W^{(t)},\mathcal{W}^{(t)}).$

We thus obtain $|\mathcal{T}|$ samples for each posterior moment (typically in
the thousands). It is often beneficial to summarize these quantities using the
law of total variance,

$\mu_{Y}=\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}\mu_{Y}^{(t)}\quad\textrm{with}\quad\Sigma_{Y}=\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}\Sigma_{Y}^{(t)}+\mathbb{C}\mathrm{ov}(\mu_{Y}^{(t)})\quad\textrm{and}\quad\sigma_{Y}=\sqrt{\mathrm{diag}(\Sigma_{Y})}.$
(8)

To provide a concrete visual, we fit a two-layer Bayesian DGP to the data from
Figure 2. The right panel of that figure shows the DGP fit over the 1d slice.
Compared to the GP (middle panel), the DGP is able to more accurately capture
the plateau shape of the surface and the sloping drop between the two flat
regions. The full 2d DGP mean surface is shown in the left panel of Figure 3,
with training data locations indicated by open circles. The estimated failure
contour is indicated by the black dashed line; the true failure contour is
overlayed in blue for reference. The center panel displays a heat map of the
predicted standard deviation, $\sigma_{Y}$ from Eq. (8). Notice that
uncertainty is high in areas near the transition and far from training data.
The equivalent figure for the one-layer GP is relegated to Appendix B since
posterior uncertainty for a stationary GP is less interesting, being solely a
function of distance to training data.

Figure 3: Two-layer DGP predicted mean (left), standard deviation (center),
and post hoc entropy (right) for the plateau function (10): white/high,
red/low. Open circles indicate training data. Purple triangles and green
diamonds represent potential Pareto front acquisitions (Section 3.3) under a
gridded candidate set and “tricands” (Section 3.2), respectively.

Although the full DGP posterior is no longer strictly Gaussian, its components
are conditionally Gaussian, allowing us to utilize the failure probability
calculation of Eq. (4) on component pieces to evaluate entropy for CL. Using
Eq. (6–7) with $\mathcal{X}=x$ and taking expectation over $t$, we obtain

$H(x)=\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}\left[-p_{x}^{(t)}\log(p_{x}^{(t)})-(1-p_{x}^{(t)})\log(1-p_{x}^{(t)})\right]\quad\textrm{where}\quad
p_{x}^{(t)}=1-\Phi\left(\frac{g-\mu_{Y}^{(t)}(x)}{\sigma_{Y}^{(t)}(x)}\right).$
(9)

It may be more expedient to plug the summarized moments $\mu_{Y}$ and
$\Sigma_{Y}$ from Eq. (8) directly into the Gaussian CDF of Eq. (4). We refer
to this as post hoc entropy as it uses post-processed DGP moments. The right
panel of Figure 2 visualizes both entropy calculations along the 1d slice as
red and dashed black lines. They are equivalent for the one-layer GP since
there was no warping at play. The DGP post hoc entropy calculation has a
higher magnitude peak, but the relative shapes of the curves are very similar.
Both identify a single high-entropy region where the predicted surface crosses
the failure boundary. We have noticed this in all of our experiments (not
shown here), and so for simplicity and ease-of-use (the post hoc calculation
is more “plug-and-play”), we utilize post hoc entropy here on out. However our
software supports both options.

The right panel of Figure 3 displays a heat map of the post hoc entropy
surface for the two-layer DGP fit to the plateau function. Notice how it is
very peaky, with a ridge of maxima along the entire predicted contour. For
now, ignore the purple triangles and green diamonds, which are discussed later
in Section 3.3. Suppose we sought a new acquisition utilizing the entropy
criterion, based on this DGP surrogate. There are two clear obstacles. First,
finding optima in the entropy surface comes with a hefty computational price-
tag. The ridge highlighted in Figure 3 was found by evaluating the DGP
surrogate on a dense grid, a method that is infeasible in higher dimensions.
Numerical optimizers, on the other hand, require many successive evaluations
of singleton predictive locations. In a DGP, each predictive location must be
mapped through the latent layer (6), through to posterior moments (7), with
this entire process repeated $|\mathcal{T}|$-many times before final aggregate
predictive mean and variances are available (8). Doing this evaluation one-at-
a-time (without closed-form gradients), is incredibly inefficient. A common
workaround in AL (e.g., Gramacy and Lee,, 2009) and BO (e.g., Eriksson et
al.,, 2019) literatures involves deploying a limited, discrete set of
candidates. With candidates, the loop over $t\in\mathcal{T}$ need only be done
once. Although simple to implement, the fidelity of such candidate searches is
low in situations where the goal is to identify a small target (e.g., a
contour) in a big input space. Space-filling candidates, like LHS’s, are
likely to miss peaky areas of high entropy, unless the size of the candidate
set is prohibitively large (even in modest input dimension).

The second obstacle involves the entropy criterion itself; it is too myopic
for CL (see, e.g., Marques et al.,, 2018; Cole et al.,, 2022). Even if we are
able to find optima in the entropy surface, this criterion only minimally
accounts for posterior predictive uncertainty. In Figure 3, a portion of the
high entropy ridge is actually very close to existing training data locations.
Even though uncertainty is low here, entropy is still high. Acquisitions near
the outer edge of the predicted contour where uncertainty is higher, or in the
upper-right corner where there is no training data, would be better. Although
there is a local maxima in the entropy surface in the upper-right corner, this
local maxima will never overcome the global maxima near the predicted contour.
This is why the hybrid-local-optimization scheme of Cole et al., (2022), which
successfully identifies local optima in the entropy surface, outperforms
globally-optimized competitors. But further mitigation of entropy’s
limitations is warranted.

###  3.2 Triangulation candidates

Here we address the first issue, and to a lesser extent the second (which is
more squarely the subject of Section 3.3), by adapting a new candidate scheme
from the BO literature: triangulation candidates (Gramacy et al.,, 2022), or
“tricands”. Rather than allocating candidates at-random or in a space-filling
scheme, we allocate candidates strategically by “in-filling” the gaps in the
existing design $X_{n}$. Start with a Delaunay triangulation of $X_{n}$: a set
of line segments between points that divy up the space so no lines cross. The
left panel of Figure 4 provides a visual of this in 2d (for the same training
data we used in the 2d plateau function example). The other panels will be
narrated in Section 3.3. Open circles indicate observed training data
locations, and black solid lines form the Delaunay triangulation. In two
dimensions, this triangulation forms triangles; in higher dimension it forms
tetrahedra although we will still refer to them as triangles. Each triangle is
defined by a $(d+1)\times d$ matrix. Denote each triangle by $T_{j}$ for
$j=1,\dots,n_{T}$. The number of triangles ($n_{T}$) is affected by the size
of the training data ($n$) and the dimension of the input space ($d$).

Figure 4: Left: Delaunay triangulation with internal and fringe ($\alpha=0.9$)
“tricands”. Middle: A Pareto front set (blue dashed lines) for two
hypothetical variables. Right: Entropy vs. standard deviation for the DGP of
Figure 3 two ways: on a grid in $x$-space (gray dots), with the Pareto front
as purple triangles; via “tricands” (black open circles) with Pareto front as
green diamonds.

Tricands are allocated based on the triangulation in two realms: internal and
fringe. Internal candidates are located at the barycenter of each triangle:

$\tilde{x}_{j}=\frac{1}{d+1}\sum_{i=1}^{d+1}T_{j}[i,\;].$

Notice in the figure how they are spread out between existing design points.
For fringe candidates, denote $F_{j}$ for $j=1,\dots,n_{F}$ as the $d\times d$
matrix representing a facet (outer edge of the triangulation). Together, these
facets form the convex hull. Calculate the midpoint of each facet as
$\bar{F}_{j}=\frac{1}{d}\sum_{i=1}^{d}F_{j}[i,\;]$ for $j=1,\dots,n_{F}$. Then
take the normal vector $\vec{v}_{j}$ extending perpendicularly from $F_{j}$ at
the midpoint $\bar{F}_{j}$ (dashed arrows in Figure 4, left). Allocate a
fringe candidate along each $\vec{v}_{j}$, with the exact placement determined
by a tuning parameter $\alpha\in[0,1]$ representing the proportional distance
from the facet midpoint to the boundary of the space. For BO, Gramacy et al.,
(2022) found $\alpha=0.5$ to perform well. For CL, we choose $\alpha=0.9$ to
encourage acquisitions near the boundary – contours are often found near
boundaries, especially in reliability applications. Fringe candidates are
represented by the tips of dashed arrows in the figure. Crucially, tricands
(both internal and fringe) are agnostic to surrogate models; they rely only on
the geometry depicted by existing training locations $X_{n}$, spreading out to
encourage exploration while still offering a higher concentration of
candidates in more densely observed areas to allow for exploitation.

Denote the combined set of internal and fringe tricands as $\mathcal{X}_{N}$,
where $N=n_{T}+n_{F}$. As input dimension ($d$) and data size ($n$) increase,
the total number of tricands will grow. When $N$ is large, it is helpful to
limit the computational burden by setting a maximum size of the candidate set,
$N_{\textrm{max}}\ll N$. While a rudimentary approach would simply sub-sample
$N_{\textrm{max}}$ points from all $N$ tricands, we find it helpful to employ
a targeted sub-sampling approach to guarantee retention of some exploitave
candidates near the failure contour. This approach is inspired by the BO-based
sub-sampling of Gramacy et al., (2022), but it requires modification to suit
the CL setting. Rather than retaining candidates around a single point as in
BO, we desire candidates along the entire contour. To this end, we first order
the data observations based on their relative closeness to the known contour
limit $g$, i.e.,

     R> ordering <- order(abs(y - g))

Starting with the observation closest to $g$, we identify all triangles
$T_{j}$ and/or facets $F_{j}$ that were generated with this point, and
randomly select one of their candidates to retain. We repeat this process,
working down the ordering and retaining one candidate adjacent to each
location, until we have collected at least 10% of $N_{\textrm{max}}$. Then we
select remaining candidates, filling out $N_{\textrm{max}}$, completely at
random. By iteratively working down the ordering of closeness to the contour,
we are able to target the entire contour set. As with ordinary tricands, this
targeted sub-sampling does not rely on any surrogate information.

###  3.3 Pareto-front acquisitions

Tricands circumvent a continuous, derivative-based optimization (troublesome
with MCMC-based surrogates), by replacing it with a geometric calculation.
They also, to a degree, encourage exploration of the acquisition space by
placing candidates where predictive uncertainty is likely to be high. However,
we have found that this is not enough to correct the myopia of entropy. Toward
that end, we propose an acquisition criteria that strikes a deliberate balance
between entropy (exploitation) and predictive uncertainty (exploration).
Although approaches similar to the ones we describe have been suggested in
other AL contexts, combining other multitudes of acquisition criteria, we are
not aware of any such attempts for CL, or with tricands or similar candidate-
based schemes.

Consider entropy $H(x)$ and predictive standard deviation $\sigma_{Y}(x)$
evaluated over $x\in\mathcal{X}$ in the input space. When using either as a
lone acquisition criterion, higher is better. Now let
$c(x)=(H(x),\sigma_{Y}(x))\in\mathbb{R}^{2}$ represent a two dimensional
“criteria” and write $c(x)\prec c(x^{\prime})$, for
$x^{\prime}\in\mathcal{X}$, if $c(x^{\prime})$ strictly dominates $c(x)$,
meaning that it is better in both dimensions simultaneously:
$H(x)<H(x^{\prime})$ and $\sigma_{Y}(x)<\sigma_{Y}(x^{\prime})$. The Pareto
frontier, defined as
$P(\mathcal{X})=\\{x\in\mathcal{X}:\\{x^{\prime}\in\mathcal{X}:c(x)\prec
c(x^{\prime}),\,x\neq x^{\prime}\\}=\emptyset\\}$, is the subset of the input
space $\mathcal{X}$ separating dominated from non-dominated points. The
definition is similar for higher-dimensional criteria $c(x)\in\mathbb{R}^{m}$,
for $m\geq 3$; for more details see Goodarzi et al., (2014, Chapter 4).222Note
that most textbook definitions of Pareto fronts/non-dominated sets do not
involve an input, $x$-variable, working instead on generic $c$-criteria
values. However, we have introduced $x$ to make an explicit link to criteria
evaluated over an input space $\mathcal{X}$, as relevant to our CL enterprise.
To connect to classical definitions one could instead write
$c(P(\mathcal{X}))$.

The Pareto frontier is a compact set, and in surrogate modeling/AL contexts it
is possible to estimate and quantify uncertainty around it (e.g., Binois et
al.,, 2015; Luo et al.,, 2018) for arbitrary $x\in\mathcal{X}$. When
restricting to a discrete subset of the input space, like to tricands
$\mathcal{X}_{N}$, the Pareto front $P(\mathcal{X}_{N})$ is also discrete, and
is known as the non-dominated set. For example, consider random samples on
$c(x)\in\mathrm{R}^{2}$, displayed as Var 1 and Var 2 in the middle panel of
Figure 4. We generated these directly in $c-$space, without generating $x$’s
first, for a simple illustration. But we shall return to $x$-space
momentarily. The blue filled circles connected by the dashed staircase
indicate the points on the Pareto front, $P$ (assuming larger is preferred).
These are not beaten by any other point with respect to both variables: other
points may be higher on one or the other, but not both.

Identifying the Pareto front from a discrete sample in 2d is straightforward:
(a) order the observations from highest to lowest for one variable (points
that yield the maximum of either variable are always in the Pareto front set);
and (b) work down the ordering one-by-one, adding points to $P$ only if they
are higher on the second variable than all preceding points (which were higher
on the first variable). Our implementation is 5–6 lines of R code, available
in our Git repo. However we note that the rPref R-package on CRAN (Roocks,,
2016) offers faster computations for larger datasets in higher dimensions.

For a particular example in our CL setting, let us return to the plateau
surface of Figure 3, which involves two input dimensions, and now two
acquisition criteria $c(x)=(H(x),\sigma_{Y}(x))$. The right panel of Figure 4
combines the visuals from the center and right panels of Figure 3 by plotting
$H$ vs. $\sigma_{Y}$, evaluated over a dense grid of locations (gray dots),
and also over tricands $\mathcal{X}_{N}$ (open black circles). An ideal
acquisition would be in the upper right corner of this plot, with both high
entropy and high uncertainty. In both Figures 3 and 4, candidates along the
Pareto front are highlighted by purple triangles for the candidate grid, and
green diamonds for tricands. Figure 3 is providing the $x$-location of the
Pareto front(s), whereas Figure 4 shows the corresponding $c(x)$-value(s).

Observe in the Figure(s) that tricands are providing similar information at a
fraction of the cost compared to a dense grid. Any of the points on either
Pareto front (green diamonds or purple triangles) represents a sensible
acquisition according to the hybrid, two-dimensional criterion. Observe in
Figure 3 that all of these choices are far from the training data $X_{n}$,
thus correcting the mypoia of entropy. They are either very close to the
estimated contour, but away from $X_{n}$. Or, if they are not near the contour
they are very far from $X_{n}$. We prefer to select an acquisition completely
at random from $P(\mathcal{X}_{N})$. Going forward we use only tricands
$\mathcal{X}_{N}$ because grids that are dense enough for good resolution
become prohibitively large in higher input dimension.

It is worth noting that tricands have the extra advantage, over grids, that
the members of $P(\mathcal{X}_{N})$, the green diamonds, are located far apart
from one another and from the training design $X_{N}$ in the input, $x$-space,
by virtue of their geometric nature. This is not true for Pareto fronts
constructed from a dense grid. This means a tricands–Pareto set could be used
to select a whole batch of new runs at once (Cole et al.,, 2022). We provide
some further discussion on this idea in Section 6.

To summarize: our candidate-based contour location scheme boils down to three
steps: (1) build tricands $\mathcal{X}_{N}$ from $X_{n}$; (2) fit a DGP
surrogate to evaluate posterior predictive standard deviation
$\sigma_{Y}(\mathcal{X}_{N})$ and predictive entropy $H(\mathcal{X}_{N})$ for
these candidates; (3) select an acquisition uniformly at random on the Pareto
front of uncertainty and entropy:
$x_{n+1}\sim\mathrm{Unif}[P(\mathcal{X}_{N})$]. While this Pareto–tricands CL
scheme is technically independent of surrogate modeling choices and may be
deployed with any surrogate – and also any multi-dimensional acquisition
criteria – we find DGPs are crucial for the non-stationary response surfaces
that motivate this work. Empirical evidence is presented in Sections 4–5.

##  4 Implementation and Benchmarking

In this section we first describe the publicly available implementation of our
method, then validate it on a variety of synthetic test problems.

###  4.1 Implementation details

We utilize the deepgp package for R on CRAN (Booth,, 2023) for Bayesian DGP
surrogate training and prediction. That package has recently been updated to
provide entropy calculations following Eq. (9), i.e., using all posterior
samples $\mu_{Y}^{(t)}$ and $\sigma_{Y}^{(t)}$, via an entropy_limit argument
to the predict function.

     R> fit <- fit_two_layer(x, y, nmcmc = 10000, true_g = 1e-6)
     R> fit <- predict(fit, x_cand, entropy_limit = g)      ## g indicates failure threshold

All experiments in this paper utilize package defaults aside from fixing the
noise parameter $\eta=1\times 10^{-6}$ (true_g above) for interpolation of
deterministic computer models. Even though full entropy calculation is offered
in the package, we opt for the post hoc calculations using summarized
posterior moments in the simulation studies to follow. Specifically, using
posterior moments stored in the S3 class fit object, this calculation requires
only a short script:

     R> mu <- fit$mean                                         ## predicted mean for x_cand
     R> sigma <- sqrt(fit$s2)                                  ## predicted s.d. for x_cand
     R> p_x <- 1 - pnorm((g - mu) / sigma)                     ## failure probability
     R> ent <- -p_x * log(p_x) - (1 - p_x) * log(1 - p_x)      ## post hoc entropy

We use the geometry (Habel et al.,, 2019) R-package for Delaunay
triangulation, which is in-turn based on a legacy C library called quickhull
(Barber et al.,, 1996). We offer convenient R and python wrappers which
include calculation of internal and fringe “tricands” candidates as well as
the targeted sub-sampling approach described in Section 3.2 in our public git
repository.333https://bitbucket.org/gramacylab/deepgp-ex/ Identification of
the Pareto front requires a simple for loop over the ordered observations and
an if statement to check values of the second variable. Reproducible code,
including selection of the Pareto front, for all simulations (Sections 4.2-5)
is provided in the same repository. Our main performance metric, here and in
Section 4, is sensitivity (App. A) on a hold-out LHS testing set over
iterations of acquisition, or for a single space-filling batch of fixed size.
Higher sensitivity is better, with 1.0 indicating detection of all failures.

We entertain the following surrogates as competitors:

  * •

DGP ESS: two-layer Bayesian ESS-based DGP of Sauer et al., (2023b) using
deepgp in R

  * •

GP MCMC: stationary/one-layer GP with MCMC-sampled separable lengthscales,
also via deepgp

  * •

GP MLE: stationary GP with maximum likelihood estimated (MLE) separable
lengthscales using scikit-learn in Python (Pedregosa et al.,, 2011)

  * •

DGP VI: doubly stochastic VI DGP of Salimbeni and Deisenroth, (2017) using
gpytorch in Python (Gardner et al.,, 2018)

All utilize the Matèrn $\nu=5/2$ kernel, and all training data are observed
without noise. For each model, we compare a static space-filling LHS to a
strategic CL sequential design (when feasible). We deploy our tricands-Pareto-
front acquisition scheme with the MCMC-based surrogates (DGP ESS pareto and GP
MCMC pareto). Since these two differ only in the choice of surrogate (two-
layer DGP v. stationary GP), they offer a glimpse into the direct impact of
non-stationary flexibility on CL performance. We deploy the hybrid entropy-
based numerical optimization scheme of Cole et al., (2022) with the MLE-based
GP (GP MLE hyb ent), which is the state-of-the-art at the time of publication,
and thus our main competitor.

The DGP VI surrogate employs approximate VI through a stochastic optimization
of the evidence lower-bound (ELBO; minimizing the KL-divergence between the
intractible DGP posterior and the target distributional family). This
optimization entails several tuning parameters including the learning rate,
number of optimization steps (epochs), and the number of Gaussian mixtures in
the target posterior. There exist automated tuning methods involving inner-
optimizations (e.g., Zimmer et al.,, 2021); however, these approaches incur a
significant computational expense. In the context of sequential design for CL,
having an internal sequential design (i.e., a BO) to fine-tune optimization
parameters seems to us like a rabbit hole. Rather, we fix these values at what
we consider to be good defaults (mimicking the other surrogates which all
utilize package defaults): $1\times 10^{-3}$ learning rate, 2500 epochs, and
32 Gaussian mixtures. Even without fine-tuning, the stochastic optimization
required for VI does not mesh well into a sequential design. In an MCMC
surrogate, consecutive model updates may be initialized at previously burned-
in values to minimize the computational burden of subsequent MCMC; we know of
no effective equivalents for ELBO optimizations. Thus we only entertained the
DGP VI competitor in the static setting.

###  4.2 Synthetic examples

#### Plateau Function.

Consider the following function,

$f(x)=2\cdot\Phi\left(\sqrt{2}\left(-4-3\sum_{i=1}^{d}x_{i}\right)\right)-1\quad\textrm{for}\;x\in[-2,2]^{d}\;\left(\textrm{but
scaled to}\;x\in[0,1]^{d}\right),$ (10)

which we have adapted from Izzaturrahman et al., (2022) to be defined in
arbitrary dimension. This test function formed the basis of our earlier
illustrations in Figures 2–4. We define the contour at $g=0$ such that
$f(x)>0$ indicates a failure. Starting with $d=2$ and an initial LHS design of
size $n_{0}=5$, we conduct sequential designs targeting this contour up to an
ending design size of $n=30$ (25 acquisitions). Median sensitivity across 50
Monte Carlo (MC) repetitions with re-randomized starting designs is displayed
in the left panel of Figure 5. The distribution of sensitivities across MC
repetitions after the last acquisition are displayed in the left side of the
right panel. As an additional benchmark, we fit each surrogate on static LHS
designs of equivalent size ($n=30$) and reported the sensitivities in the same
panel.

Figure 5: Sensitivity (higher is better) for the 2d plateau function.

Although the DGP surrogate starts at a similar place as the MCMC-based GP, it
quickly jumps into the top position after a couple of acquisitions and
maintains this superiority for the extent of the design phase. All methods
benefited from CL acquisition when compared to their static LHS counterparts.
The DGP VI model struggled to cope with this small-data setting; it was
designed for much larger problems.

Figure 6: Sensitivity for the 5d plateau function.

To make things a little trickier, we bumped the dimension up to $d=5$ and
repeated this exercise. Higher dimension requires larger data sizes, so we
started with $n_{0}=20$ and extended the designs up to $n=150$. Results are
shown in Figure 6. Performance trends are similar to those of Figure 5,
although it appears the DGP has an even bigger edge in this higher-dimensional
setting. Also noteworthy, entropy-based acquisitions behind “GP MLE hyb ent”
occasionally led the GP astray and resulted in poorer performance than their
static design counterparts, a byproduct of clustered acquisitions. DGP VI LHS
performs better, comparatively, than in the 2d setting. We attribute this to
the larger training data size.

#### Cross-in-tray function.

Next, consider the “cross-in-tray” function,

$f(x_{1},x_{2})=-0.001\left(\left|\sin(x_{1})\sin(x_{2})\mathrm{exp}\left(\left|100-\frac{\sqrt{x_{1}^{2}+x_{2}^{2}}}{\pi}\right|\right)\right|+1\right)^{0.1}$

found on the pages of the Virtual Library of Simulation Experiments
(Surjanovic and Bingham,, 2013). We use the domain $x\in[-2,2]^{2}$ (similarly
scaled down to $[0,1]^{2}$). A visual of the surface is provided in the left
panel of Figure 7 with a contour defined at $g=2$.

Figure 7: Left: Heat map of 2d cross-in-tray function with contour at $g=2$.
Center: Sequential design from two-layer DGP with Pareto-tricands
acquisitions. Right: Sequential design from stationary GP with hybrid-entropy
acquisitions. Red “x” indicate starting LHS; black circles indicate acquired
points.

We follow the same sequential design procedures as described above, but we
increase training data sizes to $n_{0}=50$ and a final design size of $n=300$
to account for the increased complexity in the surface and the contour.
Results over 20 MC repetitions are displayed in Figure 8. Again, all methods
benefit from sequential design, but the DGP has a clear advantage due to its
non-stationary flexibility. The poor performance of DGP VI LHS here is
intriguing; perhaps it could be improved with fine-tuning. But we suspect that
the built-in inducing point approximations (Snelson and Ghahramani,, 2006),
which have been shown to hinder approximation fidelity (Sauer et al.,, 2022),
are restricting the model’s ability to handle the finer scale of this failure
region.

Figure 8: Sensitivity for the 2d cross-in-tray function.

In addition, the center and right panels of Figure 7 display the resulting
design at the end of a single MC exercise for the “DGP ESS pareto” scheme and
the “GP MLE hyb ent” scheme. Red “x” indicate the starting LHS design, and
black circles indicate the acquired points. As expected, acquisitions that
rely solely on entropy tend to cluster (right panel). Pareto front
acquisitions do a better job of exploring the contour. The real “secret sauce”
though is the combination of the Pareto-tricands acquisition scheme with the
right surrogate model (the flexible two-layer DGP in this case).

##  5 RAE-2822 Airfoil Computer Experiment

Here we investigate a computer simulation of transonic flow past the RAE-2822
airfoil, with physics governed by the RANS equations. We utilize the open-
source, finite-volume-based multiphysics code provided in the SU2 suite
(Economon et al.,, 2016). A spherical domain of $100$ airfoil chord length
radius is used to model the fluid domain, which is discretized with 27,857
triangular and hexahedral mesh elements. The hexahedral elements are employed
closer to the surface of the airfoil to capture the boundary layer.
Additionally, the mesh density is made deliberately higher in the near-field
of the airfoil to resolve the shock better; see Figure 9. Conservation laws
(mass, momentum, and energy) are enforced in each mesh element forming a
system of nonlinear equations which are solved iteratively. The computational
cost of the solution scales cubically with the number of elements, thus making
them expensive particularly when queried thousands of times, as in reliability
analysis. All simulations seek the steady-state solution, with pseudo time-
stepping. In serial execution, each simulation takes 20-30 minutes of
wallclock time for convergence, and about 150 seconds with 8 parallel MPI
processes.

Figure 9: Computational mesh. Left: Full domain showing all $27,857$ mesh
elements. Right: Near-field mesh refinement for the RAE-2822 airfoil to
capture the boundary layer and shock.

We consider seven input variables: three freestream conditions and four shape
parameters. Freestream conditions include angle of attack with bounds
$[0,10]$, Mach number ranging in $[0.7,0.9]$, and Reynolds number with bounds
$[5\times 10^{6},15\times 10^{6}]$. The airfoil shape is modified by
displacing free form deformation control points around the airfoil. Each
control point is an input variable, and we consider a total of four. The
displacement of the control points smoothly deforms the airfoil shape and the
surrounding mesh. The quantity of interest is the lift-drag ratio $L/D\equiv
C_{l}/C_{d}$, where $C_{l}$ and $C_{d}$ are the lift and drag coefficients,
respectively. $L/D$ is a measure of the aerodynamic efficiency of an aircraft
and directly impacts the amount of fuel required. Strict airworthiness
standards from the Federal Aviation Administration limit the amount of fuel an
aircraft can burn per passenger mile, since more jet fuel burned results in
more greenhouse gas emissions, and hence aircraft designers aim to meet a
minimum required $L/D$. Our current goal is to identify where $L/D<3$, which
is considered failure.

DGPs are relatively new to this area, but expert knowledge suggests there is
non-stationarity in the response surfaces of these simulations. To test this,
we fit DGP ESS, GP MLE, and DGP VI surrogates to a space-filling LHS of size
$n=500$. The compute time required to obtain this 500-point LHS was nearly 20
hours; it took over a week of runtime to generate the 4,500 point testing set.
Results were shown earlier in Figure 1, where we observed that the Bayesian
DGP indeed outperformed the stationary GP. The performance of DGP VI was
underwhelming; again we suspect fine-tuning may help, but the main culprit is
a blurry inducing point approximation.

Figure 10: Sensitivity (left) and specificity (right) for “DGP ESS pareto”
sequential designs for the 7d airfoil simulation (left panels). Right panels
show results from static LHS designs (DGP VI omitted).

To improve performance with fewer simulations, we deployed our Pareto-tricands
AL scheme for CL with the DGP ESS surrogate. We initialized with a random LHS
of size $n_{0}=100$ and acquired 400 points for a resulting design of size
$n=500$. Sensitivity and specificity (App. A, higher is better) for 10 MC
repetitions (with re-randomized initial and testing LHS sets) are shown in
Figure 10. The left interior panels show progress across the sequential
designs, and the right interior panels show the results from static designs of
equivalent size (copied from Figure 1 for reference). DGP-tricands-Pareto
sequential designs outperformed static fits with as few as 200 evaluations of
the simulator, a savings of nearly 12 hours of compute time compared to the
500-point LHS. Performance saturates after several hundred acquisitions, which
is a useful metric for determining a satisfactory stopping point for the
sequential design.

To help explain/visualize how the DGP ESS Pareto-tricands scheme is gaining an
advantage, Figure 11 shows one of the resulting sequential designs, visualized
as a projection over the two most impactful inputs: angle of attack and Mach
number.

Figure 11: 2d projection of a “DGP ESS pareto” sequential design ($n=500$).
Failures at $L/D<3$.

Pareto-tricands acquisitions (filled circles and triangles) occasionally
explore, but overwhelmingly focus on low angles of attack. All failures
occurred at low angles of attack and high speeds. This makes sense, as this is
where lift is generally low. Acquisitions were often placed near this boundary
in the upper left corner, an indication that the AL procedure honed in on the
contour.

##  6 Discussion

We proposed a sequential design scheme for contour location with deep Gaussian
processes that relies on triangulation candidates and selects acquisitions on
the Pareto front of entropy and posterior predictive uncertainty. Our proposed
scheme circumvents cumbersome numerical optimizations, avoids clustered
acquisitions, and allows for CL with Bayesian DGPs (which has not been done
before). While we were specifically motivated by the application of DGPs in
aeronautic simulations, our sequential design scheme is not limited to these
cases. Pareto-tricands sequential designs may be useful for other hefty
surrogate models which are incompatible with numerical optimization of
acquisition functions, such as an MCMC-based treed Gaussian process (Gramacy
and Lee,, 2008). We believe that the Pareto front acquisition procedure is
applicable to any combination of multiple criteria. One advantage over other
hybrid schemes, such as those based on aggregation, is that it is independent
of the relative scales of the component criteria. Combining any criteria with
predictive uncertainty and selecting acquisitions on the Pareto front has
potential to encourage exploration while retaining exploitative behavior,
which could be useful for a wide range of sequential design endeavors.

While tricands proved useful in our modest dimensional cases, the
computational costs of the Delaunay triangulation may become prohibitive in
higher dimensions. In our experience, 8 dimensions pushes the boundary of
tractability for the underlying quickhull library, unless the design size $n$
is kept low. There are some opportunities for improvement, such as only
updating the portion of the triangulation near a newly acquired point. But we
suspect that new methodology will be needed to replace the Delaunay
triangulation in higher dimension while still mimicking its behavior. One
advantage to tricands/Pareto acquisition, which was not explored in this
paper, involves batch acquisition. We observed in Section 3.3 that it would be
easy to take a batch whose size matched $|P(\mathcal{X})|$, the size of the
non-dominated set. We note that one could always get more, to augment the
batch, by removing those points and recalculating the Pareto front. Exploring
the extent to which this is a productive use of resources would be an
interesting subject of future research.

Reliability analysis involves quantifying the probability of failure given
uncertain inputs (e.g., Allen and Maute,, 2004). Our goal has been to train a
surrogate to accurately identify an entire failure contour. The surrogate can
then be utilized with any probability distribution over the input variables to
estimate probabilities of failure. Nevertheless, it may be advantageous to
incorporate information about a given input distribution into the sequential
design procedure, so the surrogate can prioritize learning around certain
regions of the contour. Abdelmalek-Lomenech et al., (2022) have shown such
methods to be advantageous with GPs, albeit at great computational expense. We
suspect similar strategies will work with DGPs, as long as the computational
burden can be kept in check.

## References

  * Abdelmalek-Lomenech et al., (2022) Abdelmalek-Lomenech, R. A., Bect, J., Chabridon, V., and Vazquez, E. (2022).  “Bayesian sequential design of computer experiments to estimate reliable sets.”  arXiv preprint arXiv:2211.01008.
  * Allen and Maute, (2004) Allen, M. and Maute, K. (2004).  “Reliability-based design optimization of aeroelastic structures.”  Structural and Multidisciplinary Optimization, 27, 228–242.
  * Barber et al., (1996) Barber, C. B., Dobkin, D. P., and Huhdanpaa, H. (1996).  “The Quickhull Algorithm for Convex Hulls.”  ACM Trans. Math. Softw., 22, 4, 469–483.
  * Bichon et al., (2008) Bichon, B. J., Eldred, M. S., Swiler, L. P., Mahadevan, S., and McFarland, J. M. (2008).  “Efficient global reliability analysis for nonlinear implicit performance functions.”  AIAA journal, 46, 10, 2459–2468.
  * Binois et al., (2015) Binois, M., Ginsbourger, D., and Roustant, O. (2015).  “Quantifying uncertainty on Pareto fronts with Gaussian process conditional simulations.”  European Journal of Operational Research, 243, 2, 386–394.
  * Bitzer et al., (2023) Bitzer, M., Meister, M., and Zimmer, C. (2023).  “Hierarchical-Hyperplane Kernels for Actively Learning Gaussian Process Models of Nonstationary Systems.”  In International Conference on Artificial Intelligence and Statistics, 7897–7912. PMLR.
  * Blei et al., (2017) Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2017).  “Variational inference: A review for statisticians.”  Journal of the American statistical Association, 112, 518, 859–877.
  * Booth, (2023) Booth, A. S. (2023).  deepgp: Deep Gaussian Processes using MCMC.  R package version 1.1.1.
  * Bremer, (2019) Bremer, S. A. (2019).  The GLOBUS model: Computer simulation of worldwide political and economic developments.  Routledge.
  * Bryan et al., (2005) Bryan, B., Nichol, R. C., Genovese, C. R., Schneider, J., Miller, C. J., and Wasserman, L. (2005).  “Active learning for identifying function threshold boundaries.”  Advances in neural information processing systems, 18.
  * Bui et al., (2016) Bui, T., Hernández-Lobato, D., Hernandez-Lobato, J., Li, Y., and Turner, R. (2016).  “Deep Gaussian processes for regression using approximate expectation propagation.”  In International conference on machine learning, 1472–1481. PMLR.
  * Chevalier et al., (2014) Chevalier, C., Bect, J., Ginsbourger, D., Vazquez, E., Picheny, V., and Richet, Y. (2014).  “Fast parallel kriging-based stepwise uncertainty reduction with application to the identification of an excursion set.”  Technometrics, 56, 4, 455–465.
  * Cohn, (1994) Cohn, D. (1994).  “Neural network exploration using optimal experiment design.”  In Advances in Neural Information Processing Systems, vol. 6, 679–686. Morgan-Kaurmann.
  * Cole et al., (2022) Cole, D. A., Gramacy, R. B., Warner, J. E., Bomarito, G. F., Leser, P. E., and Leser, W. P. (2022).  “Entropy-based adaptive design for contour finding and estimating reliability.”  Journal of Quality Technology, 1–18.
  * Damianou and Lawrence, (2013) Damianou, A. and Lawrence, N. D. (2013).  “Deep gaussian processes.”  In Artificial intelligence and statistics, 207–215. PMLR.
  * Dunlop et al., (2018) Dunlop, M. M., Girolami, M. A., Stuart, A. M., and Teckentrup, A. L. (2018).  “How deep are deep Gaussian processes?”  Journal of Machine Learning Research, 19, 54, 1–46.
  * Economon et al., (2016) Economon, T. D., Palacios, F., Copeland, S. R., Lukaczyk, T. W., and Alonso, J. J. (2016).  “SU2: An open-source suite for multiphysics simulation and design.”  Aiaa Journal, 54, 3, 828–846.
  * Eriksson et al., (2019) Eriksson, D., Pearce, M., Gardner, J., Turner, R. D., and Poloczek, M. (2019).  “Scalable Global Optimization via Local Bayesian Optimization.”  In Advances in Neural Information Processing Systems, vol. 32. Curran Associates, Inc.
  * Gardner et al., (2018) Gardner, J., Pleiss, G., Weinberger, K. Q., Bindel, D., and Wilson, A. G. (2018).  “Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration.”  Advances in neural information processing systems, 31.
  * Gneiting and Raftery, (2007) Gneiting, T. and Raftery, A. E. (2007).  “Strictly proper scoring rules, prediction, and estimation.”  Journal of the American statistical Association, 102, 477, 359–378.
  * Goodarzi et al., (2014) Goodarzi, E., Ziaei, M., and Hosseinipour, E. Z. (2014).  Introduction to optimization analysis in hydrosystem engineering.  Springer.
  * Gramacy and Lee, (2009) Gramacy, R. and Lee, H. (2009).  “Adaptive design and analysis of supercomputer experiments.”  Technometrics, 51, 2, 130–145.
  * Gramacy, (2020) Gramacy, R. B. (2020).  Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences.  Boca Raton, Florida: Chapman Hall/CRC.  http://bobby.gramacy.com/surrogates/.
  * Gramacy and Apley, (2015) Gramacy, R. B. and Apley, D. W. (2015).  “Local Gaussian process approximation for large computer experiments.”  Journal of Computational and Graphical Statistics, 24, 2, 561–578.
  * Gramacy and Lee, (2008) Gramacy, R. B. and Lee, H. K. H. (2008).  “Bayesian treed Gaussian process models with an application to computer modeling.”  Journal of the American Statistical Association, 103, 483, 1119–1130.
  * Gramacy et al., (2022) Gramacy, R. B., Sauer, A., and Wycoff, N. (2022).  “Triangulation candidates for Bayesian optimization.”  Advances in Neural Information Processing Systems, 35, 35933–35945.
  * Gray, (2011) Gray, R. M. (2011).  Entropy and information theory.  Springer Science & Business Media.
  * Habel et al., (2019) Habel, K., Grasman, R., Gramacy, R. B., Mozharovskyi, P., and Sterratt, D. C. (2019).  geometry: Mesh Generation and Surface Tessellation.  R package version 0.4.5.
  * Higdon et al., (1999) Higdon, D., Swall, J., and Kern, J. (1999).  “Non-stationary spatial modeling.”  Bayesian statistics, 6, 1, 761–768.
  * Izzaturrahman et al., (2022) Izzaturrahman, M. F., Palar, P. S., Zuhal, L., and Shimoyama, K. (2022).  “Modeling Non-Stationarity with Deep Gaussian Processes: Applications in Aerospace Engineering.”  In AIAA SCITECH 2022 Forum, 1096.
  * Jones et al., (1998) Jones, D. R., Schonlau, M., and Welch, W. J. (1998).  “Efficient global optimization of expensive black-box functions.”  Journal of Global optimization, 13, 4, 455–492.
  * Kennedy and O’Hagan, (2001) Kennedy, M. C. and O’Hagan, A. (2001).  “Bayesian calibration of computer models.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63, 3, 425–464.
  * Lippe et al., (2019) Lippe, M., Bithell, M., Gotts, N., Natalini, D., Barbrook-Johnson, P., Giupponi, C., Hallier, M., Hofstede, G. J., Le Page, C., Matthews, R. B., et al. (2019).  “Using agent-based modelling to simulate social-ecological systems across scales.”  GeoInformatica, 23, 269–298.
  * Luo et al., (2018) Luo, J., Gupta, A., Ong, Y.-S., and Wang, Z. (2018).  “Evolutionary optimization of expensive multiobjective problems with co-sub-Pareto front Gaussian process surrogates.”  IEEE transactions on cybernetics, 49, 5, 1708–1721.
  * Marques et al., (2018) Marques, A., Lam, R., and Willcox, K. (2018).  “Contour location via entropy reduction leveraging multiple information sources.”  Advances in neural information processing systems, 31.
  * McKay et al., (2000) McKay, M. D., Beckman, R. J., and Conover, W. J. (2000).  “A comparison of three methods for selecting values of input variables in the analysis of output from a computer code.”  Technometrics, 42, 1, 55–61.
  * Mehta et al., (2014) Mehta, P. M., Walker, A., Lawrence, E., Linares, R., Higdon, D., and Koller, J. (2014).  “Modeling satellite drag coefficients with response surfaces.”  Advances in Space Research, 54, 8, 1590–1607.
  * Ming et al., (2022) Ming, D., Williamson, D., and Guillas, S. (2022).  “Deep gaussian process emulation using stochastic imputation.”  Technometrics, 1–12.
  * Murray et al., (2010) Murray, I., Adams, R. P., and MacKay, D. J. C. (2010).  “Elliptical slice sampling.”  In The Proceedings of the 13th International Conference on Artificial Intelligence and Statistics, vol. 9 of JMLR: W&CP, 541–548. PMLR.
  * Oakley, (2004) Oakley, J. (2004).  “Estimating percentiles of uncertain computer code outputs.”  Journal of the Royal Statistical Society: Series C (Applied Statistics), 53, 1, 83–93.
  * Paciorek and Schervish, (2003) Paciorek, C. J. and Schervish, M. J. (2003).  “Nonstationary Covariance Functions for Gaussian Process Regression.”  In Proceedings of the 16th International Conference on Neural Information Processing Systems, NIPS’03, 273–280. Cambridge, MA, USA: MIT Press.
  * Pamadi et al., (2004) Pamadi, B., Covell, P., Tartabini, P., and Murphy, K. (2004).  “Aerodynamic characteristics and glide-back performance of langley glide-back booster.”  In 22nd Applied Aerodynamics Conference and Exhibit, 5382.
  * Pedregosa et al., (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al. (2011).  “Scikit-learn: Machine learning in Python.”  the Journal of machine Learning research, 12, 2825–2830.
  * Picheny et al., (2010) Picheny, V., Ginsbourger, D., Roustant, O., Haftka, R. T., and Kim, N.-H. (2010).  “Adaptive designs of experiments for accurate approximation of a target region.”  Journal of Mechanical Design, 132.
  * Ranjan et al., (2008) Ranjan, P., Bingham, D., and Michailidis, G. (2008).  “Sequential experiment design for contour estimation from complex computer codes.”  Technometrics, 50, 4, 527–541.
  * Renganathan et al., (2022) Renganathan, A., Rao, V., and Navon, I. (2022).  “Multifidelity Gaussian processes for failure boundary andprobability estimation.”  In AIAA SCITECH 2022 Forum, 0390.
  * Renganathan et al., (2021) Renganathan, S. A., Larson, J., and Wild, S. M. (2021).  “Lookahead Acquisition Functions for Finite-Horizon Time-Dependent Bayesian Optimization and Application to Quantum Optimal Control.”
  * Renganathan et al., (2023) Renganathan, S. A., Rao, V., and Navon, I. M. (2023).  “CAMERA: A method for cost-aware, adaptive, multifidelity, efficient reliability analysis.”  Journal of Computational Physics, 472, 111698.
  * Roocks, (2016) Roocks, P. (2016).  “Computing Pareto Frontiers and Database Preferences with the rPref Package.”  The R Journal, 8, 2, 393–404.
  * Salimbeni and Deisenroth, (2017) Salimbeni, H. and Deisenroth, M. (2017).  “Doubly stochastic variational inference for deep Gaussian processes.”  arXiv preprint arXiv:1705.08933.
  * Sampson and Guttorp, (1992) Sampson, P. D. and Guttorp, P. (1992).  “Nonparametric estimation of nonstationary spatial covariance structure.”  Journal of the American Statistical Association, 87, 417, 108–119.
  * Santner et al., (2018) Santner, T., Williams, B., and Notz, W. (2018).  The Design and Analysis of Computer Experiments, Second Edition.  New York, NY: Springer–Verlag.
  * Sauer et al., (2022) Sauer, A., Cooper, A., and Gramacy, R. B. (2022).  “Vecchia-approximated deep Gaussian processes for computer experiments.”  Journal of Computational and Graphical Statistics, 1–14.
  * Sauer et al., (2023a) — (2023a).  “Non-stationary Gaussian Process Surrogates.”  arXiv:2305.19242.
  * Sauer et al., (2023b) Sauer, A., Gramacy, R. B., and Higdon, D. (2023b).  “Active learning for deep Gaussian process surrogates.”  Technometrics, 65, 1, 4–18.
  * Sauer, (2023) Sauer, A. E. (2023).  “Deep Gaussian Process Surrogates for Computer Experiments.”  Ph.D. thesis, Virginia Tech.
  * Schmidt and O’Hagan, (2003) Schmidt, A. M. and O’Hagan, A. (2003).  “Bayesian inference for non-stationary spatial covariance structure via spatial deformations.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 65, 3, 743–758.
  * Snelson and Ghahramani, (2006) Snelson, E. and Ghahramani, Z. (2006).  “Sparse Gaussian processes using pseudo-inputs.”  Advances in Neural Information Processing Systems, 18, 1259–1266.
  * Stanford et al., (2022) Stanford, B., Sauer, A., Jacobson, K., and Warner, J. (2022).  “Gradient-Enhanced Reliability Analysis of Transonic Aeroelastic Flutter.”  In AIAA SCITECH 2022 Forum, 0632.
  * Stein, (1999) Stein, M. L. (1999).  Interpolation of spatial data: some theory for kriging.  Springer Science & Business Media.
  * Surjanovic and Bingham, (2013) Surjanovic, S. and Bingham, D. (2013).  “Virtual library of simulation experiments: test functions and datasets.”  http://www.sfu.ca/~ssurjano.
  * Tietze, (2015) Tietze, N. (2015).  “Model-based calibration of engine control units using gaussian process regression.”  Ph.D. thesis, Technische Universität.
  * Vassberg et al., (2008) Vassberg, J., Dehaan, M., Rivers, M., and Wahls, R. (2008).  “Development of a common research model for applied CFD validation studies.”  In 26th AIAA applied aerodynamics conference, 6919.
  * Zimmer et al., (2021) Zimmer, L., Lindauer, M., and Hutter, F. (2021).  “Auto-pytorch: Multi-fidelity metalearning for efficient and robust autodl.”  IEEE Transactions on Pattern Analysis and Machine Intelligence, 43, 9, 3079–3090.

APPENDICES

##  Appendix A Evaluation metrics

#### Pass/Fail Classification Performance.

For a given observation $f(x)$ with failures defined at $f(x)>g$ and surrogate
posterior predicted mean $\mu_{Y}(x)$, denote the classification of $x$ as

$C(x)=\begin{cases}\textrm{TP (true
positive)},&f(x)>g\;\;\textrm{and}\;\;\mu_{Y}(x)>g\\\ \textrm{FP (false
positive)},&f(x)\leq g\;\;\textrm{and}\;\;\mu_{Y}(x)>g\\\ \textrm{TN (true
negative)},&f(x)\leq g\;\;\textrm{and}\;\;\mu_{Y}(x)\leq g\\\ \textrm{FN
(false negative)},&f(x)>g\;\;\textrm{and}\;\;\mu_{Y}(x)\leq g\\\ \end{cases}$

Alternatively, swap all $>$ and $<$ signs if failures are instead defined at
$f(x)<g$. For an entire set of locations $\mathcal{X}$, let $|\textrm{TP}|$
represent the size of the set of rows of $\mathcal{X}$ classified as TP (and
so on). Then

$\mathrm{Sensitivity}=\frac{|\textrm{TP}|}{|\textrm{TP}|+|\textrm{FN}|}\quad\quad\textrm{and}\quad\quad\mathrm{Specificity}=\frac{|\textrm{TN}|}{|\textrm{TN}|+|\textrm{FP}|}.$

#### General Prediction Performance.

For $n_{p}$ testing locations with true observed values $y^{\mathrm{true}}$,
denote $\mu^{\star}$ as the predictive mean and $\sigma^{\star}$ as the
predictive standard deviation. Root mean squared error (RMSE) is defined as

$\mathrm{RMSE}=\sqrt{\frac{1}{n_{p}}\sum_{i=1}^{n_{p}}\left(\left(\mu_{i}^{\star}-y^{\mathrm{true}}_{i}\right)^{2}\right)}.$

Given a Gaussian posterior predictive distribution, continuous rank
probability score (CRPS; Gneiting and Raftery,, 2007) is defined as

$\mathrm{CRPS}\left(y^{\mathrm{true}}\mid\mu^{\star},\sigma^{\star}\right)=\frac{1}{n_{p}}\sum_{i=1}^{n_{p}}\left[\sigma^{\star}_{i}\left(2\phi(z_{i})+z_{i}\left(2*\Phi(z_{i})-1\right)-\frac{1}{\sqrt{\pi}}\right)\right]\quad\textrm{for}\quad
z_{i}=\frac{y_{i}^{\mathrm{true}}-\mu_{i}^{\star}}{\sigma_{i}^{\star}}$

where $\phi$ is the standard Gaussian pdf and $\Phi$ is the standard Gaussian
cdf.

##  Appendix B GP for plateau function

Figure 12 re-creates Figure 3 with a stationary GP surrogate instead of a two-
layer DGP.

Figure 12: One-layer GP predicted mean (left), standard deviation (center),
and corresponding entropy (right) for plateau function (10). Open circles
indicate training data.

The GP surrogate is parameterized by separable lengthscales sampled by
Metropolis-Hastings-within-Gibbs and was fit using the deepgp R-package
(Booth,, 2023). Compared to the non-stationary DGP, the stationary GP
struggles to balance the steep drop between regimes and offers a much softer
slope in the predictive surface that does not capture the plateaus as
accurately (left panel). The posterior variance of the GP (center panel) is
solely a function of distance to observed training data. Both the mean and
variance are involved in the calculation of the entropy surface (right panel).