**Lemma 36** Using the notations from the proof of Proposition 18 and writing $v_n = \sup_{x \in \mathcal{X}} \rho_{n,t_n}(x)$, we have $\liminf v_n = 0$.

**Proof.** This is an adaptation of Lemma 12 from Vazquez and Bect (2010b).

Let $y^*$ be a cluster point of $(x_n)_{n \ge 1}$ and let $(x_{\phi(n)})_{n \ge 1}$ be a subsequence converging to $y^*$. We are going to prove that $v_{\phi(n)-1} \to 0$. Define

$$z_{\phi(n)-1} = m_{\phi(n)-1} - s_{\phi(n)-1}(x_{\phi(n)}).$$

Then, Lemma 35 gives<fnref n="1" /> $\liminf s_{\phi(n)-1}(x_{\phi(n)}) \ge \min(f(y^*), t_\infty)$. Moreover we have

$$m_{\phi(n)-1} \le \min(f(x_{\phi(n-1)}), t_{\phi(n-1)}),$$

since $\phi(n-1) \le \phi(n) - 1$, so $\lim m_{\phi(n)-1} \le \min(f(y^*), t_\infty)$ holds because $(m_{\phi(n)-1})_{n \ge 1}$ is non-increasing. The previous arguments show that $\limsup z_{\phi(n)-1} \le 0$.

Moreover, one can use Proposition 10 from Vazquez and Bect (2010b, (i) ⇒ (ii)) similarly to show that $\sigma^2_{\phi(n)-1}(x_{\phi(n)}) \to 0$ and therefore

$$v_{\phi(n)-1} = \gamma(z_{\phi(n)-1}, \sigma^2_{\phi(n)-1}(x_{\phi(n)})) \le \gamma\left(\sup_{k \ge n} z_{\phi(k)-1}, \sigma^2_{\phi(n)-1}(x_{\phi(n)})\right) \to 0,$$

since $\gamma$ is non-decreasing with respect to its first argument and continuous.

<figure><img src="image_17.png" /></figure>

**Proof of Proposition 18.** This is an adaptation of Theorem 6 from Vazquez and Bect (2010b). Write $s_n = s_{R_n, n}$ for the reGP predictor at step $n$ to avoid cumbersome notations. Then, for $x \in \mathcal{X}$, write $\rho_{n,t_n}(x) = \gamma(m_n - s_n(x), \sigma_n^2(x))$ for the expected improvement under the reGP predictive distribution, with $\gamma$ the function defined in Proposition 1.

Suppose that there exists some $x_0 \in \mathcal{X}$ such that $\sigma_n^2(x_0) \ge C_1 > 0$. The sequence $(m_n)_{n \ge 1}$ converges and the reproducing property (6) yields

$$|s_n(x_0)| \le \sqrt{k(x_0, x_0)} \|s_n\|_{\mathcal{H}(\mathcal{X})} \le \sqrt{k(x_0, x_0)} \|f\|_{\mathcal{H}(\mathcal{X})},$$

so the sequence $\|m_n - s_n(x_0)\|_{n \ge 1}$ is bounded by, say $C_2$. We then have

$$v_n = \sup_{x \in \mathcal{X}} \rho_{n,t_n}(x) \ge \gamma(m_n - s_n(x_0), \sigma_n^2(x_0)) \ge \gamma(-C_2, C_1) > 0$$

by Proposition 1. But this yields a contradiction with Lemma 36, so the decreasing sequence $(\sigma_n^2)_{n \ge 1}$ converges pointwise on $\mathcal{X}$ to zero. Proposition 10 from Vazquez and Bect (2010b) then implies that every $x \in \mathcal{X}$ is adherent to $\{x_n\}$. ■

**Lemma 37** Assume that $b$ is finite, and that either $a$ is finite too or $\int_{-\infty}^0 F_P(u) du = \int_{-\infty}^0 |u| P(du)$ is finite. Then

$$S_{a,b}^{\text{CRPS}}(P, z) = (b - a \vee z)_+ + R_2(a, b) - 2R_1(a \vee z, b),$$

where

$$R_q(a, b) = \int_{-\infty}^{+\infty} 1_{a \le u \le b} F_P(u)^q du.$$

<footnote n="1">1. Lemma 35 yields $\liminf s_{n-1}(x_{\psi(n)}) \ge \min(f(y^*), t_\infty)$ with $\psi(n) = \max\{\phi(k), \phi(k) \le n\}$, and the claim follows by extracting a $\phi$-subsequence.</footnote>
C. Chevalier and D. Ginsbourger. Fast computation of the multi-points expected improvement with applications in batch selection. In *International Conference on Learning and Intelligent Optimization*, pages 59–69. Springer, 2013.

C. Chevalier, J. Bect, D. Ginsbourger, E. Vazquez, V. Picheny, and Y. Richet. Fast parallel kriging-based stepwise uncertainty reduction with application to the identification of an excursion set. *Technometrics*, 56(4):455–465, 2014.

W. Chu and Z. Ghahramani. Gaussian processes for ordinal regression. *Journal of machine learning research*, 6(7), 2005.

S. Da Veiga and A. Marrel. Gaussian process modeling with inequality constraints. In *Annales de la Faculté des sciences de Toulouse: Mathématiques*, volume 21, pages 529–555, 2012.

A. Damianou and N. D. Lawrence. Deep Gaussian processes. In *Artificial intelligence and statistics*, pages 207–215. PMLR, 2013.

E. Di Nezza, G. Palatucci, and E. Valdinoci. Hitchhiker's guide to the fractional Sobolev spaces. *Bulletin des sciences mathématiques*, 136(5):521–573, 2012.

O. Dubrule. Cross validation of kriging in a unique neighborhood. *Journal of the International Association for Mathematical Geology*, 15:687–699, 1983.

M. M. Dunlop, M. A. Girolami, A. M. Stuart, and A. L. Teckentrup. How deep are deep Gaussian processes? *Journal of Machine Learning Research*, 19(54):1–46, 2018.

D. Eriksson, M. Pearce, J. R. Gardner, R. C. Turner, and M. Poloczek. Scalable global optimization via local Bayesian optimization. In *NeurIPS*, 2019.

P. Feliot, J. Bect, and E. Vazquez. A Bayesian approach to constrained single-and multi-objective optimization. *Journal of Global Optimization*, 67(1-2):97–133, 2017.

P. I. Frazier, W. B. Powell, and S. Dayanik. A knowledge-gradient policy for sequential information collection. *SIAM Journal on Control and Optimization*, 47(5):2410–2439, 2008.

M. N. Gibbs. *Bayesian Gaussian processes for regression and classification*. PhD thesis, Citeseer, 1998.

T. G. Gneiting and A. E. Raftery. Strictly proper scoring rules, prediction, and estimation. *Journal of the American Statistical Association*, 102:359–378, 2007.

T. G. Gneiting and R. Ranjan. Comparing density forecasts using threshold-and quantile-weighted scoring rules. *Journal of Business & Economic Statistics*, 29(3):411–422, 2011.

T. G. Gneiting, A. E. Raftery, A. H. Westveld III, and T. Goldman. Calibrated probabilistic forecasting using ensemble model output statistics and minimum crps estimation. *Monthly Weather Review*, 133(5):1098–1118, 2005.
**Proof.**

$$
\begin{align*}
S_{a,b}^{\text{tCRPS}}(P, z) &= \int_{-\infty}^{+\infty} 1_{a \le u \le b} (1_{z \le u} - F_P(u))^2 \, du \\
&= \int_{-\infty}^{+\infty} 1_{a \le u \le b} (1_{z \le u} + F_P(u)^2 - 21_{z \le u} F_P(u)) \, du \\
&= \int_{-\infty}^{+\infty} (1_{a \lor z \le u \le b} + 1_{a \le u \le b} F_P(u)^2 - 21_{a \lor z \le u \le b} F_P(u)) \, du \\
&= (b - a \lor z)_+ + R_2(a,b) - 2R_1(a \lor z, b).
\end{align*}
$$

**Lemma 38** Let $a,b \in \mathbb{R}$ with $a \le b$. Let $q \ge 1$. Then

$$
R_q(a,b) = b - a + \mathrm{EI}_q^{\uparrow}(P,b) - \mathrm{EI}_q^{\uparrow}(P,a).
$$

**Proof.**

$$
\begin{align*}
R_q(a,b) &= \int 1_{a \le u \le b} F_P(u)^q \, du \\
&= \int 1_{a \le u \le b} \prod_{j=1}^{q} E(1_{N_j \le u}) \, du \quad \text{with } N_j \stackrel{\text{iid}}{\sim} P \\
&= E\left(\int 1_{a \lor N_1 \lor \cdots \lor N_q \le u \le b} \, du\right) \\
&= E\left((b - a \lor N_1 \lor \cdots \lor N_q)_+\right) \\
&= b - a + \mathrm{EI}_q^{\uparrow}(P,b) - \mathrm{EI}_q^{\uparrow}(P,a).
\end{align*}
$$

**Proof of Proposition 22.** The first result is given by Lemma 37 and Lemma 38.

Then using the dominated convergence theorem, it is easy to see that, when $a \to -\infty$

