Evaluating Gaussian process metamodels and sequential designs for noisy level set estimation

We consider the problem of learning the level set for which a noisy black-box function exceeds a given threshold. To efficiently reconstruct the level set, we investigate Gaussian process (GP) metamodels. Our focus is on strongly stochastic simulators, in particular with heavy-tailed simulation noise and low signal-to-noise ratio. To guard against noise misspecification, we assess the performance of three variants: (i) GPs with Student-t observations; (ii) Student-t processes (TPs); and (iii) classification GPs modeling the sign of the response. In conjunction with these metamodels, we analyze several acquisition functions for guiding the sequential experimental designs, extending existing stepwise uncertainty reduction criteria to the stochastic contour-finding context. This also motivates our development of (approximate) updating formulas to efficiently compute such acquisition functions. Our schemes are benchmarked by using a variety of synthetic experiments in 1–6 dimensions. We also consider an application of level set estimation for determining the optimal exercise policy of Bermudan options in finance.


Statement of problem
Metamodeling has become widespread for approximating black-box functions that arise in applications ranging from engineering to environmental science and finance (Santner et al. 2013).Rather than aiming to capture the precise shape of the function over the entire region, in this article we are interested in estimating the level set where the function exceeds some particular threshold.There is also research on this problem under the name "regression level sets" (Scott and Davenport 2007;Willett and Nowak 2007;Yang et al. 2014).Level set estimation is common in contexts where we need to quantify the reliability of a system or its performance relative to a benchmark.It also arises intrinsically in control frameworks where one wishes to rank the payoff from several available actions (Hu and Ludkovski 2017).
We consider a setup where the latent f : D → R is a continuous function over a d-dimensional input space D ⊆ R d .The level set estimation problem consists in classifying every input x ∈ D = S ∪ C according to (1.1) Without loss of generality, the threshold is taken to be zero, so that the level set estimation is equivalent to learning the sign of the response function f .For later use, we also define the corresponding zero contour of f , namely the partition boundary ∂ S = ∂C = {x ∈ D : f (x) = 0}.(Note that ∂ S and ∂C do not mean the topological boundary of the sets.)For any x i ∈ D, we have access to a simulator Y i that generates noisy samples of f (x i ): where i are realizations of independent mean zero random variables with variance τ 2 (x), the distribution of which depends on x i only.The simulator represents realizations of a stochastic system, capturing settings where the output is a function not just of the inputs, but also of some intrinsic random shocks.Those shocks might represent random awards available to a decision maker, or random disturbances that might tip the system into a failed state.The shocks typically have a simple structure (i.i.d uniform or Gaussian) but are fed through a nonlinear black-box transformation to produce observable system outputs with unknown statistical properties.Such stochastic simulators, where different runs with same x yield distinct Y i 's, are ubiquitous across scientific, engineering and financial domains that rely on stochastic models (Baker et al. 2020).
To assess a level set estimation algorithm, we compare the resulting estimate S with the true S in terms of their symmetric difference.Let μ be a probability measure on the Borel σ -algebra B(D) (e.g., μ = Leb D the Lebesgue measure on D).Then, our loss function is L(S, S) = μ(SΔ S), (1.3) where S 1 ΔS 2 := (S 1 ∩ S c 2 ) (S c 1 ∩ S 2 ), and S c is the complement of S: S c = {x ∈ D : x / ∈ S}.Frequently, the inference is carried out by first producing an estimate f of the response function; in that case, we take S = {x ∈ D : f (x) ≥ 0} and rewrite the loss as (1.4) where I(•) is the indicator function.

Motivation
As a concrete example of level set estimation, consider the problem of evaluating the probability of failure, determined via the positive level set S of a performance function f (•) (Vazquez and Bect 2009;Picheny et al. 2010;Bect et al. 2017).The system is safe when f (x) ≤ 0 and fails otherwise.In the context where the performance function can be evaluated via deterministic experiments, the estimation of the safe zone (more precisely its volume μ(S)) was carried out in Bect et al. (2012) and Mukhopadhyay et al. (2005) employing a Gaussian process metamodel with a sequential design.
A related case study dealing with the probability of failure in a nuclear fissile chain reaction appeared in Chevalier et al. (2014).
Another application, which motivated this present investigation, comes from simulation-based algorithms for valuation of Bermudan options (Gramacy and Ludkovski 2015;Ludkovski 2018).In that context, described in Sect.5, there is a sequence of problems indexed by t.For each t, S t is the continuation region which is characterized as the zero level set of the timing value f (x; t).The stochastic interpretation links f (x; t) to a conditional expectation of a path-dependent functional of a Markov process X , and the loss function (1.3) corresponds to the quality of the estimated stopping rule relative to the underlying distribution μ(•; t) of X t .

Design of experiments for contour finding
Reconstructing S via a metamodel can be divided into two steps: the construction of the response model and the development of methods for efficiently selecting the simulation inputs x 1:n , with n as the design size, known as design of experiments (DoE).Since the level set is intrinsically defined in terms of the unknown f , an adaptive DoE approach is needed to select x n 's sequentially.
For the response modeling aspect, GP regression, or kriging, has emerged as the most popular nonparametric approach for both deterministic and stochastic black-box functions (Bect et al. 2012;Gramacy and Lee 2009;Picheny et al. 2013a;Jalali et al. 2017).GPs have also been widely used for the level set estimation problem; see Bryan and Schneider (2008), Gotovos et al. (2013), Hu and Ludkovski (2017), Picheny et al. (2010) and Ranjan et al. (2008).In a nutshell, at step n the GP paradigm constructs a metamodel f (n) that is then used to guide the selection of x n+1 and also to construct the estimate S (n) .To this end, GPs are well suited for sequential design by offering a rich uncertainty quantification aspect that can be (analytically) exploited to construct information-theoretic DoE heuristics.The standard framework is to develop an acquisition function I n (x) that quantifies the value of information from a new sample at input x conditional on an existing dataset (x 1:n , y 1:n ) and then to maximize I n : (1.5) Early level set sampling criteria were proposed by Bryan et al. (2006), Vazquez and Martinez (2006), Bichon et al. (2008), Picheny et al. (2010), and Ranjan et al. (2008) based on modifications to the expected improvement criterion (Jones et al. 1998) for response function optimization.A criterion more targeted to reduce the uncertainty about S itself was developed by Bect et al. (2012) using the concept of stepwise uncertainty reduction (SUR).Specifically, the SUR strategy aims to maximize the global learning rate about S; see also Chevalier et al. (2014) for related computational details.Further criteria using tools from random set theory were developed in Chevalier et al. (2013) and Azzimonti et al. (2016) using the notions of Vorob'ev expectation.