$$
\mathrm{EI}_q^\uparrow(P, a) = \mathrm{E}(N_1 \vee \cdots \vee N_q) - a + o(1),
$$

and therefore

$$
R_q(-\infty, b) = \lim_{a \to -\infty} R_q(a, b) = b + \mathrm{EI}_q^{\uparrow}(P, b) - \mathrm{E}(N_1 \vee \cdots \vee N_q),
$$

which gives the second statement.

Finally, a change of variable gives

$$
\int_a^{+\infty} (F_P(u) - 1_{z \le u})^2 du = \int_{-\infty}^{-a} (F_P(u) - P(U = -u) - 1_{-z < u})^2 du,
$$

and the last statement follows by observing that a probability measure admits at most a countable
number of atoms.

■
R. B. Gramacy and H. K. H. Lee. Bayesian treed Gaussian process models with an application to computer modeling. *Journal of the American Statistical Association*, 103(483):1119–1130, 2008.

P. Grisvard. *Elliptic problems in nonsmooth domains*. Pitman, 1985.

A. Hebbal, L. Brevault, and N. Melab. Bayesian optimization using deep Gaussian processes with applications to aerospace system design. *Optimization and Engineering*, 22(1):321–361, 2021.

D. Higdon. A process-convolution approach to modelling temperatures in the north atlantic ocean. *Environmental and Ecological Statistics*, 5:173–190, 1998.

D. Higdon. Space and space-time modeling using process convolutions. *Quantitative Methods for Current Environmental Issues*, 04 2002.

K. Jakkala. Deep Gaussian processes: A survey. arXiv preprint arXiv:2106.12135, 2021.

D. Jones, M. Schonlau, and W. Welch. Efficient global optimization of expensive black-box func-
tions. *Journal of Global Optimization*, 13:455–492, 12 1998.

M. Lázaro-Gredilla. Bayesian warped Gaussian processes. *Advances in Neural Information Processing Systems*, 25:1619–1627, 2012.

S. Lerch and T. L. Thorarinsdottir. Comparison of non-homogeneous regression models for probabilistic wind speed forecasting. *Tellus A: Dynamic Meteorology and Oceanography*, 65(1):21206, 2013.

J. R. Lewis, S. N. MacEachern, and Y. Lee. Bayesian restricted likelihood methods: Conditioning on insufficient statistics in Bayesian regression. *Bayesian Analysis*, 1(1), 2021.

A. F. López-Lopera, F. Bachoc, N. Durrande, and O. Roustant. Finite-dimensional Gaussian approximation with linear inequality constraints. *SIAM/ASA Journal on Uncertainty Quantification*, 6 (3):1224–1255, 2018.

H. Maatouk and X. Bay. Gaussian process emulators for computer experiments with inequality constraints. *Mathematical Geosciences*, 49(5):557–582, 2017.

S. Marmin, D. Ginsbourger, J. Baccou, and J. Liandrat. Warped Gaussian processes and derivative-based sequential designs for functions with heterogeneous variations. *SIAM/ASA Journal on Uncertainty Quantification*, 6(3):991–1018, 2018.

G. Matheron. The theory of regionalized variables and its applications. Technical Report Les cahiers du CMM de Fontainebleau, Fasc. 5, Ecole des Mines de Paris, 1971.

G. Matheron. Splines and kriging; their formal equivalence. 1981.

J. E. Matheson and R. L. Winkler. Scoring rules for continuous probability distributions. *Management science*, 22(10):1087–1096, 1976.
**Proof of Proposition 25.** The proof follows closely those of Theorem 4.2 from Narcowich et al. (2006), Theorem 4.2 from Karvonen et al. (2020), and Theorem 1 from Wynne et al. (2021). It is therefore not entirely repeated for brevity. Note that Theorem 4.1 from Arcangéli et al. (2007) must be used with $q = +\infty$ and $q = 2$ but requires $\mathbb{X}$ to be open. Under Assumption 9, the adaptation for $q = +\infty$ is given by (35) in the proof of Proposition 10. For $q = 2$, use the fact that a locally Lipschitz boundary has zero measure.

The proofs in the references give the result for $h_n$ lower than some $h_0$ (not dependent of $(x_n)_{n\ge1}$ and $f$). Removing this condition by possibly increasing the constant makes no difficulty.

The fact that $h_n \lesssim n^{-1/d}$ for a quasi-uniform sequence (see, e.g., Wendland, 2004, Proposition 14.1) gives the second statement. ■

**Proof of Proposition 27.** Assume there are observations $(x_n^B, z_n^B)$ within $B$ and let $s_{x_n^B, z_n^B}$ be the associated (interpolation) predictor. We have:

$$ \|f - s_{R,n}\|_{L^{\infty}(B)} \le \|f - s_{x_n^B, z_n^B}\|_{L^{\infty}(B)} + \|s_{x_n^B, z_n^B} - s_{R,n}\|_{L^{\infty}(B)}. $$

Use the notations of Lemma 33, let $q_{n,B}$ be the separation distance of points within $B$ (with the convention $q_{n,B} = \text{diam}(\mathbb{X})$ if $\mathbb{X}_{n,B}$ has less than two distinct elements), and let $\rho_{n,B} = h_{n,B}/q_{n,B}$. It holds that $f|_B \in W_2^{\beta+d/2}(B)$ since $B \subset \mathbb{X}_0$ and $f$ is an $R$-element of $W_2^{\beta+d/2}(\mathbb{X})$. It can be shown that the restriction of the kernel to $B$ also has regularity $\alpha$ using (15) and Lemma 31. Therefore, we can apply Proposition 25 to $f|_B$ to get:

$$ \|f - s_{x_n^B, z_n^B}\|_{L^{\infty}(B)} \lesssim h_{n,B}^{\beta}(1 + \rho_{n,B}^{\alpha-\beta}) \|f\|_{W_2^{\beta+d/2}(B)} \lesssim h_{n}^{\beta}(1 + \rho_{n}^{\alpha-\beta}) \|s_{R,X}^{(\beta)}\|_{W_2^{\beta+d/2}(X)}, $$

where we used Lemma 33, the inequality $q_n \le q_{n,B}$, the definition (15), and the fact that $s_{R,X}^{(\beta)}$ coincides with $f$ on $B$.

Consider the restrictions on $B$ and observe that $s_{R,n} \in \mathcal{H}(B)$ and that $s_{x_n^B, z_n^B}$ is the associated interpolation predictor. One can combine (8) and Proposition 10 to get

$$ \|s_{x_n^B, z_n^B} - s_{R,n}\|_{L^{\infty}(B)} \lesssim h_{n,B}^{\alpha} \|s_{R,n}\|_{\mathcal{H}(B)} \lesssim h_{n}^{\alpha} \|s_{R,n}\|_{\mathcal{H}(X)} $$

using Lemmata 31 and 33.

Observe that the proof of Lemma A.1. from Karvonen et al. (2020) is valid without restrictions on the domain $\mathbb{X}$. Therefore, there exists a constant $C_\beta > 0$ such that, for each $n$, there exists a function $u_n^{(\beta)} \in W_2^{\alpha+d/2}(\mathbb{X})$ coinciding with $s_{R,X}^{(\beta)}$ on the design $x_n$ and satisfying

$$ \|u_n^{(\beta)}\|_{W_2^{\alpha+d/2}(X)} \le C_\beta q_n^{\beta-\alpha} \|s_{R,X}^{(\beta)}\|_{W_2^{\beta+d/2}(X)}. $$

To conclude, observe that $u_n^{(\beta)} \in \mathcal{H}_{R,n}$ and thus

$$ \|s_{R,n}\|_{\mathcal{H}(X)} \le \|u_n^{(\beta)}\|_{\mathcal{H}(X)} \lesssim \|u_n^{(\beta)}\|_{W_2^{\alpha+d/2}(X)} \lesssim q_n^{\beta-\alpha} \|s_{R,X}^{(\beta)}\|_{W_2^{\beta+d/2}(X)}. \quad (38) $$

The previous display is also true when $B$ does not contain observations. Coordination with (6) and Lemma 33 shows that (29) also holds in this case.
M. D. McKay, R. J. Beckman, and W. J. Conover. 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, 2000.

E. Meeds and S. Osindero. An alternative infinite mixture of Gaussian process experts. In *Advances in Neural Information Processing Systems 18 (NIPS 2005)*, pages 883–890. MIT Press, 2006.

E. Merrill, A. Fern, X. Z. Fern, and N. Dolatnia. An empirical study of Bayesian optimization: Acquisition versus partition. *The Journal of Machine Learning Research*, 22:4–1, 2021.

J. Mockus, V. Tiesis, and A. Zilinskas. *The application of Bayesian methods for seeking the extremum*, volume 2, pages 117–129. 1978.

S. Nadarajah and S. Kotz. Exact distribution of the max/min of two Gaussian random variables. *IEEE Transactions on Very Large Scale Integration (VLSI) Systems*, 16(2):210–212, 2008.

J. Nocedal and S. J. Wright. Quadratic programming. *Numerical optimization*, pages 448–492, 2006.

D. J. Nott and William T. M. Dunsmuir. Estimation of nonstationary spatial covariance structure. *Biometrika*, 89(4):819–829, 2002.