Summary of contributions
The main goal of this article is to present a comprehensive assessment of GP-based surrogates for stochastic contour-finding.Learning S requires blending the exploitation objective to locally estimate the contour ∂ S with the exploration objective of visiting less-sampled regions.To do so, we maintain the sequential design paradigm and GP-based surrogates described above which boil down to accurate inference of the mean response and sampling noise that in turn drive the posterior mean f and the posterior GP variance s(x) 2 .Our analysis focuses on the effect of observation noise on contour-finding algorithms and complements Picheny et al. (2013b) and Jalali et al. (2017), who benchmarked GP metamodels for Bayesian optimization (BO) where the objective is to evaluate max x f (x).The latter study observed the strong impact of on performance of BO; in our own analysis, we confirm the need for specialized metamodeling frameworks and sequential design strategies in order to strike the best balance in carrying out uncertainty quantification and constructing a robust surrogate.
Motivated by our application settings, we analyze the joint role of the acquisition function I n (•) and the metamodel f .On the latter front, we seek approaches that are not too swayed by the simulation noise structure.This issue is fundamental to any realistic stochastic simulator where there is no justification for assuming Gaussian-distributed (as opposed to the physical experimental setup where represents observation noise and is expected to be Gaussian thanks to the central limit theorem).This motivates us to study alternative GP-based metamodels for learning S that are resistant to non-Gaussian in (1.2).
First, we investigate two ways to handle heavy-tailed simulation noise: t-observation GPs (Vanhatalo et al. 2009;Jylänki et al. 2011) and Student-t processes (TPs) (Shah et al. 2014;Wang et al. 2017).In particular, we document the strong performance of t-observation GPs; to our knowledge, this is the first use of either tool in sequential analysis and contour-finding.Second, to target the classification-like objective underlying (1.3), we consider the use of classification GPs (Rasmussen and Williams 2006;Williams and Barber 1998) that model the sign of the response Y (x) via a probit model driven by a latent GP Z (•): P(Y (x) > 0|x) = probit(Z (x)).Deployment of the probit regression is expected to "wash out" non-Gaussian features in beyond its effect on the sign of the observations.This context offers an interesting and novel comparison between regression and classification approaches benchmarked against a shared loss function.Third, we present an original use of monotonic GP metamodels (Riihimäki and Vehtari 2010) for level set estimation.This idea exploits a structure commonly encountered in applications where the level set S is connected, suggest-ing to force f to be monotone in the specified coordinates in order to improve upon agnostic black-box strategies.
We then combine all the above metamodels with specialized acquisition functions targeting level set estimation, cf.Sect.3. We consider four choices (Contour Upper Confidence Bound (cUCB), targeted mean squared error (tMSE), stepwise uncertainty reduction (SUR) and gradient stepwise uncertainty reduction (gSUR)) that include heuristics that depend only on the posterior standard deviation s (n) (•), as well as those that anticipate information gain from sampling at x n+1 via the look-ahead standard deviation s (n+1) (•).For the latter, because in the GPs with non-Gaussian noise s (n+1)  depends on Y (x n+1 ), we develop and implement approximate look-ahead formulas s (n+1) for all our metamodels.To our knowledge, this is the first presentation of such formulas for GPs with non-Gaussian noise, as well as TPs.
After setting up the methodological toolbox, we then provide a detailed comparison among the proposed acquisition functions, identifying the best-performing combinations of I(•) and metamodel f and documenting the complex interplay between design geometry and surrogate architecture.To this end, we benchmark across a range of synthetic experiments that aim to stress-test each scheme in terms of model misspecification, in particular in terms of handling non-Gaussian .This benchmarking reveals new insights relative to most of the cited papers that consider only deterministic contour-finding.
The rest of the article is organized as follows.Section 2 describes the metamodels we employ.Section 3 develops the sequential designs for the level set estimation problem.Section 4 compares the models using synthetic data where ground truth is known.Two case studies from derivative pricing are investigated in Sect. 5.In Sect.6, we summarize our conclusions.The look-ahead variance formulas for GPs with non-Gaussian noise are discussed in Section C of Supplementary Material.

Gaussian process regression with Gaussian noise
We begin by discussing regression frameworks for contourfinding that target learning the latent f (•).The Gaussian process paradigm treats f as a random function whose posterior distribution is determined from its prior and the collected samples as a realization of a Gaussian process specified by its mean function m(x) := E[ f (x)] and covariance function In the classical case (Rasmussen and Williams 2006), the noise distribution is homoscedastic Gaussian (x) ∼ N (0, τ 2 ), and the prior mean is zero, m(x) = 0. Given observations y 1:n = [y 1 , . . ., y n ] at inputs x 1:n = [x 1 , . . ., x n ] , the conditional distribution f |A n is then another Gaussian process, with posterior marginal mean f given by (throughout we use subscripts to indicate the metamodel type, e.g., Gsn for Gaussian noise) and The posterior mean f Gsn (x * ) is treated as a point estimate of f (x * ) and the posterior variance s (n)  Gsn (x * ) 2 = v (n) Gsn (x * , x * ) as the uncertainty of this estimate.We use f to denote the random posterior vector f (x 1:n )|A n .
Model Fitting: In this article, we model the covariance between the values of f at two inputs x and x with the squared exponential (SE) function: defined in terms of the hyperparameters ϑ = {σ 2 se , θ 1 , ..., θ d , τ } known as the process variance and lengthscales, respectively.Simulation variance τ is also treated as unknown and part of ϑ.To estimate the hyperparameters ϑ, we use the marginal likelihood p(y 1:n |x 1:n , ϑ) = p(y 1:n |x 1:n , f) p(f|ϑ)df. (2.4) One may similarly express the posterior over the hyperparameters ϑ, where p(y 1:n |x 1:n , ϑ) plays the role of the likelihood.
To avoid expensive MCMC integration, we use the maximum likelihood (ML) estimate ϑ which maximizes the likelihood (2.4).Given the estimated hyperparameters ϑ, we take the posterior of f as p(f|y 1:n , x 1:n , ϑ).The same method is used for all metamodels throughout this article.Remark: In the present article, we target the non-Gaussian aspects, in particular the likely heavy-tailed property of simulation noise.A complementary strand of the literature focuses on heteroscedastic simulation variance; see the stochastic kriging approach of Ankenman et al. (2008) and the earlier works by two of the authors (Binois et al. 2018(Binois et al. , 2019)).