C. Park and D. Apley. Patchwork kriging for large-scale Gaussian process regression. *The Journal of Machine Learning Research*, 19(1):269–311, 2018.

E. Parzen. Statistical inference on time series by Hilbert space methods, I. Technical report, STAN-FORD UNIV CA APPLIED MATHEMATICS AND STATISTICS LABS, 1959.

S. J. Petit, J. Bect, P. Feliot, and E. Vazquez. Model parameters in Gaussian process interpolation: an empirical study of selection criteria. 2021. URL https://arxiv.org/abs/2107.06006.

V. Picheny, S. Vakili, and A. Artemev. Ordinal Bayesian optimisation. *arXiv preprint arXiv:1912.02493*, 2019.

L. Pronzato and M. J. Rendas. Bayesian local kriging. *Technometrics*, 59(3):293–304, 2017.

C. E. Rasmussen and Z. Ghahramani. Infinite mixtures of Gaussian process experts. *Advances in neural information processing systems*, 2:881–888, 2002.

C. E. Rasmussen and C. K. I. Williams. *Gaussian Processes for Machine Learning*. Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, USA, 2006.

R. G. Regis. Constrained optimization by radial basis function interpolation for high-dimensional expensive black-box problems with infeasible initial points. *Engineering Optimization*, 46(2): 218–243, 2014.

J. Rivoirard and T. Romary. Continuity for kriging with moving neighborhood. *Mathematical Geosciences*, 43(4):469–481, 2011.
The fact that $h_n \lesssim n^{-1/d}$ for a quasi-uniform sequence (see, e.g., Wendland, 2004, Proposition 14.1) gives the second statement. ■

**Proof of Proposition 28.** Let $g \in \mathcal{H}(X)$ make $f$ an $R$-element of $\mathcal{H}(X)$. Applying Proposition 12 to $g$ proves one implication.

Conversely, the sequence $(\|s_{R,n}\|_{\mathcal{H}(X)})_{n \ge 1}$ is non-decreasing. Suppose that it is bounded. One can then extract a weakly convergent subsequence from $(s_{R,n})_{n \ge 1}$. The sets $\mathcal{H}_{R,n}$ are closed and convex and hence weakly closed. The weakly convergent subsequence is ultimately in each of these sets so its limit lies in their intersection. It can be checked that the existence of a continuous function coinciding with $f$ on $X_0$ and defining the same $X_j$s ensures that this intersection is equal to $\mathcal{H}_{R,X}$. ■

**Proof of Proposition 29.** The result is given by (38) with $\alpha = \nu$ and $q_n \lesssim n^{-1/d}$ by quasi-uniformity (see, e.g., Wendland, 2004, Proposition 14.1). ■

**Proof of Proposition 30.** The proof follows the ideas of Karvonen (2023, Theorem 3.11 and Theorem 3.12). Write $K_{n,\nu}$ to underline the dependence of the covariance matrix on the parameter $\nu$. Let $\varepsilon > 0$ be small enough and write

$$ \mathcal{M}(\nu) = \min_{z \in C_{R,n}} \mathcal{L}(\nu; z) \propto \text{constant} + \log(\det(K_{n,\nu})) + \min_{z \in C_{R,n}} z^T K_{n,\nu}^{-1} z. $$

For $\tau = \beta_R(f) + d/2 - \varepsilon/2$, we are going to prove that $\inf_{\nu_{\min} \le \nu \le \beta_R(f)+d/2-\varepsilon} \mathcal{M}(\nu) - \mathcal{M}(\tau)$ is ultimately positive.

First, observe that Proposition 3.7 of Karvonen (2023) applies in our setting (by potentially transforming $X$ to deal with an isotropic covariance function). Our Proposition 10 adapts Proposition 3.6 from the previous reference in the setting of Assumption 9. Proceed as in the reference by using these two propositions and Stirling's formula to get

$$ \inf_{\nu_{\min} \le \nu \le \beta_R(f) + d/2 - \varepsilon} \log(\det(K_{n,\nu})) - \log(\det(K_{n,\tau})) \ge \mathcal{O}(n) + \frac{\varepsilon}{d} \log(n!) \sim \frac{\varepsilon}{d} n \log(n). $$

To conclude, observe that $f$ is an $R$-element of $W_2^{\beta_R(f)-\varepsilon/2+d/2}(X)$ by definition, so (38) and quasi-uniformity yields:

$$ \min_{z \in C_{R,n}} z^T K_{n,\tau}^{-1} z \lesssim n. $$

■
I. Rychlik, P. Johannesson, and M. R. Leadbetter. Modelling and statistical analysis of ocean-wave data using transformed Gaussian processes. *Marine Structures*, 10(1):13–47, 1997.

E. Snelson, C. E. Rasmussen, and Z. Ghahramani. Warped Gaussian processes. *Advances in neural information processing systems*, 16:337–344, 2004.

N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. In *International Conference on Machine Learning*, pages 1015–1022, 2010.

M. L. Stein. Nonstationary spatial covariance functions. *Unpublished technical report*, 2005.

M. L. Stein. 2010 Rietz lecture: When does the screening effect hold? *The Annals of Statistics*, 39 (6):2795–2819, 2011.

M.L. Stein. *Interpolation of Spatial Data: Some Theory for kriging*. Springer Series in Statistics. Springer New York, 1999.

I. Steinwart, D. Hush, and C. Scovel. An explicit description of the reproducing kernel Hilbert spaces of Gaussian rbf kernels. *IEEE Transactions on Information Theory*, 52(10):4635–4643, 2006.

R. S. Strichartz. Multipliers on fractional Sobolev spaces. *Journal of Mathematics and Mechanics*, 16(9):1031–1060, 1967.

S. Surjanovic and D. Bingham. Virtual library of simulation experiments: Test functions and datasets. Retrieved October 13, 2020, from http://www.sfu.ca/~ssurjano/branin.html, 2013.

V. Tresp. Mixtures of Gaussian processes. *Advances in neural information processing systems*, pages 654–660, 2001.

E. Vazquez and J. Bect. Pointwise consistency of the kriging predictor with known mean and covariance functions. In *mODa 9–Advances in Model-Oriented Design and Analysis*, pages 221–228. Springer, 2010a.

E. Vazquez and J. Bect. Convergence properties of the expected improvement algorithm with fixed mean and covariance functions. *Journal of Statistical Planning and inference*, 140(11):3088–3095, 2010b.

E. Vazquez and J. Bect. A new integral loss function for Bayesian optimization. *arXiv preprint arXiv:1408.4622*, 2014.

J. Ver Hoef, N. Cressie, and R. Barry. Flexible spatial models for kriging and cokriging using moving averages and the fast fourier transform (fft). *Journal of Computational and Graphical Statistics - J COMPUT GRAPH STAT*, 13:265–282, 2004.
**Appendix D. EI Benchmark Results**

The results are provided in Figures 11–15, for all the test functions from Table 1, using the same format as in Figure 8.

Observe that all the EGO-R methods yield (sometimes very) substantial improvements on Beale, Goldstein-Price, Six-hump Camel, Dixon-Price (4), all instances of Perm and Rosenbrock, Three-hump Camel, and Zakharov (4) and (6). Focusing on these test functions, the concentration heuristic appears to give the best results, while the constant heuristic brings the least improvement. The spatial heuristic gives intermediate results between the two, sometimes giving results similar to the concentration heuristic and, at other times, results similar to the constant heuristic.

Taking variability (colored areas) into account, EGO-R with the constant and spatial heuristics performs only slightly better than EGO on Cross-in-Tray and Dixon-Price (6) and (10). The concentration heuristic, however, brings substantial improvements on these test functions, though. All EGO-R methods perform as EGO on Branin, Ackley (4) and (6), Log-Goldstein-Price, Michalewicz (4) and (10), all the Shekel instances, Hartman (3), and Zakharov (10).

The performance of the EGO-R methods is rarely worse than EGO, and never dramatically so. Nevertheless, the concentration heuristic yields substantially higher median empirical best minima on Ackley (10), and Michalewicz (6) and (10). This is also the case of the spatial heuristic on Ackley (10).

**Remark 39** The results are spread out for the Hartman (6) function. Investigations have shown this is due to a fraction of runs getting trapped in local minima. The proportions of runs reaching, for *n* = 300, an empirical best evaluation close to the global minimum are 72%, 68%, 67%, and 66% for EGO and EGO-R with the constant, spatial, and concentration heuristics, respectively. (To measure the closeness to the global minimum, a list $x_1^*, \dots, x_N^*$ of local minima was estimated by manually analyzing the output of repetitions of a random-start local optimizer. Then, a function value $f(x)$ is said to be close to the global minimum if $|f(x) - f(x_i^*)|$ is minimized by $\arg\min_i f(x_i^*)$.)
J. Villemonteix, E. Vazquez, and E. Walter. An informational approach to the global optimization of expensive-to-evaluate functions. *Journal of Global Optimization*, 44(4):509–534, 2009.

P. Virtanen, R. Gommers, T. E. Oliphant, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods*, 17:261–272, 2020.

G. Wahba. *Spline Models for Observational Data*. Society for Industrial and Applied Mathematics, Philadelphia, 1990.

H. Wendland. *Scattered data approximation*, volume 17. Cambridge university press, 2004.

Y. Xiang, D. Y. Sun, W. Fan, and X. G. Gong. Generalized simulated annealing algorithm and its application to the Thomson model. *Physics Letters A*, 233(3):216–220, 1997.

Y. Xiong, W. Chen, D. Apley, and X. Ding. A non-stationary covariance-based kriging method for metamodelling in engineering design. *International Journal for Numerical Methods in Engineering*, 71(6):733–756, 2007.

Y. Yang and J. Ma. An efficient em approach to parameter learning of the mixture of Gaussian processes. In *International Symposium on Neural Networks*, pages 165–174. Springer, 2011.

C. Yuan and C. Neubauer. Variational mixture of Gaussian process experts. In *Advances in Neural Information Processing Systems*, pages 1897–1904, 2009.

S. E. Yuksel, J. N. Wilson, and P. D. Gader. Twenty years of mixture of experts. *IEEE transactions on neural networks and learning systems*, 23(8):1177–1193, 2012.

H. Zhang and Y. Wang. Kriging and cross-validation for massive spatial data. *Environmetrics*, 21 (3-4):290–304, 2010.
<figure><img src="image_18.png" /><figcaption>Figure 11: EGO and EGO-R results for a subset of the test functions from Table 1. Same legend as in Figure 8.</figcaption></figure>
Cross-in-Tray

<figure><img src="image_19.png" /></figure>

Dixon-Price (4)

<figure><img src="image_20.png" /></figure>

Dixon-Price (6)

<figure><img src="image_21.png" /></figure>

Dixon-Price (10)

<figure><img src="image_22.png" /></figure>

Log-Goldstein-Price

<figure><img src="image_23.png" /></figure>

Goldstein-Price

<figure><img src="image_24.png" /><figcaption>Figure 12: EGO and EGO-R results for a subset of the test functions from Table 1. Same legend as in Figure 8.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Hartman (6)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Michalewicz (4)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Michalewicz (6)</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Michalewicz (10)</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Perm (4)</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure 13: EGO and EGO-R results for a subset of the test functions from Table 1. Same legend as in Figure 8.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure 14: EGO and EGO-R results for a subset of the test functions from Table 1. Same legend as in Figure 8.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Shekel (10)</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Three-hump Camel</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Zakharov (4)</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Zakharov (6)</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Zakharov (10)</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Figure 15: EGO and EGO-R results for a subset of the test functions from Table 1. Same legend as in Figure 8.</figcaption></figure>
Appendix E. UCB Benchmark Results

The Upper Confidence Bound (UCB, Cox and John, 1992; Srinivas et al., 2010) is an other sampling criterion for function minimization. Given a parameter $\beta > 0$, the function $f$ is minimized by choosing sequentially $x_{n+1}$ minimizing

$$ \mu_n(x) - \beta \sigma_n(x), \qquad (39) $$

with respect to $x \in \mathcal{X}$. In other words, the UCB algorithm minimizes an optimistic quantile of the Gaussian predictive distribution to select the next observation. Large values of $\beta$ favor exploration.

Varying $\beta$ slowly with $n$ helps to prove convergence (Srinivas et al., 2010). Nevertheless, a few experiments (not shown here) suggest that, overall, on the test functions from Table 1, the best vanilla GP baselines were obtained with a fixed value. Consequently, this finite-sample study was restricted to a fixed value of $\beta$ as in Cox and John (1992) and Picheny et al. (2013). The value was chosen so that (39) corresponds to the 10%-quantile of the predictive distribution.

As with the EI criterion, the application of reGP to the UCB algorithm simply involves plugging the reGP predictive distribution into the criterion (39). For choosing the validation threshold $t_n^{(0)}$, we will again compare the three heuristics introduced in Section 5.3.1, namely constant, concentration, and spatial. (Given $t_n^{(0)}$, the relaxation range candidates will also be defined as in Section 5.3.1.) More generally, the assessment was carried out using the methodology described in Section 5.3.1. The criterion (39) was optimized using a subset simulation algorithm. The results are shown in Figures 16–20.

When applied to UCB, all reGP variants outperform the standard UCB algorithm on Beale, Six-hump Camel, Dixon-Price (4) and (10), Goldstein-Price, all instances of Perm and Rosenbrock, Three-hump Camel, Zakharov (4) and (6). Furthermore, the concentration heuristic yields the best results on these test functions.

The spatial and constant heuristics perform as standard GPs on Branin, and Dixon-Price (6). Only the concentration heuristic brings a substantial benefit. However, GPs and all reGP variants give similar results on Ackley (4), Michalewicz (6) and (10), Hartman (3), and Zakharov (10). The concentration heuristic underperforms on Ackley (6) and (10).

The results are spread out for all algorithms on Cross-in-Tray, Log-Goldstein-Price, Hartman (6), Michalewicz (4), and Shekel (5), (7), and (10). As with the Hartman (6) function with EI, this is caused by local minima. This is a known feature of the standard UCB algorithm with fixed $\beta$ (see Jones, 2001, Section 5). (However, as mentioned earlier, this produced the best baselines, overall, in this study.) Table 2 shows the fractions of runs getting close to the global minima on these test functions. Observe that the reGP variants are not routinely more often trapped in local minima. Furthermore, the figures show that the concentration heuristic still benefits runs reaching the global minimum.
<table><caption>Table 2: Proportion of UCB runs reaching, for $n = 300$, an empirical best evaluation close to the global minimum. These proportions are estimated as described in Appendix D.</caption><thead><tr><th>Problem</th><th>GP</th><th>constant</th><th>spatial</th><th>concentration</th></tr></thead><tbody><tr><td>Hartman (6)</td><td>66%</td><td>65%</td><td>65%</td><td>67%</td></tr><tr><td>Shekel 5</td><td>31%</td><td>32%</td><td>36%</td><td>36%</td></tr><tr><td>Shekel 7</td><td>41%</td><td>44%</td><td>46%</td><td>47%</td></tr><tr><td>Shekel 10</td><td>19%</td><td>16%</td><td>16%</td><td>16%</td></tr><tr><td>Log-Goldstein-Price</td><td>87%</td><td>84%</td><td>81%</td><td>94%</td></tr><tr><td>Cross-in-Tray</td><td>97%</td><td>88%</td><td>86%</td><td>88%</td></tr><tr><td>Michalewicz (4)</td><td>26%</td><td>22%</td><td>16%</td><td>24%</td></tr></tbody></table>
<figure><img src="image_38.png" /><figcaption>Figure 16: UCB results for a subset of the test functions from Table 1. Same legend as in Figure 8. The red color stands for the UCB algorithm with vanilla GPs and the blue, green, and black colors for UCB with reGP, using the constant, concentration, and spatial heuristics respectively.</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>Figure 17: UCB results for a subset of the test functions from Table 1. Same legend as in Figure 16.</figcaption></figure>
Hartman (3)

<figure><img src="image_40.png" /></figure>

Hartman (6)

<figure><img src="image_41.png" /></figure>

Michalewicz (4)

<figure><img src="image_42.png" /></figure>

Michalewicz (6)

<figure><img src="image_43.png" /></figure>

Michalewicz (10)

<figure><img src="image_44.png" /></figure>

Perm (4)

<figure><img src="image_45.png" /><figcaption>Figure 18: UCB results for a subset of the test functions from Table 1. Same legend as in Figure 16.</figcaption></figure>
<figure><img src="image_46.png" /><figcaption>Figure 19: UCB results for a subset of the test functions from Table 1. Same legend as in Figure 16.</figcaption></figure>
<figure><img src="image_47.png" /><figcaption>Figure 20: UCB results for a subset of the test functions from Table 1. Same legend as in Figure 16.</figcaption></figure>
Appendix F. Application to the Estimation of an Excursion Set

Consider the problem of estimating the set {$x \in \mathbb{X}, f(x) \le u$}, for some $u \in \mathbb{R}$, as illustrated in Section 4.4. When $f$ is expensive to evaluate, a common approach is to use active learning with Gaussian process surrogates. A large variety of sampling criteria are available in the literature (see, e.g., Bryan et al., 2005; Bect et al., 2012; Chevalier et al., 2014; Bogunovic et al., 2016). To illustrate the benefits of reGP for the estimation of an excursion set, we choose the straddle heuristic of Bryan et al. (2005), which consists in choosing $x_{n+1}$ maximizing

$$1.96\sigma_n(x) - |\mu_n(x) - u|, \quad x \in \mathbb{X}. \qquad (40)$$

This simple (but effective) criterion rewards locations where the prediction is both close to $u$ and uncertain.

As with the EI and UCB criteria, the coordination of reGP with the straddle heuristic simply involves plugging the reGP predictive distribution into the criterion (40). Consequently, with observations $(\underline{x}_n, \underline{z}_n)$ and using the framework described by Algorithm 1, we are left to specify a range of interest $Q_n$ and a list $\mathbb{R} \setminus Q_n = R_n^{(0)} \supset \dots \supset R_n^{(G-1)}$ of relaxation range candidates. As in Section 4.4, symmetric ranges of interest $Q_n = (u - t_n^{(0)}, u + t_n^{(0)})$ are considered. For $t_n^{(0)}$, the three heuristics proposed in Section 5.3.1 are adapted by considering $\alpha$-quantiles of $|u - z_i|$. Then, given $t_n^{(0)}$, we propose to use $R_n^{(g)} = (-\infty, u - t_n^{(g)}] \cup [u + t_n^{(g)}, +\infty)$ with $t_n^{(g)}$ ranging logarithmically from $t_n^{(0)}$ to $\max_i |f(x_i) - u|$. As was done for optimization, we use the convention that $t_n^{(G-1)} > \max_i |f(x_i) - u|$ so that the G-th model is a vanilla GP.