Gaussian process regression with Student t noise
Taking the noise term i as Gaussian is widely used since the marginal likelihood is then analytically tractable.In a stochastic simulation setting, however, the exact distribution of the outputs relative to their mean is unknown and often is clearly non-Gaussian.A more robust choice is to assume that i has a Student-t distribution (Jylänki et al. 2011).In particular, this may work better when the noise is heavy-tailed by making inference more resistant to outliers (O'Hagan 1979).In the resulting t-GP formulation, i is assumed to be tdistributed with variance τ 2 and ν > 2 degrees of freedom.(The latter is treated as another hyperparameter.)Using the Gamma function Γ , the marginal likelihood of observing y 1:n can be written as . (2.5) The likelihood p tGP (y 1:n |x 1:n , f) is no longer Gaussian, and integrating (2.5) against the Gaussian prior p( f |ϑ) is intractable; we therefore use the Laplace approximation (LP) method (Vanhatalo et al. 2009) to calculate the posterior.
A second-order Taylor expansion of log p tGP (f|x 1:n , y 1:n ) around its mode, f(n) tGP := arg max f p tGP (f|x 1:n , y 1:n ), gives a Gaussian approximation where ˝tGP is the Hessian of the negative conditional log posterior density at f(n) tGP : and tGP is diagonal, since the likelihood factorizes over observations.Using (2.6), the approximate posterior distribution is also Gaussian (2.9) Note the similarity to (2.1)-(2.2):With Student-t likelihood, the mode f(n) tGP plays the role of y 1:n and W −1 tGP replaces the noise matrix τ 2 I. Critically, the latter implies that the posterior variance is a function of both designs x 1:n and observations y 1:n .

Gaussian process classification
Our target in (1.1) is to learn where the mean response is positive, which is equivalent to classifying each x ∈ D as belonging either to S or to N .Assuming that This motivates us to consider the alternative of directly modeling the response sign (rather than overall magnitude) via a classification GP model (Cl-GP) (Williams and Barber 1998;Rasmussen and Williams 2006).The idea is to model the probability of a positive observation Y (x) by using a probit regression: P(Y (x) > 0|x) = (Z (x)), with (•) the standard normal cdf.We follow the notations used by Rasmussen and Williams (2006) to formulate the Cl-GP.Other link functions can also be used to model the probability: For example, Williams and Barber (1998) used logistic regression with The latent classifier function is taken as the GP Z ∼ G P(0, K (•, •)).After learning Z , we then set S = {x ∈ D : Z (x) > 0}.
Remark: Cl-GP fundamentally recovers the "medianzero" level set where {x : P(Y (x) ≥ 0) = 0.5}.In cases where the noise is skewed, the recovered Cl-GP would therefore be biased relative to the level set {x : f (x) ≥ 0}.In that situation, use of Cl-GP could be viewed as a bias-variance trade-off.
To compute the posterior distribution of Z conditional on A n , we use the fact that for an observation (x i , y i ) and conditional on z i = Z (x i ), the likelihood of y i > 0 is To simplify notation, we use Y (x) = sign Y (x) ∈ {−1, 1} to represent the signed responses driving Cl-GP, leading to p Cl ( y1:n |z, (2.10) Similar to t-GP, we follow the implementation described in Williams and Barber (1998) and Rasmussen and Williams (2006) and use a Laplace approximation for the non-Gaussian p Cl (z|x 1:n , y1:n ) in Eq. (2.10) (details to be found in Supplementary Material Section B).The posterior mean for Z (•) at x * is then expressed by using the GP predictive mean equation (2.1) and LP approximation (B.1 of Supplementary Material): (2.12) We see the same algebraic structure, with z(n) a stand-in for y 1:n in (2.1) and V −1 a stand-in for τ 2 I in (2.2).

Student-t process regression with Student t noise
Instead of just adding Student-t likelihood to the observations, Shah et al. (2014) proposed t-processes (TPs) as an alternative to GPs, deriving closed-form expressions for the marginal likelihood and posterior distribution of the tprocess by imposing an inverse Wishart process prior over the covariance matrix of a GP model.They found the tprocess to be more robust to model misspecification and to be particularly promising for BO as TPs retain most of the appealing properties of GPs, including analytical expressions, with increased flexibility.
Dealing with noisy observations is less straightforward with TPs, since the sum of two independent Student-t distributions has no closed form.Shah et al. (2014) showed that this drawback can be circumvented by incorporating the noise directly in the kernel.The corresponding data-generating mechanism is taken to be multivariate-t y 1:n ∼ T ν, m(x 1:n ), K + τ 2 I , where the degrees of freedom are ν ∈ (2, ∞).The posterior predictive distribution is then , where (Shah et al. 2014) Comparing with the regular GPs, we have the same posterior mean f Gsn (x * ), but the posterior covariance now depends on observations y 1:n and is inflated: Gsn (x * , x * ).Moreover, the latent function f and the noise are uncorrelated but not independent.Assuming the same hyperparameters, as n goes to infinity, the above predictive distribution becomes Gaussian.
Inference of TPs can be performed similarly as for a GP, for instance, based on the marginal likelihood: (2.15) One issue is estimation of ν, which plays a central role in the TP predictions.It is found in Shah et al. (2014) that restricting ν to be small when maximizing (2.15) is important in order to avoid degenerating to the plain GP setup.

Metamodel performance for level set inference
To evaluate the performance of different metamodels, we consider several metrics.The first statistic is the error rate ER based on the loss function L defined in Eq. ( 1.3), measuring the distance between the level set S and its estimate S: (2.16) The error rate ER evaluates the accuracy of the estimated S when the ground truth is known.In a realistic case study when the latter is unavailable, we replace R by its empirical counterpart, based on quantifying the uncertainty in S through the associated uncertainty of f .Following Azzimonti et al. (2016), we define the integrated posterior error E as the expected distance in measure between the random set S|A and S: (2.17) with local posterior error Ē(x) calculated by using (2.1) and (2.2): (2.18) The local posterior error Ē(x) is the posterior probability of wrongly classifying x conditional on the training dataset A. It is intrinsically tied to the point estimate f (x) and the associated posterior variance s(x) 2 through the Gaussian uncertainty quantification.For TP, the predictive distribution is Student-t, so that the Gaussian cdf is replaced with survival function.For Cl-GP, we replace f (x) with z(x) in (2.18).Uncertainty Quantification: To quantify the overall uncertainty about S (rather than local uncertainty about f (x)), a natural criterion is the volume of the credible band C I α ∂ S that captures inputs x whose sign remains ambiguous given A. Given a credibility level α (e.g., α = 0.05), we define where Noting that the posterior error is Ē(x) = 0.5 for x ∈ ∂ S, we have the approximation D Ē(x)dx We have confirmed empirically that the area of C I α ∂ S is roughly proportional to the integrated posterior error E. We also mention that E is equivalent to the Vorob'ev deviation Chevalier et al. (2013) with a median Vorob'ev threshold.

Sequential design
We estimate the level set S in a sequential design setting that assumes that f is expensive to evaluate, for example, because of the complexity of the underlying stochastic simulator.Therefore, efficient selection of the inputs x 1:n is important.In sequential design, at each step the next sampling location x n+1 is selected given all previous measurements.We follow the standard approach to sequential design that is based on greedily optimizing a posterior-based acquisition function I as in (1.5).These strategies got popularized thanks to the success of the expected improvement (EI) criterion and the associated efficient global optimization (EGO) algorithm (Jones et al. 1998).The basic loop for sequential design is: and sample We now propose several acquisition functions I n (x) in Eq. (1.5).The key plan is to target regions close to the boundary ∂ S. A second strategy is to use the look-ahead posterior standard deviation s (n+1) conditional on sampling at x, in order to assess the corresponding information gain.This links the constructed design to the metamodel for f , since different surrogate architectures quantify uncertainty differently.

Contour upper confidence bound
The first metric, dubbed Contour Upper Confidence Bound (cUCB), stems from the Upper Confidence Bound (UCB) strategies proposed by Srinivas et al. (2012) for Bayesian optimization.The idea of UCB is to express the exploitationexploration trade-off through the posterior mean f (x) and standard deviation s(x).Following the spirit of UCB, cUCB blends the minimization of | f (n) (x)| (exploitation) with maximization of the posterior uncertainty s (n) (x) (exploration): where γ (n) is a step-dependent sequence of weights.Thus, cUCB targets inputs with high uncertainty (large and close to the boundary ∂ S (small | f (n) |), additionally weighted by μ(•).Small γ (n) leads to aggressive sampling concentrated along the estimated ∂ S; large γ (n) leads to space-filling sampling that effectively minimizes the ultimate integrated mean-squared error.Thus, the choice of γ 's is critical for the performance; in particular, γ (n) should be increasing to avoid being trapped in local minima of Based on our experiments (see Supplementary Material Section A), we recommend to adapt γ (n) to the relative ratio between f (n) (x) (for steeper response surfaces, γ should be larger) and s (n) (x).(γ needs to rise as posterior uncertainty decreases.).Our recipe is where Ave(s (n) ) denotes the average of posterior standard deviation and IQR is the inter-quantile range of f (n) over D. This keeps both terms in (3.1) comparable as n changes.
In Supplementary Material Section A, we investigate the performance of cUCB with different γ (n) 's, see Table 1 in Supplementary Material.We observe that (1) there is no sin-gle choice that consistently performs best; (2) the adaptive γ (n) is most frequently the best, achieving the smallest error rate in roughly half the cases.Note that since s (n) decreases in n (signal-to-noise ratio increases over time), (3.2) makes γ (n) increase in n, which is consistent with the theoretical results of Srinivas et al. (2012).
Remark 1 The local posterior error Ē(x) as defined in Eq. (2.18) could be directly used as an acquisition function, i.e., (x) = 0.5 for any x where f (n) (x) = 0), so does not possess a unique maximizer as soon as ∂ S is non-trivial.One potential solution could be to maximize (3.3) over a finite candidate set, which, however, requires significant fine-tuning.

Integrated stepwise uncertainty reduction
In order to take into account the spatial structure of D, we next consider a criterion that targets the global reduction in the uncertainty of S. The integrated Stepwise Uncertainty Reduction metric (SUR), first proposed by Bect et al. (2012), is linked to the posterior error E from Sect.2.5: (3.4) Crucially, I SU R ties the selection of x n+1 to the lookahead mean f (n+1) (x n+1 ) and look-ahead standard deviation s (n+1) (x n+1 ) that appear on the right-hand side of (3.4).
To compute the integral over Y (x), we replace Similarly, we plug in the approximate onestep-ahead standard deviation s (n+1) , cf.Supplementary Material Section C (especially Equations (C.1), (C.13) and (C.15)) for s (n+1) (x).Lastly, we replace the integral over D with a sum over a finite subset D of size M leading to Then, I SU R (x) can be viewed as measuring the overall information gain about S from sampling at x.The motivation behind SUR is to minimize the expected one-step-ahead posterior error E, which would correspond to 1-step Bayesoptimal design.
Remark 2 In the case where i ∼ N (0, τ 2 ), Chevalier ( 2013) calculated an analytical expression for E Y (x) E (n+1) , see Supplementary Material Section E. That computation does not work with non-Gaussian noise, or for Cl-GP and TP.
In Supplementary Material Section E, we compare the two resulting expressions and find that the approximate (3.4) with s (n+1) (x) performs generally better, being less sensitive to the noise structure.

Gradient SUR
Our third strategy focuses on quickly reducing the posterior error by comparing the current Ē(x n+1 ) given A n and the expected Ē(x n+1 ) conditional on the one-step-ahead sample, A n ∪ {x n+1 , y n+1 }.This is conceptually similar to SUR above but considers just the local information gain, dropping the computationally intensive integral over D and only integrating out the effect of Y (x n+1 ) on Ē(x n+1 ): The name gSUR is because (3.6) is related to the knowledge gradient strategy of Frazier et al. (2008), modified to target contour-finding.We apply the same approximation as for SUR to simplify the expectation over Y (x) (3.7) Note that if x is such that f (n) (x) = 0, then both terms above are 1/2 and I gSUR n (x) = 0. Thus, the gSUR criterion will not place samples directly on ∂ S, but will aim to bracket the zero contour.To our knowledge, gSUR is a new criterion for sequential level set estimation that interpolates between the local cUCB and the integrated SUR while still making use of the predictive local gain in information about the contour.As for SUR, one may utilize Equation (E.2) in Supplementary Material for I gSUR under the implicit assumption that the simulation noise is Gaussian.

Targeted mean squared error
As a last alternative, we utilize the targeted mean squared error (tMSE) criterion.tMSE is a localized form of the targeted IMSE criterion in Picheny et al. ( 2010) and is defined as where (3.10) The tMSE criterion upweighs regions close to the zero contour through the weight function W tMSE n (x) which measures the distance of x to ∂ S (n) using the Gaussian posterior density N ( f (n) , s (n) (x) 2 ).Like cUCB, tMSE is based only on the posterior at step n and does not integrate over future Y (x)'s.The tMSE is implemented in the KrigInv R package.Similar to cUCB and gSUR, tMSE picks the input in the region close to the zero contour at step n and favors exploitation, while SUR and tIMSE pick the input that will minimize the integrated one-step-ahead uncertainty over the entire space and is explorative.Computationally, the localized criteria are more efficient since they skip the expensive integration.In this article, we use SUR as a representative of the exploratory criteria and compare its performance with the other localized criteria.
Remark 3 In Picheny et al. (2010), an additional parameter σ was added to the definition of W tMSE n (x) by replacing s (n) (x) everywhere with s (n) (x) 2 + σ 2 .Larger σ yields more space-filling designs as W tMSE n (x) becomes flatter.Since Picheny et al. (2010) dealt with deterministic experiments, σ was necessary to ensure that W tMSE n (x) is well defined at existing x 1:n and the recommendation was for σ to be 5% of the range of f .In our case, s (n) (x) is intrinsically bounded away from zero and (3.10) works well as is.Additional experiments (available upon request) indicate that the performance of (3.9) is not sensitive to σ , so to minimize the number of tuning parameters we stick to σ = 0 in (3.10).
In the TP case, for cUCB, gSUR and SUR, we replace the standard normal cdf (•) appearing in the formulas by its Student-t counterpart (with the estimated degrees of freedom ν n ).For tMSE, to maintain tractability, we keep the same expression (3.10) for the weights W tMSE .

Illustration
For instructive purposes, we consider a one-dimensional case where we use the Gaussian observation GP to learn the sign of the quadratic f (x) = x 2 − 0.75 2 on D = [0, 1], where S = [0, 0.75] and with the unique zero contour at ∂ S = 0.75.The initial design x 1:10 consists of n = 10 inputs drawn according to Latin hypercube sampling (LHS).The observations are Y (x) = f (x) + , where ∼ t 3 (0, 0.1 2 ).In the top plot in Fig. 1, we plot the true f (•), the posterior mean f (100) (•), and the associated 95%-CI.We also show the credible band for ∂ S; in the respective bottom panel, we plot the acquisition functions I cUCB ∂ S for the boundary (2.19).Both cUCB and tMSE select a location very close to the boundary f (n) (x n+1 ) 0. We note that SUR has the flattest acquisition function among all the criteria, consistent with the idea that it will tend to be more explorative.After using the various acquisition functions to select x n+1 at n = 11, . . ., 100, we show in Fig. 2 the resulting designs x 1:n and the final estimate f (100) with a Gaussian observation GP metamodel.As desired, all methods target the true zero contour at ∂ S = 0.75.As a result, the posterior variance s (n) (x) 2 is much lower in this neighborhood; in contrast, especially for tMSE and cUCB, few samples are taken far from x = 0.75, and the posterior uncertainty there remains high.The true zero contour is within the estimated posterior CI for all the criteria.
The bottom row in Fig. 2 shows the sampled location x n as a function of step n.We observe that cUCB and tMSE heavily concentrate their search around the zero contour, leading to few samples (and consequently relatively large posterior errors E (n) ) in other areas, and much wider posterior CIs for these two criteria although the overall error rate ER is comparable.All criteria exhibit an "edge" effect; that is, besides the desired zero contour x = 0.75, multiple samples are taken close to the edges of the input space at x = 0 and x = 1.This occurs due to the relatively large posterior variance s 2 (•) in those regions (which arises intrinsically with any spatiallybased metamodel) that in turn strongly influences I in (3.1), (3.7), (3.5) and (3.9).Inputs sampled by the gSUR criterion bracket the contour ∂ S from both directions, matching the two-hill-and-a-valley shape of I gSUR in Fig. 1.We note that the two sampling "curves" get closer as n grows, indicating a gradual convergence of the estimated zero contour ∂ S (n) , akin to a shrinking credible interval of S (n) .The SUR cri- terion generates a much more diffuse design: It engages in more exploration and is less dependent on the current levels of the posterior error E. This eventually creates a flatter profile for Ē(x).
The preceding discussion considered a single metamodel choice for f .Although Figs. 1 and 2 only present results of one run for each design with plain GP, the features observed are generic.Other metamodels will generate different design features; in particular, sensitivity to will lead to a different mix of exploration (x n 's far from the zero contour) and exploitation even for the same choice of a I n criterion.Figures 6 and 8, as well as Table 3, emphasize our message that one must jointly investigate the combinations of I(•) and f when benchmarking the performance of the algorithm.

Benchmark construction
As synthetic experiments, we consider four benchmark problems in dimension d = 1, 2 and 6.For the latter two, we borrow the widely used Branin-Hoo 2-D, Michalewicz 2-D and Hartman 6-D functions; see, for example, Picheny et al. (2013b).The original functions have been rescaled to map their sample space D onto [0, 1] d ; see Table 1.
The latent functions are chosen to cover a variety of problem properties.The quadratic f in 1-D is strictly monotonically increasing, yielding a single boundary ∂ S. The original Branin-Hoo function (Picheny et al. 2013b) is modified so that f is increasing in x 1 and the zero-level set has a non-trivial shape in x 2 .The Hartman is a multimodal function with a complex zero contour.The parameters in the original Hartman function described, for instance, in Picheny et al. (2013b) are adjusted to reduce the "bumps" in the zero contour and make the problem more appropriate for the sign classification task.These three functions are used to provide a comprehensive comparison between all metamodels.
The Michalewicz function features two intersecting "valleys" surrounded by plateaus and is known to be challenging for GP Bayesian optimization because the spatial correlation is high along the plateaus, but very low (response almost discontinuous) around the local minima.This implies a spatially non-stationary response covariance which makes a Gaussian-noise GP highly misspecified.By construction, Cl-GP focuses only on the sign of f (•) and so should outperform level set estimation for such a case study.We used this function to compare performance of classification and regression GPs.
A large number of factors can influence the performance of metamodels and designs.In line with the stochastic simulation perspective, we concentrate on the impact of the simulation noise and consider four observation setups.These cover a variety of noise distributions and signal-to-noise ratio (SNR), measured through the proportion of standard deviation σ τ to the range R f of the response.The first two settings use Student-t distributed noise, with (i) low σ τ and (ii) high σ τ .The third setting uses (iii) Gaussian mixture noise to further test misspecification of .Lastly, we consider the challenging case of (iv) a heteroscedastic Student-t noise with state-dependent degrees of freedom.In total, we have 3 × 4 × 4 × 6 experiments (indexed by dimensionality, noise setting, design heuristic and metamodel type).
Besides the noise distribution, we fix all other metamodeling aspects.All schemes are initialized with n 0 = 10d inputs drawn from an LHS design on [0, 1] d and use the SE kernel (2.3) for the covariance matrix K .To analyze for the variability due to the initial design and the noise realizations, we 8.00 0.50 3.00 10.00 3.00 8.00 3.50 6.00 10.00 10.00 1.70 0.50 3.50 1.00 8.00 8.00 1.70 6.00 10.00 1.00 6.00 9.00 6.00 9.00 perform 20 macro-runs of each design/acquisition function combination.For each run, the same initial inputs are used across all GP metamodels and designs; however, the initial x 1:n 0 vary across runs.
Optimization of the Improvement Metric: We employed the cUCB, SUR, tMSE and gSUR criteria to maximize the improvement metric I and select the next input x n+1 .This maximization task is non-trivial in higher dimensions because I is frequently multimodal and can be flat around its local maxima.We use a genetic optimization approach as implemented in the ga (genetic algorithm) function in MATLAB using the default parameter settings.(We considered increasing the number of generations with no material impact on results.)This is a global, gradient-free optimizer that uses an evolutionary algorithm to explore the input space D.
Evaluation of Performance Metrics: Recall that evaluating the quality of ∂ S is based on ER and E from (2.16) and (2.17) that require integration over D. In practice, these are computed based on a weighted sum over a finite D, E s(x m ) μ(x m ) for a space-filling sequence D ≡ x 1:M ∈ D of test points.In 1-D experiments, D was an equispaced grid of size M = 1000.In higher dimensions, to avoid the use of a lot of test points that are required to ensure an accurate approximation, we adaptively pick D that targets the critical region close to the zero contour.To do so, we replace the integral with a weighted sum: where M = M 1 + M 2 and the test locations x 1:M 1 and x 1:M 2 are subsampled from a large space-filling (scrambled Sobol) sequence on D. The weight p c determines the relative volume of D 1 and D 2 = D\D 1 , where on D 1 = {x : f (x) 0} we are close to the zero contour.In the experiments below, we use M 1 = 0.8M, M 2 = 0.2M, and p c = 0.4, so that the density of test points close to ∂ S is six times relative to those far from the zero contour.We employ the same strategy for speeding the evaluation of the posterior error E.
Surrogate Inference: Values of hyperparameters ϑ are crucial for a good performance of GP metamodels.We estimate ϑ using maximum likelihood.Except for TP, all models are fitted with the open source package GPstuff (Vanhatalo et al. 2013) in MATLAB.TPs are fitted with the hetGP (Binois et al. 2018) package in R, with procedures matching the MATLAB ones as much as possible.Auxiliary tests did not reveal any significant effects from using other available tools for plain GPs and t-GP, such as GPML (Rasmussen and Nickisch 2010).
In principle, the hyperparameters ϑ change at every step of the sequential design, in other words, whenever A n is augmented with (x n+1 , y n+1 ).To save time, however, we do not update ϑ at each step.Instead, we first estimate the hyperparameters ϑ based on the initial design A n 0 and then freeze them, updating their values only every few steps.Specifically, ϑ is re-estimated at steps n 0 + 1, n 0 + 2, n 0 + 4, n 0 + 8, n 0 + 16, . ... This is done solely for computational efficiency and is driven by the fact that in later stages (n 100), inferred hyperparameters tend to change minimally step to step.We observed minimal impact on accuracy from doing so.
The lengthscales θ i are the most significant for surrogate goodness of fit.A too-small lengthscale will make the estimated f look "wiggly" and might lead to overfitting, while θ i too large will fail to capture an informative shape of the true f and hence S. Since our input domain is always [0, 1] d , we restrict θ i ∈ [0.3, 2] ∀i to be on the order of the length of the sample space D.
Computational Overhead: All the considered metamodels are computationally more demanding than the baseline Gaussian GP.For t-GP and Cl-GP, additional cost arises due to the Laplace approximation.TP necessitates estimation of the parameter ν and also the computation of β in (2.14).
In the experiments considered, the respective computation times were roughly double to triple relative to the Gaussiannoise GP.In terms of sequential design, cUCB, tMSE and gSUR have approximately equal overhead; SUR is significantly more expensive because it requires evaluating the sum in (3.5).Note that all heuristics include two expensive steps: optimization for x n+1 and computation of f (n) and s (n)  (and/or s (n+1) ).
Overall timing of the schemes is complicated because of the combined effects of n (design budget), M (size of test set), and the use of different software (some schemes run in R and others in MATLAB).Most important, the ultimate computation time is driven by the simulation cost of generating Y (x)-samples, which is trivial in the synthetic experiments but assumed to be large in the motivating context.
All of the computer programs that are used to produce results in this section and Sect. 5 are uploaded on Zenodo https://zenodo.org/record/4584456#.

Comparison of regression GP metamodels
Figure 3 shows the boxplots of the error rate ER of S (N )  at the final design (N = 100 in 1-D; N = 150 in 2-D; N = 1000 in 6-D).The plots are sorted by noise settings and design strategies, facilitating comparison between the discussed metamodels.In Table 3, we list the best metamodel and design combination in each case.Several high-level observations can be made.First, we observe the limitations of the baseline Gaussian GP metamodel, which cannot tolerate too much model misspecification.As the noise structure gets more complex, the classical GP surrogate begins to show increasing strain; in the last t/hetero setup, it is both unstable (widely varying performance across runs) and inaccurate, with error rates upward of 30% on "bad" runs.In addition, according to results shown in Table 3, across all of the twelve cases, besides the 1-D experiment with t/small noise, the Gaussian-noise GP never performs as the best metamodel.This result is not surprising but confirms that the noise distribution is key for the contour-finding task and illustrates the non-robustness of the Gaussian observation model, due to which outliers strongly influence the inference.
Second, we document that the simple adjustment of using Student-t observations significantly mitigates the above issue.t-GP performs consistently and significantly better than Gaussian-noise GP in essentially all settings.This result is true even when both models are misspecified (the Gsn/mi x and t/hetero cases).The performance of t-GP was still better (though not statistically significantly so) when we tested it in the setting of homoscedastic Gaussian noise (see Table 3 in Supplementary Material Section F).The latter fact is not surprising-t-GP adaptively learns the degreesof-freedom parameter ν and hence can "detect" Gaussian noise by setting ν to be large.Conversely, in heavy-tailed noise cases, the use of Student-t likelihood will effectively ignore outliers (O'Hagan 1979) and thus produce more accurate predictions than working with a Gaussian observation assumption.We find that t-GP can handle complex noise structures and offers a good choice for all-around performance, making it a good default selection for applications.It brings smaller error rate ER, more stable hyperparameter estimation, less contour bias and tighter contour CI.Moreover, t-GP is significantly better than all the other GPs in eight of the twelve setups, indicating that t-GP is essentially the best out of all GP metamodels in most cases.Figure 4 shows the mean error rate ER (2.16) as a function of step n in the 6-D t/large and t/hetero experiments.This illustrates the learning rates of different models and schemes as data are collected.In both cases, t-GP starts from more accu-Fig.3 Boxplots of final error rate ER (n)   Third, we also inspect the performance of the TP metamodel.As shown in Table 3, TP is the best in two cases out of the twelve, both of which are with the t/small noise.We note that TP works worst in t/hetero cases, having both large error rate ER and posterior error E. Therefore, TP does not work well in cases with low SNR or greatly misspecified noise.This may be related to the parameterization of TPs, with noise handled in the kernel, which seems less robust to misspecification.Also, since TPs revert to GPs as n increases, the advantage of flexibility offered by the modeling decreases as iterations goes and thus may not last enough for low SNRs, which require more samples.It is apparent, for instance, in Fig. 4, where the learning rate at early stage (step < 250) is larger than for its counterparts.
Table 3 shows that there is no one overall "best" design for all metamodels across all cases.However, it does suggest some design/metamodel "combos" that work better than others.The classification GPs seem to prefer more aggressive designs, such as cUCB and tMSE, while t-GP prefers more exploratory designs, such as SUR, especially in higher dimension.Additional discussion regarding the relative performance of design heuristics is in Sect.4.4.

Classification GPs for level set estimation
We generally find that classification GPs are not the best metamodel.Namely, Cl-GP with cUCB design has the smallest error rate only in one (t/hetero in 1-D) out of 12 cases shown in Table 3.However, Cl-GP is often competitive and is better than Gaussian-noise GP in some cases with tMSE and cUCB designs (except for the 6-D cases, where the error rate ER of cUCB is not significantly different from that of SUR, although mean of SUR is slightly smaller).There is significant improvement for models with low SNR; the only exception is for the low-noise setup where Cl-GP underperforms baseline GP. Figure 4 shows that in the 6-D t/hetero case, Cl-GP achieves error rate ER which is two-thirds of ER for GP and TP at step 0. Also, it enjoys faster reduction in ER (n) , while the design size n is small.This matches the intuition that employing classification "flattens" the signal by removing outliers.By considering only the sign of the response, the classification model disregards its magnitude, simplifying the noise structure at the cost of some information loss.The net effect is helpful when the noise is misspecified or too strong so as to interfere with learning the mean response.It is detrimental if the above gain is outweighed by the information loss, as apparently happens in the 6-D experiments when the design size n is large.
We also observe that the stability of Cl-GP is highly dependent on the design: Some designs create large across-run variations in performance.We hypothesize that this is due to a more complex procedure for learning the hyperparameters of Cl-GP; therefore, designs that are not aggressive enough to explore the zero contour region (such as gSUR) face difficulties in estimating ϑ.As a result, relative to t-GP, Cl-GP has larger sampling variances.
To provide better intuition regarding the settings where Cl-GP works best, we test the performance of classification GP and regression GP on the 2-D modified Michalewicz function with Gsn/mi x noise.Table 4 and Fig. 5 present the results comparing Cl-GP with cUCB against t-GP with SUR.The latter scheme was chosen as being the typical bestperforming combination.We evaluate ER with budgets of n = 150 and n = 500.In both cases, Cl-GP performs better than t-GP, achieving significantly smaller error rate ER.This is visualized in Fig. 5 that shows the estimated contours ∂ S at n = 500.While both models detect the vertical zero contour, only Cl-GP identifies the horizontal boundary and samples inputs in that neighborhood.Figure 5 also shows the error rate ER as a function of step for Cl-GP and t-GP.Cl-GP starts with a smaller ER and higher learning rate for both cUCB and SUR in the early stage (step < 200) compared Results are based on 20 macro-replications of each scheme with t-GP.Cl-GP with cUCB ends up with the most accurate estimate, echoing the results shown in Table 4.
In conclusion, the modified Michalewicz function illustrates the potential gain of using Cl-GP over regression GPs where the covariance is not spatially stationary.The nonstationary space also benefits the localized criteria over the global ones, since the global change might be neutralized and may not provide an exact measure of uncertainty.This also explains why Cl-GP performs significantly better with cUCB than SUR.

Designs for contour-finding
A key goal of our study is qualitative insights about experimental designs most appropriate for noisy level set estimation.Through identifying the best-performing heuristics, we get an inkling regarding the structure of near-optimal designs for (1.1).In this section, we illustrate the latter within a 2-D setup that can be conveniently visualized.Taking the t/large experiment as an example, in Fig. 6 we plot the fitted zero contour ∂ S at N = 150 together with the chosen inputs x 1:150 across the six metamodels and the four I heuristics.As expected, most of the designs are around the contour ∂ S, which is the intuitive approach to minimize the error ER.Nevertheless, we observe significant differences in designs produced by different I's.The cUCB criterion places most of the samples close to the estimated zero contour ∂ S, reflecting its aggressive exploitation nature.It leads to little information collected about other regions than the contour ∂ S, especially around the boundaries of sample space D and hence relatively large Ē(x) there, inflating E for cUCB.As visualized in Fig. 7, the posterior error E tends to be the largest for cUCB.For tMSE, the samples tend to cluster at several subregions of ∂ S and on the edges of D. For gSUR, x 1:n cover a band along ∂ S, resembling the shape of the cUCB design but more dispersed.Both approaches are better at reducing E compared with cUCB but are not directly aimed at this.For SUR, the design is much more exploratory, covering a large swath of D. All these findings echo the 1-D setting in Fig. 2.  Another noteworthy feature is replication of some inputs, that is, repeated selection of the same x sites.This does not occur for cUCB, but happens for SUR, tMSE and gSUR that frequently (across algorithm runs) sample repeatedly at the vertices of D (indicated by the size of the corresponding marker in Fig. 6).The replication is typically mild.(We observe 145+ unique designs among a total of 150 x n 's.)This finding echoes (Binois et al. 2019) on the importance of replication to distinguish between signal and noise, which is a key distinction with the noise-free setting ≡ 0.
Given the above discussion and the relative overhead of the different heuristics, we conclude that in lower-dimensional problems, there is little benefit to using the more sophisticated SUR criterion, while for higher-dimensional problems, SUR criterion is significantly better than the others.Beyond that, tMSE appears to be an adequate and cheaper choice.However, as the input space becomes more complicated (with a higher dimension or lower SNR observations), we need more exploration over the input space and the explorative criteria like SUR start to shine.

Application to optimal stopping problems in finance
In our next case study, we consider contour-finding for determining the optimal exercise policy of a Bermudan financial derivative.The underlying simulator is based on a d-dimensional geometric Brownian motion (X t ) that rep-resents prices of d assets with interest rate r , dividend yield δ and volatility σ and follows the log-normal dynamics where ΔW t ∼ N (0, ΔtΣ) i.i.d.across t, and Σ is a given d × d covariance matrix.Let h(t, x) be the option payoff from exercising when X t = x ∈ R d .The Bermudan option pricing problem consists of maximizing the expected reward h(τ, X τ ) over all stopping times τ ∈ {0, Δt, 2Δt, . . ., T } (exercising is allowed every Δt time units) bounded by the specified horizon T : (5.2) The approach in the so-called regression Monte Carlo (RMC) methods (Longstaff and Schwartz 2001;Tsitsiklis and Van Roy 2001) is to convert the decision of whether to exercise the option τ (t, x) = t or continue τ (t, x) > t, when X t = x at intermediate step t, into comparing the immediate reward h(t, x) vis à vis the reward-to-go C(t, x).The optimal strategy is to dynamically pick the action with the higher expected payoff and is equivalent to determining the zero level set (known as the stopping region) S t = {x ∈ D : f (x; t) ≤ 0} of the timing value f (x; t) := C(t, x)−h (t, x).
During backward dynamic programming, we iterate over t = T , T − Δt, . . ., 0, recursively estimating the respective S t .To do so, we simulate trajectories X x t:T emanating from input x; the simulator Y (x; t) returns the difference between the pathwise payoff based on the exercise strategy summarized by the forward-looking { S s , s > t} along the trajectory (X x t:T ), and h(t, x).
We refer to Ludkovski (2018) for the full details of employing a GP metamodel for learning the timing value f (•; t); as noted there this setting implies a skewed, non-Gaussian, heteroscedastic distribution of the simulation noise and is a challenging stochastic contour-finding problem.Note that in order to reflect the underlying distribution of X t at time t (conditional on the given initial value X 0 = x 0 ), the weighting measure μ(x) = p X t (x|x 0 ) is used.Thus, μ(•) is log-normal based on (5.1).In line with the problem context, we no longer directly measure the accuracy of learning {S t } but instead focus on the ultimate output of RMC, which is the estimated option value in (5.2).The latter must itself be numerically evaluated via an out-of-sample Monte Carlo simulation that averages realized payoffs along a large database of M paths x 1:M 0:T : (5.3) Since our goal is to find the best exercise strategy, higher V 's indicate a better approximation of the ground truth {S t }.
To allow a direct comparison, we set parameters matching the test cases in Ludkovski (2018), considering a 2-D and 3-D example.In both cases, the volatility matrix Σ = σ I in (5.1) is diagonal with constant terms; that is, the coordinates X 1 1:n , . . ., X d 1:n are independently and identically distributed.As a first example, we consider a 2-D basket Put option with parameters r = 0.06, σ = 0.2, δ = 0, Δt = 0.04, K = 40, T = 1.The payoff is h(t, x) = e −rt (K − x 1 +x 2 2 ) + with K = 40.Here, it is known that stopping becomes optimal once both asset prices x 1 and x 2 become sufficiently low, so the level set S t is toward the bottom-left of D; see Fig. 8.In contrast, stopping is definitely suboptimal when h(t, x) = 0 ⇔ (x 1 +x 2 )/2 > K. Consequently, the input sample space is taken to be In this first case study, the timing value f (x; t) is known to be monotonically increasing in the asset price x.To incorporate this constraint, we augment the four main metamodels (GP, t-GP, Cl-GP and TP) with two monotonic versions, M-GP and MCl-GP.By constraining the fitted f to be monotone, we incorporate structural knowledge about the ground truth, which in turn reduces posterior uncertainty and thus might produce more accurate estimates of S. Monotonicity of the metamodel for f is also one sufficient way to guarantee that the outputted level set S is a connected subset of D.
Our monotone GPs are based on Riihimäki and Vehtari (2010).In general, any infinite-dimensional Gaussian pro-cess is intrinsically non-monotone, since the multivariate Gaussian distribution is always supported on the entire R d , rather than an orthant.Nevertheless, local monotonicity in f can be enforced by considering the gradient ∇ f of f which is also a Gaussian process.Specifically, Riihimäki and Vehtari (2010) proposed to adaptively add virtual observations for ∇ f ; we employ the resulting implementation in the public GPstuff library (Vanhatalo et al. 2013) to build our own version dubbed M-GP.We employ the same strategy to restrict the coordinates z j of the latent probit GP Z to be increasing (decreasing) across D. Implementation details are included in Supplementary Material Section D.
Because of the iterative construction of the simulator, the SNR gets low for small t's.The variance τ 2 (x) is also highly state dependent, tending to be smaller for sites further from the zero contour.While in this case study, the simulations are very fast, low SNR requires many hundreds of observations to reliably detect the level set.Since the expense of sequential design of GP metamodels comes mainly from choosing the new input at each step, it is impractical to have such a large n.To reduce metamodel overhead, we employ batched designs (Ludkovski 2018;Ankenman et al. 2008), reusing x ∈ D for r replications to collect observations y (1) (x), . . ., y (r ) (x) from the corresponding simulator Y (x).Then, we treat the mean of the r observations, as the response for input x and use (x, ȳ(x)) as a single design entry.Such reduction by a factor of r in the number of unique inputs n = N /r significantly speeds fitting and updating.Moreover, the statistical properties of ȳ are improved thanks to the central limit theorem (CLT): Noise variance τ 2 (x) = τ 2 (x)/r is much smaller, and its distribution is more Gaussian.In the case studies below, we observe raw skewness in the range of [−2, 1] and raw kurtosis in the range of [2,20]; after batching skewness decreases by a factor of √ r and excess kurtosis mostly disappears for r ≥ 10.Nevertheless, unless r is in the dozens, the distribution of ȳ For the 2-D Put case study, we then test a total of three budget settings: (i) r = 3, n = 80 (low budget of N = 240 simulations); (ii) r = 15, n = 80 (high budget N = 800 with moderate replication); and (iii) r = 48, n = 25 (high N = 800 with high replication).Comparing (ii) and (iii) shows the competing effects of having non-Gaussian noise (for lower r ) and small design size (low n).The initial design size is n 0 = 10.In this example, taking n 80 gives only marginally better performance but significantly raises the computation time and hence is ruled out as impractical.Three setups are investigated for the 3-D example: r = 3, n = 100 (lowbudget of N = 300), r = 20, n = 100 (moderate-budget of N = 2000) and r = 20, n = 200 (high budget N = 4000), both with n 0 = 30.In all examples, the results are based on 25 runs of each scheme and are evaluated through the resulting expected reward V (0, x 0 ) (5.3) on a fixed out-ofsample testing set of M = 160,000 paths of X 0:T .
The results for Gaussian-noise GP, t-GP and TP can be reproduced via the publicly available mlOSP R package at http://github.com/mludkov/mlOSP, which in particular implements all the discussed acquisition functions.

Results
Tables 5 and 6 compare the different designs and metamodels.To assess the sequential design gains, we also report the results from using a baseline non-adaptive LHS design on D. At low budget, we observe the dramatic gains of using adaptive designs for level set estimation, which allow us to obtain the same performance with an order-of-magnitude smaller simulation budget.The tMSE and gSUR criteria work best for the 2-D Put, while SUR is the best for the 3-D max-Call, indicating that the exploratory designs start to win out in more complex settings with higher d.
Regarding the metamodels, in the low-budget setups, the monotonic GP metamodel works best for the 2-D Put and t-GP for the 3-D max-Call.For the higher budget, which also coincides with higher r ∈ {10, 50}, the metamodel performance is similar, with t-GP slightly better than the other GP variants.In particular, once the SNR is high, classical Gaussian GP is effectively as good as any alternative.
In both examples, TP and classification metamodels do not work well, possibly because of being more sensitive to the heteroscedastic aspect.We note that TP as well as the classification metamodels suffers from instability, so that lower V (0, x 0 ) is matched with a high sampling standard devia-  Figure 8 shows the estimated exercise boundary ∂ S t with its 95% CI at t = 0.4 for the 2-D Put, for each of the five metamodels, each with the design yielding the highest payoff.We observe that all the best-performing designs look similar, placing about a dozen x n 's (some of which are from the initial design x 1:n 0 ) throughout D and the rest tightly along the zero contour.The results suggest that the criteria are largely interchangeable and that simpler I n heuristics are able to reproduce the features of the more sophisticated or expensive SUR.The heuristics do differ in their uncertainty quantification; t-GP and GP generate tightest CI bands, while those of classification GPs and TP are too wide, indicating lack of confidence in the estimate.Of note, the regression GP metamodels (GP, t-GP and M-GP) also generate the lowest sampling variance for V (0, x 0 ).
Based on these results, our take-aways are threefold.First, similar to Ludkovski (2018) we document significant gains from sequential design.Second, we find that while using SUR is helpful for more complicated settings with higher dimension d and larger budget, tMSE is the recommended DoE heuristic for lower-dimensional cases, achieving excellent results with minimal overhead (in particular without requiring look-ahead variance).Third, we find that for applications with thousands of simulations, the Gaussian observation model is sufficient, since the underlying design needs to be replicated r 1 in order to avoid excessively large K-matrices.Therefore, there is little need for more sophisticated metamodels, although useful gains can be realized from enforcing the monotonic structure, if available.

Conclusion
We have carried a comprehensive comparison of five metamodels and four design heuristics on 19 case studies (4×3+1 synthetic, plus six real worlds).In sum, the considered alternatives to standard Gaussian-observation GP do perform somewhat better.In particular, t-GP directly nests plain GP and hence essentially always matches or exceeds the performance of the latter.We also observe gains from using Cl-GP when SNR is low or the response covariance is not spatially stationary and from monotonic surrogates when the underlying response is monotone.That being said, final recommendation regarding the associated benefit depends on computational considerations, as the respective overhead becomes larger (and exact updating of the metamodel no longer possible).
In terms of design, we advocate the benefits of tMSE in low-dimensional simulations, which generates high-performing experimental designs without requiring expensive acquisition function (or even look-ahead variance).The tMSE criterion does sometimes suffer from the tendency to put many designs at the edge of the input space but otherwise tends to match the performance of more complex and computationally intensive I n 's.For complex simulations, SUR is probably still the best choice (although in that case, randomset-based heuristics should also be considered).Particularly in higher dimensions with misspecified noise, SUR is the best choice among all designs for t-GP.We also stress that the user ought to thoughtfully pick the combination of sequential design and metamodel, since cross-dependencies are involved (e.g., classification metamodels generally do not work well with the SUR criterion in lower dimension).
In the original application to BO,Srinivas et al. (2012) proposed γ(n) = C log n and showed that for a certain choice of C, one can then control with high probability the respective cumulative regret.A constant choice of γ(n) = 1.96 corresponds to the Straddle scheme inBryan et al. (2006) and leads to I n (x) ≥ 0 ⇔ x ∈ C I 0.95 (∂ S).Gotovos et al. (2013) employed γ(n) = 3 for a confidence region instead of confidence interval as in Straddle and (3.1), andBogunovic et al. (2016) suggested γ(n) = log(|D|n 2 ); both papers mention that the recommendation inSrinivas et al. (2012) is too conservative and tends to over-exploration.
as defined in Eqs.(3.1), (3.7), (3.5), and (3.9).Comparing the acquisition functions of the four criteria, we find that, besides SUR, all of the others have maxima within the shaded credible interval of the boundary C I 0.95 ∂ S .In practice, we care only about the maximizer of the acquisition function, rather than its full shape, since the former drives the selection of the next sample x n+1 .The x n+1 's selected by cUCB and tMSE criteria are close.For the gSUR criterion, because I gSUR n (x) = 0 at ∂ S, there are two local maxima with a "valley" between them.The interval between the two local maxima is roughly the credible interval C I0.95

Fig. 1
Fig. 1 Comparison of acquisition functions.Upper panel: true function f = (x + 0.75)(x − 0.75) (black dashed line), the posterior mean f (•) (solid line) and 95% C I α f (shaded area) based on observed samples (x 1:100 , y 1:100 ) (blue dots).Along the x-axis, we also show the credible interval of the partition boundary C I α ∂ S (gray solid line) relative the true zero contour S = [0, 0.75] (red triangle).Lower panel: acquisition functions I n (•) for cUCB, gSUR, SUR and tMSE criteria, with vertical lines marking the respective maxima arg max x I n (x)

Fig. 2
Fig. 2 Top row: Fitted metamodel f (100) (solid red line) and its 95%-CI (shaded region) versus the true f = (x + 0.75)(x − 0.75) (dashed black), for each of the four design strategies.The estimated 95% CI for the zero contour ∂ S is marked on the x-axis with a gray interval; red triangle indicates the true zero contour ∂ S = 0.75.Bottom row: sampled inputs x n (on the x-axis to match the top row) as a function of step n = 1, . . ., 100 (on the y-axis, moving from top to bottom) for cUCB, tMSE, gSUR and SUR criteria.The rug plots at the bottom visualize the overall distribution of x 1:n at n = 100.The first ten inputs are selected using a (fixed-across schemes) LHS design on D = [0, 1]

Fig. 5 a
Fig. 5 a Michalewicz function; b error rate ER (n) (2.16) as a function of step n in the 2-D modified Michalewicz experiments; c and d: The estimated zero contour ∂ S (red solid line) and its 95% credible band (dot dashed lines) for Cl-GP and t-GP at n = 500.Blue dots are the samples selected by cUCB.Black dashed line is the true zero contour ∂ S

Fig. 8
Fig. 8 The estimated exercise boundary ∂ S (solid line with 95% CI as dashed lines) at t = 0.6 for 2-D Bermudan Put from Sect. 5. Shading, which varies panel to panel, indicates the point estimate for the latent f (x) or z(x).We also show the design (x 1:n , y 1:n ) with positive y n 's

Table 3
Mean (w/standard deviation) error rate ER and corresponding best-performing sequential design heuristic for the 1-D exponential function, 2-D Modified Branin-Hoo function and 6-D Modified Hartman6 function synthetic case studies Results are based on 20 macro-replications of each scheme.Best combinations for each column are indicated in bold especially in t/hetero case where the initial error rate ER for t-GP is approximately two-thirds of that for GP.It learns the contour of interest over iterations, has a fast decreasing error rate ER and ends up with the smallest error rate in both cases.

Table 4
Mean (w/standard deviation) error rate ER for Cl-GP and t-GP with cUCB and SUR in 2-D modified Michalewicz synthetic experiments

Table 5
Results are the mean (standard deviation) payoff of 25 runs of experiments evaluating on the same out-ofsample testing set of M = 160,000 X 0:T -paths at each run.Best combinations indicated in bold.Best designs highlighted in gray remains heteroscedastic and highly state dependent in both x and t, making a Gaussian-noise GP strongly misspecified.

Table 6
Performance of different designs and models on the 3-D Bermudan max-Call in