As in Section 5.3.1, the experiments are conducted using Matern covariance functions with regularity $\nu = 5/2$ and constant mean functions. We also use $G = 10$ and initial designs of size $n_0 = 10d$. The criterion (40) is optimized using a subset simulation algorithm.

Four test cases are considered for the straddle heuristic in this work. First, the two functions $c_6$ and $c_6^7$ introduced in Section 4.4, with $u = 0$. Then, as in Section 5.3.1, we also consider the Goldstein-Price function, as well as its log-version. The vanilla Goldstein-Price function ranges from 3 to about $10^6$. To depart from optimization, the intermediate value $u = 1000$ is chosen, which is a spatial quantile with a level of about 25.4%. The corresponding value $u = \log(1000)$ is also used for the log-version. For each of these four test cases, we ran $n_{rep} = 100$ (random) repetitions, each with a budget of $n_{tot} = 300$ evaluations.

As often done for the estimation of excursion sets, the performance of the sequential strategies is measured using model predictions. More specifically, writing $P_n$ for a predictive distribution, we use a Sobol' sequence $\tilde{x}_1, \dots, \tilde{x}_N$ to compute the estimate

$$\frac{1}{N} \sum_{i=1}^{N} P_n(\xi(\tilde{x}_i) > u) 1_{f(\tilde{x}_i) \le u} + P_n(\xi(\tilde{x}_i) \le u) 1_{f(\tilde{x}_i) > u} \qquad (41)$$

of the (relative) expected volume of the symmetric difference between {$x \in \mathbb{X}, f(x) \le u$} and the random set {$x \in \mathbb{X}, \xi(x) \le u$}.

Figure 21 shows the evolution of quantiles of (41) across the $n_{rep}$ repetitions for the four test cases. Observe the large improvement brought by reGP on the $c_6$ test case. Using GPs, the evaluation metric (41) does not improve as $n$ increases, which is consistent with the observations of Feliot et al. (2017). However, the metric drops using reGP. Using the constant heuristic, the metric eventually stagnates at a value substantially lower than that obtained with vanilla GPs, whereas it continues
to decrease using the concentration and spatial heuristics. Figure 7 illustrates the relaxation obtained with $n = 300$ observations after a run of the concentration heuristic.

The performances of all the methods are comparable on $c_6^7$. The straddle heuristic with vanilla GPs is more effective than on $c_6$. All reGP instances are a little behind, with slightly higher medians but overlapping colored areas (except at the very beginning). Curiously, the concentration and spatial heuristics are less effective on $c_6^7$ than on $c_6$ (while the level sets are the same, since $t \mapsto t^7$ is increasing).

The three reGP instances give very similar results on the Goldstein-Price test case, with a substantial improvement over vanilla GPs. Finally, for the Log-Goldstein-Price test case, the reGP instances begin by giving very similar results to those obtained with GPs. Then, from around $n = 150$, the performances of the concentration and spatial heuristics improve to become substantially better. Observe that this is also the case for a fraction of runs with the constant heuristic.

In conclusion, we observe that reGP is a useful modeling technique for the estimation of excursion sets with the straddle heuristic. It can give clear improvements in the best cases, while giving results comparable to those obtained with GPs in the worst cases.
<figure><img src="image_48.png" /><figcaption>Figure 21: Evolution of quantiles of (41) against *n*. The solid lines represent the medians and, for each *n*, the shaded areas are eliminated by the 10% and 90% quantiles. The red color stands for the straddle heuristic with vanilla GPs and the blue, green, and black colors for straddle with reGP, using the constant, concentration, and spatial heuristics respectively.</figcaption></figure>
# Appendix G. Extension to the case of noisy observations

## G.1 Gaussian process regression and kernel ridge regression

We now consider the case where observations are corrupted by additive noise, a common setting when dealing with stochastic simulators (Baker et al., 2020) or when optimizing validation loss in machine learning tasks with mini-batch training (Snoek et al., 2012).

Let $x_n = (x_1, \dots, x_n) \in \mathbb{X}^n$ be the input locations and let $\underline{Z}_n = (Z_1, \dots, Z_n)^\mathrm{T} \in \mathbb{R}^n$ denote the noisy observations, modeled as

$$Z_i = f(x_i) + \varepsilon_i,$$

where $f: \mathbb{X} \to \mathbb{R}$ is the latent function of interest and $\varepsilon_i \sim \mathcal{N}(0, \eta^2)$ are independent Gaussian noise terms with common variance $\eta^2 > 0$.

We place a Gaussian process prior $\xi \sim \mathrm{GP}(\mu, k)$ on $f$, using the same notations as in Section 2. The posterior distribution of $\xi$ given the observations $\underline{Z}_n$ remains Gaussian:

$$\xi | \underline{Z}_n \sim \mathrm{GP}(\mu_n, k_n), \qquad (42)$$

where the posterior mean and covariance functions are

$$\mu_n(x) = \mu(x) + k(x, \underline{x}_n)(K_n + \eta^2 I)^{-1}(\underline{Z}_n - \mu(\underline{x}_n)), \\ k_n(x, y) = k(x, y) - k(x, \underline{x}_n)(K_n + \eta^2 I)^{-1}k(y, \underline{x}_n)^\mathrm{T}.$$

Here, $K_n$ is the $n \times n$ covariance matrix with entries $[K_n]_{ij} = k(x_i, x_j)$, and $I$ is the $n \times n$ identity matrix.

The posterior mean in the zero-mean case ($\mu = 0$) coincides with the solution of a kernel ridge regression (KRR) problem. This correspondence is a direct consequence of the representer theorem, as formalized below.

**Proposition 40 (Kimeldorf and Wahba, 1970)** *The solution to the regularized least-squares problem*

$$\left\{ \begin{array}{ll} \text{minimize} & \|h\|_{\mathcal{H}(\mathbb{X})}^2 + \eta^{-2} \sum_{i=1}^{n} (Z_i - h(x_i))^2 \\ \text{subject to} & h \in \mathcal{H}(\mathbb{X}) \end{array} \right.$$

is unique and given by

$$s_{\underline{Z}_n}^{(\eta)} = k(\cdot, \underline{x}_n)(K_n + \eta^2 I)^{-1} \underline{Z}_n.$$

Hence, when $\mu = 0$, the posterior mean $\mu_n$ coincides with the KRR solution $s_{\underline{Z}_n}^{(\eta)}$.

## G.2 Relaxed Gaussian process regression

As seen in Proposition 40, the posterior mean of a Gaussian process model does not interpolate noisy data. How can the goal-oriented modeling framework introduced in Section 3 be extended to this setting?

A natural extension is to adapt Definition 6 by optimizing the likelihood with respect to noisy observations falling in the relaxation range.

Assume that both the mean and the covariance functions belong to parametric families, indexed by $\theta \in \Theta$. Fix the noise variance $\eta^2 > 0$ (see Appendix G.3 for an estimator). Let $R = \cup_{j=1}^J R_j$
be a union of disjoint closed intervals and define the corresponding constraint set $C_{R,n} \subset \mathbb{R}^n$ as in Section 3.1.

We define the relaxed GP predictive distribution as the posterior distribution (42) of the GP conditioned on $\underline{Z}_n = \underline{z}_n^*$ and on the estimated parameter $\hat{\theta}_n$, where

$$ (\hat{\theta}_n, \underline{z}_n^*) = \underset{\theta \in \Theta, \underline{z} \in C_{R,n}}{\operatorname{arg min}} \mathcal{L}(\theta; \underline{z}), $$

and $\mathcal{L}$ is the (negative) log-likelihood of the observations written as

$$ \mathcal{L}(\theta; \underline{z}) \propto \log \det(K_n + \eta^2 I) + (\underline{z} - \mu(\underline{x}_n))^\top (K_n + \eta^2 I)^{-1} (\underline{z} - \mu(\underline{x}_n)) + \text{constant.} $$

Note that we recover Definition 6 when $\eta^2 = 0$.

When the mean function is zero and the covariance function $k$ is fixed, the resulting procedure can be formulated as a regularized regression problem:

**Proposition 41** Let $\mathcal{H}(X)$ be the RKHS associated with $k$, and let $\eta^2 > 0$ be fixed. Consider the relaxed regularized least-squares problem

$$ \left\{ \begin{array}{ll} \underset{h \in \mathcal{H}(X)}{\operatorname{minimize}} & \|h\|^2_{\mathcal{H}(X)} + \eta^{-2} \sum_{i=1}^{n} e(Z_i, h(x_i); R), \\ \text{subject to} & h \in \mathcal{H}(X), \end{array} \right. \qquad (43) $$

where the relaxed squared error $e(y_1, y_2; R)$ is defined for all $y_1, y_2 \in \mathbb{R}$ by

$$ e(y_1, y_2; R) = \begin{cases} \inf_{r \in R_j} (y_2 - r)^2 & \text{if } y_1 \in R_j \text{ for some } j, \\ (y_2 - y_1)^2 & \text{otherwise.} \end{cases} $$

Then, the unique solution of (43) is $s_{\underline{z}_n}^{(\eta)}$, where $\underline{z}_n^*$ is the unique solution of

$$ \underset{\underline{z} \in C_{R,n}}{\operatorname{arg min}} \underline{z}^\top (K_n + \eta^2 I)^{-1} \underline{z}. \qquad (44) $$

**Proof** Write $\Delta = \eta^2 I$. For $\underline{z} = (z_1, \dots, z_n)^\top \in \mathbb{R}^n$, straightforward calculations (using the identity $I - K_n(K_n + \Delta)^{-1} = \Delta(K_n + \Delta)^{-1}$) show that

$$ \|s_{\underline{z}}^{(\eta)}\|^2_{\mathcal{H}(X)} + \eta^{-2} \sum_{i=1}^{n} (z_i - s_{\underline{z}}^{(\eta)}(x_i))^2 = \underline{z}^\top (K_n + \Delta)^{-1} \underline{z}. $$

The quadratic problem (44) has a unique solution $\underline{z}_n^*$ since the matrix is positive-definite and the feasible set is closed and convex. Write $L^*$ for the corresponding minimum value.

For $h \in \mathcal{H}(X)$, the definition of $e$ yields

$$ \|h\|^2_{\mathcal{H}(X)} + \eta^{-2} \sum_{i=1}^{n} e(Z_i, h(x_i); R) = \inf_{\underline{z} = (z_1, \dots, z_n)^\top \in C_{R,n}} \|h\|^2_{\mathcal{H}(X)} + \eta^{-2} \sum_{i=1}^{n} (z_i - h(x_i))^2. $$
The infimum is reached by a $\underline{z}(h) = (z_1(h), \dots, z_n(h))^T \in C_{R,n}$. We then have

$$
\begin{align*}
\|h\|_{\mathcal{H}(\mathbf{X})}^2 + \eta^{-2} \sum_{i=1}^{n} e(Z_i, h(x_i); R) &= \|h\|_{\mathcal{H}(\mathbf{X})}^2 + \eta^{-2} \sum_{i=1}^{n} (z_i(h) - h(x_i))^2 \\
&\geq \|s_{\underline{z}(h)}^{(\eta)}\|_{\mathcal{H}(\mathbf{X})}^2 + \eta^{-2} \sum_{i=1}^{n} (z_i(h) - s_{\underline{z}(h)}^{(\eta)}(x_i))^2 \\
&= \underline{z}(h)^T (K_n + \Delta)^{-1} \underline{z}(h) \\
&\geq L^*.
\end{align*}
$$

If the previous inequalities are equalities, then $h = s_{\underline{z}_n}^{(\eta)}$ because of the unicity of the solution of (44) and by Proposition 40.

Finally, writing $\underline{z}_n^* = (z_1^*, \dots, z_n^*)$, we have

$$
\|s_{\underline{z}_n}^{(\eta)}\|_{\mathcal{H}(\mathbf{X})}^2 + \eta^{-2} \sum_{i=1}^{n} e(Z_i, s_{\underline{z}_n}^{(\eta)}(x_i); R) \le \|s_{\underline{z}_n}^{(\eta)}\|_{\mathcal{H}(\mathbf{X})}^2 + \eta^{-2} \sum_{i=1}^{n} (z_i^* - s_{\underline{z}_n}^{(\eta)}(x_i))^2 = L^*.
$$

<figure><img src="image_49.png" /></figure>

This shows that, in the fixed-kernel case, the relaxed GP regression procedure corresponds to a modified KRR problem: observations $Z_i$ outside the relaxation range $R$ are penalized by squared prediction error, while observations inside $R_j$ are only penalized by the squared distance from the prediction $h(x_i)$ to the interval $R_j$.

The resulting predictive distribution thus reflects a form of soft constraint on the latent function: fidelity to the data is enforced outside $R$, while inside $R$ the model is only constrained to remain consistent with the relaxation range.

G.3 Application to Bayesian optimization with the UCB criterion

We now illustrate the use of reGP in Bayesian optimization with noisy observations. We adopt the same UCB sampling strategy as described in Appendix E, using a fixed $\beta$ parameter. Importantly, we do not inflate the posterior variance with the estimated noise variance $\eta^2$; the UCB criterion is applied to the predictive distribution of the latent process $\xi$.

Two test functions are considered: the Goldstein-Price function and its logarithmic transformation. In both cases, Gaussian noise is added to the evaluations. For the Goldstein-Price function, we use $\eta = 100$ (standard deviation), which is large relative to the local variations of the function near its minimum. Specifically, while the function ranges from 3 to approximately $10^6$, the 25% and 5% spatial quantiles are around 1000 and 100, respectively. The noise level is thus significant near the global minimum, whereas extreme values act as outliers. This setting is favorable to the reGP strategy, which focuses on accurate modeling in target regions. For the Log-Goldstein-Price function (value range: 1.10 to 13.83), we set $\eta = 3$.

In both cases, the noise variance $\eta^2$ is assumed to be unknown and must be estimated from the data. Jointly estimating $\eta^2$, the model parameters, and the relaxed values by maximizing the likelihood is not advisable in the reGP framework, as the relaxed values tend to shrink toward the posterior mean function, potentially leading to a biased underestimation of the noise level. Instead,
we estimate $\eta^2$ using standard GP maximum likelihood, restricted to observations within the range of interest, and treat the result as a fixed value in the relaxed procedure.

For reGP, the relaxation threshold is selected from a predefined list using the LOO-tCRPS criterion, as described in Section 5.3.1. The spatial heuristic is adapted by setting the validation threshold $t_n^{(0)}$ to a quantile of the previous model's predictions at the previous design $\bar{x}_{n-1}$. Further adjustments are likely needed for the concentration heuristic.

We run $n_{\text{rep}} = 100$ independent repetitions with random initial designs of size $n_0 = 10d$ and total budgets of $n_{\text{tot}} = 300$ evaluations. Results are reported in Figure 22. We define $x_n^* = \arg\min \mu_n$ as the predicted minimizer at step $n$, and use $f(x_n^*)$ as an oracle benchmark.

On the Goldstein-Price function, reGP yields substantially lower median noisy evaluations up to approximately $n = 150$. This is confirmed by the superior oracle values $f(x_n^*)$ obtained with reGP during the same phase. After that point, the difference between noisy values and the global minimum becomes comparable to the noise level, and the two approaches yield similar performance. Notably, the estimation error $|f(x_n^*) - \mu_n(x_n^*)|$ remains consistently smaller with reGP across all steps.

On the Log-Goldstein-Price function, both methods perform similarly overall. However, standard GP models show higher 90% quantiles for $f(x_n^*)$ and $f(x_n) + \varepsilon_n$, indicating more variability in the predicted minima.

These results illustrate the robustness of reGP in noisy optimization tasks.
<figure><img src="image_50.png" /><figcaption>Figure 22: Bayesian optimization with noisy evaluations: performance of the UCB strategy. Left: Goldstein-Price function. Right: Log-Goldstein-Price function. Each plot shows medians (solid lines) and 10%/90% quantiles (shaded areas) as functions of the number of evaluations *n*. Top: noisy evaluations $f(x_n) + \epsilon_n$. Middle: oracle values $f(x_n^*)$ at predicted minimizers $x_n^* = \arg\min \mu_n$. Bottom: absolute estimation errors $|f(x_n^*) - \mu_n(x_n^*)|$.</figcaption></figure>
References

R. A. Adams and J. J. F. Fournier. Sobolev spaces. Elsevier, 2003.

R. Arcangéli, M. C. López de Silanes, and J. J. Torrens. An extension of a bound for functions in Sobolev spaces, with applications to (m, s)-spline interpolation and smoothing. *Numerische Mathematik*, 107(2):181–211, 2007.

N. Aronszajn. Theory of reproducing kernels. *Transactions of the American mathematical society*, 68(3):337–404, 1950.

S.-K. Au and J. L. Beck. Estimation of small failure probabilities in high dimensions by subset simulation. *Probabilistic engineering mechanics*, 16(4):263–277, 2001.

S. Ba and V. R. Joseph. Composite Gaussian process models for emulating expensive functions. *The Annals of Applied Statistics*, 6(4):1838–1860, 2012.

F. Bachoc and A. Lagnoux. Posterior contraction rates for constrained deep Gaussian processes in density estimation and classification. *Communications in Statistics-Theory and Methods*, 54(3): 774–811, 2025.

E. Baker, P. Barbillon, A. Fadikar, R. B. Gramacy, R. Herbei, D. Higdon, J. Huang, L. R. Johnson, A. Mondal, B. Pires, J. Sacks, and V. Sokolov. Stochastic simulators: An overview with opportunities. *arXiv preprint arXiv:2002.01321*, 2020.

J. Y. Bao, F. Ye, and Y. Yang. Screening effect in isotropic Gaussian processes. *Acta Mathematica Sinica*, 36(5), 2020.

J. Bect, D. Ginsbourger, L. Li, V. Picheny, and E. Vazquez. Sequential design of computer experiments for the estimation of a probability of failure. *Statistics and Computing*, 22:773–793, 2012.

R. Benassi, J. Bect, and E. Vazquez. Bayesian optimization using sequential Monte Carlo. In *Learning and Intelligent Optimization – 6th International Conference, LION 6*, pages 339–342, Paris, France, January 16-20 2012. Springer.

A. Benavoli, D. Azzimonti, and D. Piga. A unified framework for closed-form nonparametric regression, classification, preference and mixed problems with Skew Gaussian Processes. *Machine Learning*, 110(11):3095–3133, 2021.

A. Berlinet and C. Thomas-Agnan. *Reproducing kernel Hilbert spaces in probability and statistics*. Springer Science & Business Media, 2004.

M. J. Best and B. Ding. On the continuity of the minimum in parametric quadratic programs. *Journal of Optimization Theory and Applications*, 86(1):245–250, 1995.

E. Bodin, M. Kaiser, I. Kazlauskaitė, Z. Dai, N. Campbell, and C. H. Ek. Modulating surrogates for Bayesian optimization. In *International Conference on Machine Learning*, pages 970–979. PMLR, 2020.
I. Bogunovic, J. Scarlett, A. Krause, and V. Cevher. Truncated variance reduction: A unified approach to Bayesian optimization and level-set estimation. In *Advances in Neural Information Processing Systems*, volume 29, 2016.

Z. I. Botev. The normal law under linear restrictions: simulation and estimation via minimax tilting. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 79(1):125–148, 2017.

H. Brezis. *Functional Analysis, Sobolev Spaces and Partial Differential Equations*, volume 2. Springer, 2011.

B. Bryan, R. C. Nichol, C. R. Genovese, J. Schneider, C. J. Miller, and L. Wasserman. Active learning for identifying function threshold boundaries. In *Advances in Neural Information Processing Systems*, volume 18. MIT Press, 2005.

A. D. Bull. Convergence rates of efficient global optimization algorithms. *Journal of Machine Learning Research*, 12(10), 2011.

C. Chevalier and D. Ginsbourger. Fast computation of the multi-points expected improvement with applications in batch selection. In *International Conference on Learning and Intelligent Optimization*, pages 59–69. Springer, 2013.

C. Chevalier, J. Bect, D. Ginsbourger, E. Vazquez, V. Picheny, and Y. Richet. Fast parallel kriging-based stepwise uncertainty reduction with application to the identification of an excursion set. *Technometrics*, 56(4):455–465, 2014.

W. Chu and Z. Ghahramani. Gaussian processes for ordinal regression. *Journal of machine learning research*, 6(7), 2005.

D. D. Cox and S. John. A statistical method for global optimization. In *IEEE International Conference on Systems, Man, and Cybernetics*, pages 1241–1246. IEEE, 1992.

S. Da Veiga and A. Marrel. Gaussian process modeling with inequality constraints. In *Annales de la Faculté des sciences de Toulouse: Mathématiques*, volume 21, pages 529–555, 2012.

A. Damianou and N. D. Lawrence. Deep Gaussian processes. In *Artificial intelligence and statistics*, pages 207–215. PMLR, 2013.

E. Di Nezza, G. Palatucci, and E. Valdinoci. Hitchhiker's guide to the fractional Sobolev spaces. *Bulletin des Sciences Mathématiques*, 136(5):521–573, 2012.

O. Dubrule. Cross validation of kriging in a unique neighborhood. *Journal of the International Association for Mathematical Geology*, 15:687–699, 1983.

M. M. Dunlop, M. A. Girolami, A. M. Stuart, and A. L. Teckentrup. How deep are deep Gaussian processes? *Journal of Machine Learning Research*, 19(54):1–46, 2018.

D. Eriksson, M. Pearce, J. R. Gardner, R. C. Turner, and M. Poloczek. Scalable global optimization via local Bayesian optimization. In *Advances in Neural Information Processing Systems*, volume 32, 2019.
P. Feliot, J. Bect, and E. Vazquez. A Bayesian approach to constrained single-and multi-objective optimization. *Journal of Global Optimization*, 67(1-2):97-133, 2017.

P. I. Frazier, W. B. Powell, and S. Dayanik. A knowledge-gradient policy for sequential information collection. *SIAM Journal on Control and Optimization*, 47(5):2410-2439, 2008.

A. Gessner, O. Kanjilal, and P. Hennig. Integrals over Gaussians under linear domain constraints. In *International Conference on Artificial Intelligence and Statistics*, pages 2764-2774. PMLR, 2020.

M. N. Gibbs. *Bayesian Gaussian processes for regression and classification*. PhD thesis, Citeseer, 1998.

T. G. Gneiting and A. E. Raftery. Strictly proper scoring rules, prediction, and estimation. *Journal of the American Statistical Association*, 102:359-378, 2007.

T. G. Gneiting and R. Ranjan. Comparing density forecasts using threshold-and quantile-weighted scoring rules. *Journal of Business & Economic Statistics*, 29(3):411-422, 2011.

T. G. Gneiting, A. E. Raftery, A. H. Westveld III, and T. Goldman. Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation. *Monthly Weather Review*, 133(5):1098-1118, 2005.

R. B. Gramacy and H. K. H. Lee. Bayesian treed Gaussian process models with an application to computer modeling. *Journal of the American Statistical Association*, 103(483):1119-1130, 2008.

P. Grisvard. *Elliptic Problems in Nonsmooth Domains*. Pitman, 1985.

A. Hebbal, L. Brevault, and N. Melab. Bayesian optimization using deep Gaussian processes with applications to aerospace system design. *Optimization and Engineering*, 22(1):321-361, 2021.

D. Higdon. A process-convolution approach to modelling temperatures in the north atlantic ocean. *Environmental and Ecological Statistics*, 5:173-190, 1998.

D. Higdon. Space and space-time modeling using process convolutions. *Quantitative Methods for Current Environmental Issues*, 04 2002.

A. Iske. *Approximation Theory and Algorithms for Data Analysis*. Springer, 2018.

K. Jakkala. Deep Gaussian processes: A survey. arXiv preprint arXiv:2106.12135, 2021.

D. Jones. A taxonomy of global optimization methods based on response surfaces. *Journal of Global Optimization*, 21:345-383, 2001.

D. Jones, M. Schonlau, and W. Welch. Efficient global optimization of expensive black-box functions. *Journal of Global Optimization*, 13:455-492, 12 1998.

T. Karvonen. Error bounds and the asymptotic setting in kernel-based approximation. *Dolomites Research Notes on Approximation*, 15(3):65-77, 2022.

T. Karvonen. Asymptotic bounds for smoothness parameter estimates in Gaussian process interpolation. *SIAM/ASA Journal on Uncertainty Quantification*, 11(4):1225-1257, 2023.
T. Karvonen, G. Wynne, F. Tronarp, C. Oates, and S. Sarkka. Maximum likelihood estimation and uncertainty quantification for Gaussian process approximation of deterministic functions. *SIAM/ASA Journal on Uncertainty Quantification*, 8(3):926-958, 2020.

G. S. Kimeldorf and G. Wahba. A correspondence between Bayesian estimation on stochastic processes and smoothing by splines. *The Annals of Mathematical Statistics*, 41(2):495–502, 1970.

S. Lacoste-Julien, F. Huszár, and Z. Ghahramani. Approximate inference for the loss-calibrated Bayesian. In *Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics*, volume 15, pages 416–424. PMLR, 11–13 Apr 2011.

M. Lázaro-Gredilla. Bayesian warped Gaussian processes. *Advances in Neural Information Processing Systems*, 25:1619–1627, 2012.

S. Lerch and T. L. Thorarinsdottir. Comparison of non-homogeneous regression models for probabilistic wind speed forecasting. *Tellus A: Dynamic Meteorology and Oceanography*, 65(1):21206, 2013.

J. R. Lewis, S. N. MacEachern, and Y. Lee. Bayesian restricted likelihood methods: Conditioning on insufficient statistics in Bayesian regression. *Bayesian Analysis*, 1(1), 2021.

A. F. López-Lopera, F. Bachoc, N. Durrande, and O. Roustant. Finite-dimensional Gaussian approximation with linear inequality constraints. *SIAM/ASA Journal on Uncertainty Quantification*, 6(3):1224–1255, 2018.

M. Lukić and J. Beder. Stochastic processes with sample paths in reproducing kernel Hilbert spaces. *Transactions of the American Mathematical Society*, 353(10):3945–3969, 2001.

H. Maatouk and X. Bay. Gaussian process emulators for computer experiments with inequality constraints. *Mathematical Geosciences*, 49(5):557–582, 2017.

S. Marmin, D. Ginsbourger, J. Baccou, and J. Liandrat. Warped Gaussian processes and derivative-based sequential designs for functions with heterogeneous variations. *SIAM/ASA Journal on Uncertainty Quantification*, 6(3):991–1018, 2018.

G. Matheron. The theory of regionalized variables and its applications. Technical Report Les cahiers du CMM de Fontainebleau, Fasc. 5, Ecole des Mines de Paris, 1971.

J. E. Matheson and R. L. Winkler. Scoring rules for continuous probability distributions. *Management Science*, 22(10):1087–1096, 1976.

M. McIntire, D. Ratner, and S. Ermon. Sparse Gaussian processes for Bayesian optimization. In *UAI*, 2016.

M. D. McKay, R. J. Beckman, and W. J. Conover. 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, 2000.

E. Meeds and S. Osindero. An alternative infinite mixture of Gaussian process experts. In *Advances in Neural Information Processing Systems 18 (NIPS 2005)*, pages 883–890. MIT Press, 2006.
E. Merrill, A. Fern, X. Z. Fern, and N. Dolatnia. An empirical study of Bayesian optimization: Acquisition versus partition. *The Journal of Machine Learning Research*, 22:4–1, 2021.

J. Mockus, V. Tiesis, and A. Zilinskas. The application of Bayesian methods for seeking the extremum. *Towards Global Optimization*, 2:117–129, 1978.

M. J. Morais and J. W. Pillow. Loss-calibrated expectation propagation for approximate Bayesian decision-making. *arXiv preprint arXiv:2201.03128*, 2022.

H. B. Moss, S. W. Ober, and V. Picheny. Inducing point allocation for sparse Gaussian processes in high-throughput Bayesian optimisation. In *International Conference on Artificial Intelligence and Statistics*, pages 5213–5230. PMLR, 2023.

S. Nadarajah and S. Kotz. Exact distribution of the max/min of two Gaussian random variables. *IEEE Transactions on Very Large Scale Integration (VLSI) Systems*, 16(2):210–212, 2008.

F. J. Narcowich, J. D. Ward, and H. Wendland. Sobolev error estimates and a Bernstein inequality for scattered data interpolation via radial basis functions. *Constructive Approximation*, 24(2): 175–186, 2006.

J. Nocedal and S. J. Wright. Quadratic programming. *Numerical Optimization*, pages 448–492, 2006.

D. J. Nott and William T. M. Dunsmuir. Estimation of nonstationary spatial covariance structure. *Biometrika*, 89(4):819–829, 2002.

C. Park and D. Apley. Patchwork kriging for large-scale Gaussian process regression. *The Journal of Machine Learning Research*, 19(1):269–311, 2018.

S. J. Petit. An asymptotic study of the joint maximum likelihood estimation of the regularity and the amplitude parameters of a periodized Matérn model. *Electronic Journal of Statistics*, 19(1): 2052–2094, 2025.

S. J. Petit, J. Bect, P. Feliot, and E. Vazquez. Parameter selection in Gaussian process interpolation: an empirical study of selection criteria. *SIAM/ASA Journal on Uncertainty Quantification*, 11(4): 1308–1328, 2023.

V. Picheny, T. Wagner, and D. Ginsbourger. A benchmark of kriging-based infill criteria for noisy optimization. *Structural and Multidisciplinary Optimization*, 48:607–626, 2013.

V. Picheny, S. Vakili, and A. Artemev. Ordinal Bayesian optimisation. *arXiv preprint arXiv:1912.02493*, 2019.

L. Pronzato and M. J. Rendas. Bayesian local kriging. *Technometrics*, 59(3):293–304, 2017.

C. E. Rasmussen and Z. Ghahramani. Infinite mixtures of Gaussian process experts. *Advances in Neural Information Processing Systems*, 2:881–888, 2002.

C. E. Rasmussen and C. K. I. Williams. *Gaussian Processes for Machine Learning*. Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, USA, 2006.
R. G. Regis. Constrained optimization by radial basis function interpolation for high-dimensional
expensive black-box problems with infeasible initial points. *Engineering Optimization*, 46(2):
218-243, 2014.

J. Rivoirard and T. Romary. Continuity for kriging with moving neighborhood. *Mathematical Geosciences*, 43(4):469–481, 2011.

I. Rychlik, P. Johannesson, and M. R. Leadbetter. Modelling and statistical analysis of ocean-wave
data using transformed Gaussian processes. *Marine Structures*, 10(1):13–47, 1997.

M. Scheurer, R. Schaback, and M. Schlather. Interpolation of spatial data – a stochastic or a deterministic problem? *European Journal of Applied Mathematics*, 24(4):601–629, 2013.

K. Schittkowski. The nonlinear programming method of Wilson, Han, and Powell with an augmented Lagrangian type line search function: Part 1: Convergence analysis. *Numerische Mathematik*, 38:83–114, 1982.

B. Schölkopf, R. Herbrich, and A. J. Smola. A generalized representer theorem. In *International Conference on Computational Learning Theory*, pages 416–426. Springer, 2001.

E. Snelson, C. E. Rasmussen, and Z. Ghahramani. Warped Gaussian processes. *Advances in Neural Information Processing Systems*, 16:337–344, 2004.

J. Snoek, H. Larochelle, and R. P. Adams. Practical Bayesian optimization of machine learning
algorithms. In *Advances in Neural Information Processing Systems*, volume 25, 2012.

N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger. Gaussian process optimization in the bandit
setting: No regret and experimental design. In *International Conference on Machine Learning*,
pages 1015–1022, 2010.

M. L. Stein. Nonstationary spatial covariance functions. *Unpublished technical report*, 2005.

M. L. Stein. 2010 Rietz lecture: When does the screening effect hold? *The Annals of Statistics*, 39
(6):2795–2819, 2011.

M.L. Stein. *Interpolation of Spatial Data: Some Theory for kriging*. Springer Series in Statistics.
Springer New York, 1999.

I. Steinwart, D. Hush, and C. Scovel. An explicit description of the reproducing kernel Hilbert
spaces of Gaussian RBF kernels. *IEEE Transactions on Information Theory*, 52(10):4635–4643,
2006.

R. S. Strichartz. Multipliers on fractional Sobolev spaces. *Journal of Mathematics and Mechanics*, 16(9):1031–1060, 1967.

S. Surjanovic and D. Bingham. Virtual library of simulation experiments: Test functions and
datasets. Retrieved October 13, 2020, from http://www.sfu.ca/~ssurjano/branin.html,
2013.
S. Takeno, M. Nomura, and M. Karasuyama. Towards practical preferential Bayesian optimization with skew Gaussian processes. In *International Conference on Machine Learning*, pages 33516–33533. PMLR, 2023.

V. Tresp. Mixtures of Gaussian processes. In *Advances in Neural Information Processing Systems*, pages 654–660, 2001.

E. Vazquez and J. Bect. Pointwise consistency of the kriging predictor with known mean and covariance functions. In *mODa 9–Advances in Model-Oriented Design and Analysis*, pages 221–228. Springer, 2010a.

E. Vazquez and J. Bect. Convergence properties of the expected improvement algorithm with fixed mean and covariance functions. *Journal of Statistical Planning and Inference*, 140(11):3088–3095, 2010b.

E. Vazquez and J. Bect. A new integral loss function for Bayesian optimization. *arXiv preprint arXiv:1408.4622*, 2014.

J. Ver Hoef, N. Cressie, and R. Barry. Flexible spatial models for kriging and cokriging using moving averages and the fast Fourier transform (FFT). *Journal of Computational and Graphical Statistics*, 13:265–282, 2004.

J. Villemonteix, E. Vazquez, and E. Walter. An informational approach to the global optimization of expensive-to-evaluate functions. *Journal of Global Optimization*, 44(4):509–534, 2009.

G. Wahba. *Spline Models for Observational Data*. Society for Industrial and Applied Mathematics, Philadelphia, 1990.

W. Wang and B.-Y. Jing. Gaussian process regression: Optimality, robustness, and relationship with kernel ridge regression. *Journal of Machine Learning Research*, 23(193):1–67, 2022.

H. Wendland. *Scattered Data Approximation*, volume 17. Cambridge university press, 2004.

G. Wynne, F. X. Briol, and M. Girolami. Convergence guarantees for Gaussian process means with misspecified likelihoods and smoothness. *Journal of Machine Learning Research*, 22(123):1–40, 2021.

Y. Xiong, W. Chen, D. Apley, and X. Ding. A non-stationary covariance-based kriging method for metamodelling in engineering design. *International Journal for Numerical Methods in Engineering*, 71(6):733–756, 2007.

A. Yang, C. Li, S. Rana, S. Gupta, and S. Venkatesh. Sparse spectrum Gaussian process for Bayesian optimization. In *Advances in Knowledge Discovery and Data Mining*, pages 257–268. Springer, 2021.

Y. Yang and J. Ma. An efficient EM approach to parameter learning of the mixture of Gaussian processes. In *International Symposium on Neural Networks*, pages 165–174. Springer, 2011.

C. Yuan and C. Neubauer. Variational mixture of Gaussian process experts. In *Advances in Neural Information Processing Systems*, pages 1897–1904, 2009.
S. E. Yuksel, J. N. Wilson, and P. D. Gader. Twenty years of mixture of experts. *IEEE Transactions on Neural Networks and Learning Systems*, 23(8):1177–1193, 2012.

H. Zhang and Y. Wang. Kriging and cross-validation for massive spatial data. *Environmetrics*, 21 (3–4):290–304, 2010.