Efficient Robust Bayesian Optimization for Arbitrary Uncertain inputs
Authors: Lin Yang
Junlong Lyu
Wenlong Lyu
Zhitang Chen
Venue: NeurIPS 2023 poster (2023)
Keywords: bayesian optimization
robust optimization
TL;DR: A robust Bayesian optimization method with MMD-based kernel design andn Nystrom-approximation-based acceleration.

Abstract: Bayesian Optimization (BO) is a sample-efficient optimization algorithm widely employed across various applications. In some challenging BO tasks, input uncertainty arises due to the inevitable randomness in the optimization process, such as machining errors, execution noise, or contextual variability. This uncertainty deviates the input from the intended value before evaluation, resulting in significant performance fluctuations in the final result. In this paper, we introduce a novel robust Bayesian Optimization algorithm, AIRBO, which can effectively identify a robust optimum that performs consistently well under arbitrary input uncertainty. Our method directly models the uncertain inputs of arbitrary distributions by empowering the Gaussian Process with the Maximum Mean Discrepancy (MMD) and further accelerates the posterior inference via Nystrom approximation. Rigorous theoretical regret bound is established under MMD estimation error and extensive experiments on synthetic functions and real problems demonstrate that our approach can handle various input uncertainties and achieve a state-of-the-art performance.

Efficient Robust Bayesian Optimization for Arbitrary Uncertain inputs
Lin Yang
Junlong Lyu
Huawei Noah's 's Ark Lab China
Huawei Noah's Ark Lab Hong Kong SAR, China Jynjunlong@huawe com com
yangl wanglliin@@hawei.comm
Wenlong Lyu Huawei Noah's ss Ark Lab China
Zhitang Chen
Huawei Noah's Ark Lab Hong Kong SAR, China
www.nlong22huuweei.comm com
chenzhitang200uavei.cco
Abstract
Bayesian Optimization (BO) is a sample-efficient optimization algorithm widely employed across various applications. In some challenging BO tasks, input uncer- tainty arises due to the inevitable randomness in the optimization process, such machining errors, execution noise, contextual variability. This uncertainty deviates the input from the intended value before evaluation, resulting in signifi- cant performance fluctuations in final result. In this paper, introduce a novel robust Bayesian Optimization algorithm, AIRBO, which can effectively identify a robust optimum that performs consistently well under arbitrary input uncertainty. Our method directly models the uncertain inputs of arbitrary distributions by em- powering the Gaussian Process with the Maximum Mean Discrepancy (MMD) and further accelerates the posterior inference via Nyström approximation. Rig- orous theoretical regret bound is established under MMD estimation error and extensive experiments on synthetic functions and real problems demonstrate that our approach can handle various input uncertainties and achieve a state-of-the-art performance.
Introduction
Bayesian Optimization (BO) is a powerful sequential decision-making algorithm for high-cost black- box optimization. Owing to its remarkable sample efficiency and capacity balance exploration and exploitation, BO has been successfully applied diverse domains, including neural architec- ture search (32]. hyper-parameter tuning (4] 29 [2]. and robotic control [18] 5], among others. Nevertheless, in some real-world problems, the stochastic nature of the optimization process, such as machining error during manufacturing, execution noise of control, variability contextual factor, inevitably introduces input randomness, rendering the design parameter to deviate to z' before evaluation. This deviation produces fluctuation of function value y and eventually leads to performance instability of the outcome, In general, the input randomness is determined by the application scenario and can of arbitrary distribution, even quite complex ones. Moreover, some cases, we cannot observe the exact deviated input x' but rough estimation for the input uncertainty. This is quite common for roboties and process controls. For example, consider a robot control task shown in Figure Taj, drone is sent to target location æ to perform measurement task. However, due to the exccution noise caused by the fuzzy control or sudden wind, the drone ends up at location Instead of observing the exact value of
r' P(x) and gets noisy measurement y = f(x')
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Execution naise
Machining error: P()
Measurement noise: W'0003) Yn ((xx) + Sn
E ME x'-P(x)
Xn
x/
In + 6/2 ((x))
11 X2.
B0
Kay Yn.
(a) An example case: drone mea- surement with exccution noisc.
(b) Problem formulation.
Figure 1: Robust Bayesian optimization problem.
r', we only have coarse estimation of the input uncertainty P(z). The goal is to identify a robust location that gives the maximal expected measurement under the process randomness.
To find a robust optimum, it is crucial to account for input uncertainty during the optimization process. Existing works [24 3 T0 along this direction assume that the exact input value, iee, x' in Figure lb is observable and construct a surrogate model using these exact inputs. Different techniques are then employed identify the robust optimum: Nogueira et al. utilize the unscented transform to propagate input uncertainty to the acquisition function (24), while Beland and Nair integrate over the exact GP model to obtain the posterior with input uncertainty [3) Meanwhile, (7 designs robust Sefiiennee-bounder acquisition and applies min-max optimization to identify the robust optimum. Similarly, [1O] constructs an adversarial surrogate with samples from the exact surrogate. These methods work quite well but are constrained by their dependence on observable input values, which may not always be practical.
An alternative approach involves directly modeling the uncertain inputs. A pioneering work by Moreno al. (20) assumes Gaussian input distribution and employs symmetric Kullhack-Leibler divergence (SKL) to measure the distance of input variables. Dallaire et al. [13 implement Gaussian process model with an expected kernel and derive closed-form solution by restricting the kernel to linear, quadratic, or squared exponential kernels and assuming Gaussian inputs. Nonetheless, the applicability of these methods is limited due to their restrictive Gaussian input distribution assumption and kernel choice. To surmount these limitations, Oliveira et al. propose robust Gaussian process model that incorporates input distribution by computing an integral kernel. Although this kernel can be applied to various distributions and offers a rigorous regret bound, its posterior inference requires a large sampling and can be time-consuming.
In this work, we propose an Arbitrary Input uncertainty Robust Bayesian Optimization algorithm (AIRBO). This algorithm can directly model the uncertain input of arbitrary distribution and propagate the input uncertainty into the surrogate posterior, which can then be used to guide the search for robust optimum. To achieve this, we employ Gaussian Process (GP) as the surrogate and empower its kernel design with the Maximum Mean Discrepancy (MMD), which allows us to comprehensively compare the uncertain inputs Reproducing Kernel Hilbert Space (RKHS) and accurately quantify the target function under various input uncertainties (.I). Moreover, to stabilize the MMD estimation and accelerate the posterior inference, we utilize Nyström approximation to reduce the space complexity of MMD estimation from O(m²) to ((m)), where h & (Sec. 3.2). This can substantially improve the parallelization of posterior inference and rigorous theoretical regret bound is also established under the approximation error (Sec. 4. Comprehensive evaluations on synthetic functions and real problems in Seo demonstrate that our algorithm can efficiently identify robust optimum under complex input uncertainty and achieve state-of-the-an -rt performance.
2 Problem Formulation
In this section, we first formulize the robust optimization problem under input uncertainty then briefly review the intuition behind Bayesian Optimization and Gaussian Processes.
2.1
Optimization with Input Uncertainty
As illustrated Figure Tb
we consider an optimization of expensive black-box function: f(xx),
where is the design parameter to be tuned. At each iteration n. we select a new query point x.

---

according to the optimization heuristics. However, due to the stochastic nature of the process, such as machining error or execution noise, the query point is perturbed x', before the function evaluation. Moreover, we cannot observe the exact value of I and only have vague probability estimation of its value: Ps..- After the function evaluation, we get noisy measurement y = f(x(,) + Sn. where Go is homogeneous measurement noise sampled from (0, o²). The goal is to find an optimal design parameter x" that maximizes the expected function value under input uncertainty:
arg max x'
(z')dx' max Ep, [f]
Depending on the specific problem and randomness source, the imput distribution Pz can be arbitrary in general and even become quite complex sometimes. Here we do not place any additional assump- tion on them, except assuming we can sample from these input distributions, which can be easily done by approximating it with Bayesian methods and learning a parametric probabilistic model (T6]. Additionally, we assume the exact values of x' are inaccessible, which is quite common in some real-world applications, particularly in robotics and process control (25).
2.2 Bayesian Optimization
In this paper, we focus on finding the robust optimum with BO. Each iteration of BO involves two key steps: 1) fitting surrogate model and II) maximizing acquisition function.
Gaussian Process Surrogate: To build sample-efficient surrogate,
we choose Gaussian Process (GP) as the surrogate model in this paper. Following (34) GP can be interpreted from a weight- space view: given a set of n observations, D,, (z yi I,.., n). Denote all the inputs as X E IRD Dxn and all the output vector as y E R We first consider a linear surrogate:
f(x) - x" y f(x) + + ~~05C~~0000
where w is the model parameters and is the observation noise. This model's capacity is limited due to its linear form. To obtain a more powerful surrogate, we can extend it by projecting the input æ into feature space Φfx:). By taking Bayesian treatment and placing zero mean Gaussian prior on the weight vector: w ~ (O, Σp), its predictive distribution can be derived as follows (see Section2.1 of (34] for detailed derivation):
fe|x,, ~ No (x)) v
- (z)
where A = (X))X_X and is a identity matrix. Note the predictive distribution also Gaussian and the feature mappings are always in the form of inner product with respect to Ep-
This implies we are comparing inputs a feature space and enables us to apply kernel trick. Therefore, instead of exactly defining feature mapping <(-), we can define kernel: k(x,
<(x)" X(z' = (xx) va("). Substituting it into Eq. gives the vanilla GP posterior:
f- y, X. (u,, E.,, ).), where u. K(X,, X)(K(XXX) + " Σ, K( X.) X(X,,X) XX(KXXXX)
(4)
-(X,
From this interpretation of GP, note that its core idea to project the input x (possibly infinite) feature space (xx) and compare them in the Reproducing Kernel Hilbert Space (RKHS) defined by kernel.
Acquisition Function Optimization: Given the posterior of GP surrogate model, the next step is to decide a query point 2n, The exploitation and exploration balance is achieved by designing an acquisition and adopt the Upper Confidence Bound (UCB) acquisition:
a(x (D.). Through numerous acquisition functions exist (28), we follow 25
((xxDn)) p. (x) + Bo, (x),
(5)
where B is hyper-parameter control the level of exploration.
3 Proposed Method
To cope with randomness during the optimization process, we aim to build a robust surrogate that can directly accept the uncertain inputs of arbitrary distributions and propagate the input uncertainty into
3

---

the posterior. Inspired by the weight-space interpretation of GP, we empower GP kernel with MMD compare the uncertain inputs RKHS. In this way, the input randomness is considered during the covariance computation and naturally reflecedd in the resulting posterior which then can be used to 3.1 Modeling the Uncertain Inputs
guide the search for a robust optimum(Sec, [.T]. To further accelerate the posterior inference, we employ Nyström approximation to stabilize the MMD estimation and reduce its space complexity (Sec. 3.2p
Assume Pa Pxx C P are set of distribution densities over I&d representing the distributions of the uncertain inputs. We are interested in building a GP surrogate over the probability space P, which requires to measure the difference between the uncertain inputs.
To do so, we turn the Integral Probabilistic Metric (IPM) (23). The basic idea behind IPM is to define distance measure between two distributions P' and Q as the supremum over a class of functions G of the absolute expectation difference:
d(P, Q) sup (E)r En-)@(())
where G is class of functions that satisfies certain conditions. Different choices of G lead to various IPMs. For example, if we restrict the function class to be uniformly bounded in RKHS we can get the realizes the Wasserstein distance [14).
MMD [5). while a Lipschitz-continuoon
In this work, we choose MMD as the distance measurement for the uncertain inputs because of its intrinsic connection with distance measurement in RKHS. Given a characteristic kernel k Rd x Rd + R and associate RKHS Mg. define the mean map & P + H4 such that (eb(P).) = Eplg). Vg E H4. The MMD between P. Q defined as:
MMD Q)=
sup [g(u)(() Ev-og(v)] 1/0p- - -ell: Il
Without any additional assumption on the input distributions, except we can get m samples {u, 1m 1m from P, Q respectively, MMD can be empirically estimated as follows [1]: {m
2
MMD² (P,)) m(m -
E 1Si
(k(u,, ug) + k(v,,
Σ
k(u,, Uy),
(8)
To integrate MMD into the GP surrogate, we design an MMD-based kernel over P as follows:
k(P,Q) = exp(-aMMD
with a learnable scaling parameter a. This is a valid kernel, and universal w.r.t. C(P) under mild conditions (see Theorem 2.2. (TT).. Also, it is worth to mention that, to compute the GP posterior, we only need to sample m points from the input distributions, but do not require their corresponding function values.
With the MMD kernel, our surrogate model places prior GP k(P,,, Pa:)) and obtain dataset D., = ((i,, v1 [2 Pa,, = 2,,,.....1) The posterior is Gaussian with mean and variance:
(P.) = (P.) (K, + -01)-
(10)
(P.) = =k(P.. P.) - k,, (P.)) + o
(11)
where yn := [y). yn)² (P.) = [k(PP, P,,), ... , PT.))T and [K,] = h(Pr,, Pay).
3.2 Boosting posterior inference with Nyström Approximation
To derive the posterior distribution of our robust GP surrogate, it requires estimating the MMD between each pair of inputs. Gretton el al. prove the empirical estimator in Eq. approximaas MMD in bounded and asymptotic way (15). However, the sampling size m used for estimation greatly affects the approximation error and insufficient sampling leads to high estimation variance(ref. Figure Ba).

---

Such an MMD estimation variance causes numerical instability of the covariance matrix and propa- gates into the posterior distribution and acquisition function, rendering the search for optimal query point challenging task. Figure 3b giges an example of MMD-GP posterior with insufficient samples, which produces noisy acquisition function and impedes the search of optima. Increasing the sampling size can help alleviate this issue. However, the computation and space complexities of the empirical MMD estimator scale quadratically with the sampling size m. This leaves us with a dilemma that insufficient sampling results in highly-varied posterior while larger sample size can occupy significant GPU memory and reduce the ability for parallel computation.
To reduce the space and computation complexity while retaining stable MMD estimation, we resort to the Nyström approximation (33). This method alleviates the computational cost of kernel matrix by randomly selecting h subsamples from the m samples(h < m) and computes an approximated matrix via K = Kn K2 KT mh. Combining this with the MMD definition gives its Nyström estimator:
MMD* (P) = E x
pkk(u, u')] + Ey 2
( (ie,, v')] - 2 2E um 3 Q (kuu, v)]
2 ITW1. W1,,
T m
U1 m
+ m²
IT m
VIm.
m²
m²
2 - W mh W WtW.WWtWW W mh
1
U
mh U Un +
1 V mn V
VT mb 1m
1 (12)
m²
m²
m
where = = K(u, n', == V(,,v",, W K(,, v) are the kernel matrices, 1. represents a m-by-1 vector of ones, m defines the sampling size and h controls the sub-sampling size. Note that this Nyström estimator reduces the space complexity of posterior inference from O PMMNm' to OmM), where M and N are the numbers of training and testing samples, m the sampling size for MMD estimation while h < m is the sub-sampling size. This can significantly boost the posterior inference of robust GP by allowing more inference to run parallel on GPU.
4 Theoretical Analysis
Assume x E X' C Rd and Pe E Px C P are a set of distribution densities over R2 representing the distribution of the noisy input. Given characteristic kernel A: : TRd x IRd TR and associate RKHS He. we define the mean map v/ P + Hr such that (v²(P,, g) = Er Vgg E He.
We consider a more general case. Choosing any suitable functional L such that k(P, == L v/r) ) is positive-definite kernel over P, for example the linear kernel (Wp, w/r) and radial kernels cxp(- aalll - v Il:) using the MMD distance as a metric. Such a kernel ls is associated with RKHS H₁ containing functions over the space of probability measures P.
One important theoretical guarantee to conduct GP model is that our object function can be ap- proximated by functions in Hg. which relies on the universality of k. Let C(P) be the class of continuous functions over P endowed with the topology of weak convergence and the associated Borel olaggerra, and we define ff E C(P) such that
f(P) = Erlf,
which is just our object function, For l be radial kernels, it has been shown that is is universal w.r. C(P) given that X' compact and the mean map 1b is injective 11 22 For & be linear kernel which not universal, it has been shown in Lemma 1, 26) that f E H2 if and only if f H and further llfll; l|. Thus, in the remain of this chapter, we may simply assume f E 74.
Suppose we have an approximation kernel function i: Q) near to the exact kernel function k(P.Q). The mean P., (p.) and variance 0², (p.) are approximated by
(P.) = kP.) (KK, + -211 (P.) = -k(P.. P.) - k, (P.) (KK, + o
(13)
(14)
where Yn := [v1..... ",,'" (P.) = Pt), PI), .k()),..,,,,,,,....... and Knny k(P.,
5

---

The maximum information gain corresponding to the kernel is denoted as In det o RR,,,
In
sup RE |||
i(y.: f. (R)
Denote
- ((P, Q) as the error function when estimating the kernel k: We suppose e(P, Q) has an upper bound with high probability:
Assumption 1. For any E > 0, P.Q Px, we may choose an estimated k(P, Q) such that the error function e(P,Q) can be upper-bounded by e, with probability at least E, that is, P (le(P, Q)I S cz) >-1)
Remark. Note that this assumption is standard in our case: we may assume m || S Φ, where @ is the feature map corresponding to the k. Then when using empirical estimator, the error with probability at least - according to Lemma 1, [8]. When using the Nyström estimator, the error has a similar form as the empirical one, and under mild conditions, when
between MMDempizica and MMD is controlled by 4 / log(/E))
O(ym) m log(m)), we get the error of the order O(m log( (/)) with probability at least - e. One can check more details in Lemma
Now we restrict our Gaussian process in the subspace Px = (Pa, x E x C C. We assume the observation y; = fx;) + Ç: with the noise G. The input-induced noise is defined as A
f(x:) ) - f(P=:). Then the total noise is y; Er,, [/] C + ASp., : We can state our main result, which gives cumulative regret bound under inexact kernel calculations,
[(x))) - ( [
Theorem 1. Let 0. Hi. and the corresponding ll.fll: S b, maxxsx [f()) M. Suppose the observation noise Çi ys f(z;) is oe--uu--Gaann Gaussian, and thus with high probability [çil < A for some A > 0. Assume that both k and Ps satisfy the conditions for AfP. to be DE-sub-Gaussiae for given 0g > O. Then, under Assumption ] with E > 0 and corresponding ee, setting 0² = 1 1+ 2 running Gaussian Process with acquisition function
n
a(x(Dn)) = (Px) + (Px)
(15)
where B. = (b Vo' + oZ (f.. +1 -1n)))))
we have that the uncertain-inputs cumudative regret satisfies:
Rn E Vnnnnnnn ln + n² (hn - -nnde= + nòes)
(16)
with probability at least 1 1- 5 - NE. Here R. = Σi= Fr. and F. = m Ep, [f]- Er.,
-Ep, [f]
The proof of our main theorem D can be found in appendix B.3
The assumption that Çi is a gg-sub-Gaussian can be met when Pa uniformly bounded or Gaussian, as stated Proposition 3, [26]. Readers may check the definition of sub-Gaussian in appendix, Definition [
is standard in GP fields. The assumption that A fr, is
To achieve an regret of order R. E (V)n).., the same order as the exact Improved gP regret (3
and ensure this with high probability, we need to take = O(8/n). e, = O(n ( A n +)), 1 V n) log(n)) for MCMC approximation, or with a same sample size m and a subsample size h of order O( 9 -- v n for Nyström approximation with some > O. Note that (I6) only offers upper bound for cumulative regret, in real applications the calculated regret may be much smaller than this bound, as the approximation error e, can be fairly small even with few samples when the input noise is relatively weak.
and this requires a sample size m of order On, - (6,
To analysis the exact order of 7n could be difficult, as influenced by the specific choice of embedding kernel k and input uncertainty distributions Pa. I, e X. Nevertheless, we can deduce the following result for a wide range of cases, showing that cumulative regret is sub-linear under mild conditions. One can check the proof in appendix [B.4
Theorem 2 (Bounding the Maximum information gain). Suppose r-t differentiable with bounded derivatives and translation invariant, i.e., k(z, y) = k(x 0). Suppose the input uncertainty is i.i.d, that is, the noised input density satisfies Pa, (z) Po(a xz),, Vx, X. Then if the space X is compact in Rdd the maximum information gain 4/, satisfies
a(d+1)
fn = O(n log(n)).
Thus, when d(d + 1), the accumulate regret is sub-linear respect n, with sufficiently small e

---

OF
INOOT 23000
PADOP 30303
1e00
(a) GP posterior with (b) MMDGP posterior (c) MMDGP posterior with (d) MMDGP posterior Gaussian imput uncertainty with input uncertainty a variance-changing beta with Chi-squared distri- P P= (000000.)
P P0, 0.01). distribution.
bution of changing DoF.
Figure 2: Modeling results under different types of input uncertainties.
5
Evaluation
In this section, we first experimentally demonstrate AIRBO's ability to model uncertain inputs of arbitrary distributions, then validate the Nyström-basee based inference acceleration for GP posterior, followed by experiments on robust optimization of synthetic functions and real-world benchmark.
5.1 Robust Surrogate
Modeling arbitrary uncertain inputs: We demonstrate MMDGP's capabilities by employing an RKHS function as the black-box function and randomly selecting 10 samples from its input domain. Various types of input randomness are introduced into the observation and produce training datasets of D = ((r,, fx; + 6,))]6, Pz; 1 with different Pz configurations. Figure and and 2b compare the modeling results of a conventional GP and MMDGP under a Gaussian input uncertainty Ps = (0, 0. We observe that the GP model appears to overfit the observed samples without recognizing the input uncertainty, whereas MMDGP properly incorporates the input randomness into its posterior.
To further examine our model's ability under complex imput uncertainty, we design the input dis- tribution to follow beta distribution with input-dependent variance: P. beta(a 0.5, 0.5. a = (si 4 + 1)). The MMDGP posterior is shown in Figure 2 As the input variance changes along v, inputs from the left and right around a given location x; yield different MMD) distances, resulting in an asymmetric posterior (e.g., around x = 0.05 and x = 0.42). This suggests that MMDGP can precisely model the multimodality and asymmetry of the input uncertainty.
Moreover, we evaluated MMDGP using a step-changing Chi-squared distribution Pa x² ( 0.01), where g(x) = 0.5 if E 0.6). and g 7.0 otherwise. This abrupt change in g(x) significantly alters the input distribution from sharply peaked distribution to a flat one with a long tail. Figure 2d illustrates that our model can accurately capture this distribution shape variation, as evidenced by the sudden posterior change around x = 0.6. This demonstrates our model can thoroughly quantify the characteristics of complex input uncertainties.
Comparing with the other surrogate models: models under the step-changing Chi-squared input distribution. The results are reported in Figure and they demonstrate our model outperforms obviously under such a complex input uncertainty (see Appendix D.I for more details)
We also compare our model with the other surrogate
5.2 Accelerating the Posterior Inference
Estimation variance of MMD: We first examine the variance of MMD estimation by employing two beta distributions P = beta(a 0.4, = 0.2, o 0.1) and Q = beta(a 0.4, 8 0.2, 0.1) + c, where e is an offset value. Figure 3a shows the empirical MMDs computed via Eq. g with varying sampling sizes as Q moves away from P. We find that a sampling size of 20 is inadequate, leading to high estimation variance, and increasing the sampling size to 100 stabilizes the estimation.
We further utilize this beta distribution P as the input distribution and derive the MMDGP posterior via empirical estimator in Figure 3bb Note that the MMD instability caused by inadequate sampling subsequently engenders fluctuating posterior and culminates in noisy acquisition function, which prevents the acquisition optimizer (e.g., L-BFGG-B in this experiment) from identifying the optima. Although Figure 3c shows that this issue can be mitigated by using more samples during empirical
7

---

Petur
DAEGU TO SCO
1970 605
A
(a) The empirical MMD (b) The noisy posterior de- (c) The posterior becomes (d) The Nyström estima- and covariance values be- rived from sampling size
smoother with sampling tor with less memory
tween two beta distribu- of (upper) traps the size of 100 and acq. op- sumption also produces tions P and Q as Q moves acq. optimizer at æ = 0 timizer can easily identify smooth posterior that is away from P.
(lower).
the optima.
casy to optimize.
Figure 3: The posterior derived from the empirical and Nyström MMD approximators with varying sampling sizes.
Table Performance of Posterior inference for 512 samples.
Method
Sampling Size
Sub-sampling Size
Inference Time (seconds)
Batch Size (samples)
Empirical
20
1.143 + 0.083
512
Empirical
100
8.117 + 0.040
128
Empirical
1000
840.715 2.182 0.780 + 0.001
1
Nystrom 100
10
512
Nystrom 1000
100
21.473 0.984
128
MMD estimation, it is crucial to note that larger sampling size significantly increases GPU memory usage because its quadratic space complexity of OMNN²? (M and N are the sample number of training and testing, is the sampling size for MMD estimation). This limitation severely hinders parallel inference for multiple samples and slows the overall speed of posterior computation.
Table summarizes the inference time of MMDGP posteriors at 512 samples with different sampling sizes. We find that, for beta distribution defined this experiment, the Nyström MMD estimator with sampling size of 100 and sub-sampling size of 10 already delivers comparable result to the empirical estimator with 100 samples (as seen in the acquisition plot of Figure 3d). Also, the inference time is reduced from 1. to 0.78 seconds by enabling parallel computation for more samples. For the cases that require much more samples for MMD estimation (e.g., the input distribution is quite complex high- sigh-dimensiona), this Nyström- based acceleration can have more pronounced impact.
Effeet of Nyström estimator on optimization: optimization, we also perform an ablation study in Appendix [.2] the results in Figure [ suggest that Nyström estimator slightly degrades the optimization performance but greatly improves the inference efficiency.
investigate the effect Nyström estimator on
5.3 Robust Optimization
Experiment setup: To experimentally validate AIRBO's performance, we implement our algorithm
based on Bo'Torch (2] and employ a linear combination of multiple rational quadratic kernels (6) to compute the MMD as Fq. 9 We compare our algorithm with several baselines: 1) uGP-UCB (26) closely related work that employs an integral kernel to model the various input distributions. has a quadratic inference complexity of (M Nm²), OM where M and are the sample numbers of the training and testing set, and m indicates the sampling size of the integral kernel. GUCB is the standard GP with UCB acquisition, which represents a broad range of existing methods that
'The code will be available on https::/github.. wttp://hhaaeenoohh/teeeee noah/HEBB0 and more implementation
details can be found in Appendix C.I

---

(a) The RKHS function
(b) Robust regrets of (c) The double-pen peak func- RKHS function.
(d) Robust regrets of
tion
double-ped peake function
Figure 4: Robust optimization results on synthetic functions.
focus non-robust optimization. 3) SKL-UCB employs symmetric Kullback- divergence to measure the distance between the uncertain inputs (20). Its closed-form solution only exists if the input distributions are the Gaussians. 4) ERBF-UCE is the robust GP with the expected Radial Basis Function kernel proposed in (13). It computes the expected kernel under input distribution using the Gaussian integrals. Assuming the input distributions are sub- Gau-Gausss this method can efficiently find the robust optimum. Since all the methods use UCB acquisition, we simply distinguish them by their surrogate names in the following tests.
At the end of the optimization, each algorithm needs to decide a final outcome oy, perceived to be the robust optimum under input uncertainty at step n. For fair comparison, we employ the same outcome policy across all the algorithms: z = arg max, (2), where fú, (x)) is the posterior mean of robust surrogate at and D,, = {( fI, (; Pa.) are the observations so far. The optimization performance is measured in terms of robust regret as follows:
r(x,,) Es-P,- [f(x" + +]] Es [f 8)].
(17)
where x² is the global robust optimum and x; represents the outcome point at step n. For each algorithm, repeat the optimization process 12 times and compare the average robust regret.
We begin the optimization evaluation with an RKHS function that is widely used in previous BO works [1 24 10 Figure 4a shows its exact global optimum resides at = 0.892 while the robust optimum is around æ 0.08 when the inputs follow a Gaussian distribution NO,, 0.013). According to Figure 4b all the robust BO methods work well with Gaussian uncertain inputs and efficiently identify the robust optimum, but the GP-UCB stagnates at a local optimum due to its neglect of input uncertainty. Also, notice the regret of method decrease slightly slower than uGP works in this low- dimensinaio and Gaussian-input case, but later cases with higher dimension and more complex distribution show method is more stable and efficient.
RKHS function:
1D double-peak function: To test with more complex input uncertainty, we design a blackbox function with double peaks and set the input distribution to be multi-modal distribution P.. bela(a = 0.4, 3 0.2, a 0.1). Figure 4c shows the blackbox function (black solid line) and the corresponding function expectations estimated numerically via sampling from the input distribution (i.e., the colored lines). Note the true robust optimum is around x = 0.251 under the beta distribution, but an erroneous location at z 0.352 may be determined if the input uncertainty is incorrectly presumed to be Gaussian. This explains the results in Figure 4d] the performance of SKL-UCB and ERBF-UCB are sub-optimal due to their misidentification of inputs Gaussian variables, while our method accurately quantifies the input uncertainty and outperforms the others.
10D bumped-bow function: we also extend our evaluation to a 10D bumped-bow function (27 under concatenated circular distribution. Figure D demonstrates AIRBO scales efficiently to high dimension and outperforms the others under complex input uncertainty(see Appendix D.3).
Robust robot pushing: To evaluate AIRBO in a real-world problem, we employ robust robot pushing benchmark from (31), in which a ball is placed at the origin point of a 2D space and a robot learns push it predefined target location (gx. 9y). This benchmark takes a 3-dimensional input (rz, Ty. rt), where Tz: re E [-5, +5] are the 2D coordinates of the initial robot location and r. E [0, 30 controls the push duration. We set four targets in separate quadrants, i.e., gl = (-3, 1-3, 92 = (-3.3), 3), Os (4.3, 4.3), and a "twin" target at g = (5.1.3. 3.0). and describe the input uncertainty via a two-component Gaussian Mixture Model (defined in Appendix D.4)) Following 77, 10). this blackbox benchmark outputs the minimum distance to these 4 targets under squared and linear distances: loss = min(d (g,, 1), d(gg, I), d( d(gs, 1), d(ggg, I)), where d(gr. 1) is the Fuclidean distance between the ball's ending location I and the i-th target. This produces loss landscape as

---

nystum
nystram
ESS
megretmeen
Rogretteett
1.3
10-
0.5
20
40
00
80
100
20
40
60
00
100
stop
(a) Contour of the robot push world
(b) Robust regrets of different algorithms
Figure 5: Robust optimization of the robot push problem.
10ftcaur 24 use
officee 1801 gruLL de33 20
(DS ai 2- trn dec at ae
2420 bPh
Figure 6: The robot's initial locations and push times found by different algorithms
shown in Figure Sal Note that 92 is more robust target than Ji because of its linear-form distance while pushing the ball to quadrant is the best choice as the targets, 93 and gá. match the dual-mode pattern of the input uncertainty. According to Figure our method obviously outperforms the others because efficiently quantifies the multimodal input uncertainty. This can be further evidenced by the push configurations found by different algorithms in Figure 6 in which each dot represents the robot's initial location and its color represents the push duration. We find that AIRBO successfully recognizes the targets in quadrant as an optimal choice and frequently pushes from quadrant III to quadrant I. Moreover, the pushes started close to the origin can easily go far away under input variation, so our method learns to push the ball from a corner with long push duration, which is more robust in this case.
Discussion and Conclusion
In this work. generalize robust Bayesian Optimization to uncertain input setting. The weight- space interpretation of GP inspires us to empower the GP kernel with MMD and build robust surrogate for uncertain inputs of arbitrary distributions. We also employ the Nyström approximation to boost the posterior inference and provide theoretical regret bound under approximation error. The experiments synthetic blackbox function and benchmarks demonstrate method can handle various input uncertainty and achieve state-of-the-ard optimization performance.
There are several interesting directions that worth to explore: though we come to current MMD-baaed kernel from the weight-space interpretation of GP and the RKHS realization of MMD, our kernel design exhibits deep connection with existing works on kernel over probability measures 22 Along this direction, as our thcoretic regret analysis in Section does not assume any particular form of kernel and the Nyström acceleration can also be extended to the other kernel computation, it is possible that AIRBO can be further generalized to a more rich family of kernels. Moreover, the MMD used in our kernel is by means limited to its RKHS realization. In fact, any function class F that comes with uniform convergence guarantees and is sufficiently rich can be used, which renders different realizations of MMD. With proper choice of function class F. MMD can be expressed as the Kolmogorov metric or other Farth-mover distances (15) is also interesting to extend AIRBO with the other IPMs.
7 Acknowledgements
We sincerely thank Yanbin Zhu and Ke Ma for their help on formulating the problem. Also, heartfelt appreciation goes to Lu Kang for her constant encouragement and support throughout this work.
10

---

References
[1] John-Alexander Assael et al. "Heteroscedastic treed bayesian optimisation". In: arXiv preprint arXiv:1410.7172 (2014).
[2] Maximilian Balandat al. "BoTorch: A Framework for Efficient Monte-Carie Bayesian
Optimization". In: Advances in Neural Information Processing Systems 33. 2020. URL: http:
//arxivv.o. orgx abs/1910. 06403
[3] Justin Beland and Prasanth Nair. "Bayesian Optimization Under Uncertainty". In: (2017),
5.
[4] James Bergstra al. "Algorithms for typer-parameter optimization". In: Advances in neural information processing systems 24 (2011).
[5] Felix Berkenkamp, Andreas Krause, and Angela P Schoellig. "Bayesian optimization with safety constraints: safe and automatic parameter tuning robotics". In: Machine Learning (2021), pp. 1-35.
[6] Mikolaj Binkkowki et al. Demystifying MMD GANs. 2021. DOI: 10. 20.48550/ar// arr0v./rxiv...1 1801.
01401 arxiv: 1801.01401 01401
[7] Ilija Bogunovic et al. "Adversarially Robust Optimization with Gaussian Processes". In: Advances in Neural Information Processing Systems. Vol. 31. Curran Associates, Inc., 2018.
[8] Antoine Chatalic et al. "Nyström Kernel Mean Embeddings". In: International Conference on Machine Learning. 2022.
[9] Sayak Ray Chowdhury and Aditya Gopalan. "On Kernelized Multi-armed Bandits". In: International Conference on Machine Learning. 2017.
[10] Ryan B. Christianson and Robert B. Gramacy. Robust Expected Improvement for Bayesian Optimization. 2023. arXiv: 2302.08612
[11] Andreas Christmann and Ingo Steinwart. "Universal Kernels on Non-Standard Input Spaces".
In: NIPS. 2010.
[12] Alexander Cowen-Rivers al. "HEBO: pushing the limits of sample-efficient hyper- parameter optimisation". In: Journal of Artificial Intelligence Research 74 (2022), pp. 1269- 1349.
[13] Patrick Dallaire, Camille Besse, and Brahim Chaib-draa. "An Approximate Inference with Gaussian Process to Latent Functions from Uncertain Data". In: Neurocompuling. Adaptive Incremental Learning Neural Networks 74.11 (201 1), pp. 1945-1955. 10::|00.1066///
neucom. 2010.0000000 2010 09.00.... 024
[14] Arthur Gretton, Dougal Sutherland, and Wittawat Jitkrittum. "Interpretable comparison of distributions and models".. In: Advances in Neural Information Processing Systems [Tutorial] (2019).
[15] Arthur Gretton al. "A kernel two-sample test". In: The Journal of Machine Learning Research 13.1 (2012), pp. 723-773.
[16] Trevor Hastie et al. The elements of statistical learning: data mining. inference, and prediction.
Vol. 2. Springer, 2009.
[17] Thomas Kühn. "Eigenvalues of integral operators with smooth positive definite kernels". In: Archiv der Mathematik 49 (1987), pp. 525-534. URL: https: : //api. semanticscholar.
Leg//orppsII::: 121372638
[18] Ruben Martinez-Cantie et al. "Active policy learning for robot planning and exploration under uncertainty." In: Robotics: Science and systems. Vol. 3. 2007, pp. 321-328.
[19] Seyedali Mirjalili and Andrew Lewis. "Obstacles and difficulties for robust benchmark prob- lems: A novel penalty-based robust optimisation method". In: Information Sciences 328 (2016), pp. 485-509.
[20] Pedro Moreno, Purdy Ho, and Nuno Vasconcelos. "A Kullback- Divergence Based Ker- nel for SVM Classification Multimedia Applications". Advances Neural Information Processing Systems. Vol. 16. MIT Press, 2003.
[21] Krikamol Muandet al. "Kernel Mean Embedding of Distributions: A Review and Beyond'' In: Foundations and Trends in Machine Learning 1-.1-2 (2017), pp. 1-141. DOI: 10. 1561/ 2200000060) arXiv: 1605.09522
[22] Krikamol Muandet et al. "Learning from Distributions via Support Measure Machines". In: NIPS. 2012.
11

---

[23] Alfred Müller. "Integral probability metrics and their generating classes of functions". In: Advances in applied probability 29.2 (1997). pp. 429-443.
[24] Jose Nogueira et al. "Unscented Bayesian Optimization for Safe Robot Grasping". In: 2016 JEEE/RS/ International Conference Intelligent Robots and Systems (IROS). Daejeon, South Korea: IEFE, 2016, pp. 1967-1972. ISBN: 778-1-5090-3762-9 DOI: 10. 1109/IROS 2016.
7759310
[25] Rafacl Oliveira, Lionel Ott, and Fabio Ramos. "Bayesian Optimisation under Uncertain In- puts". In: Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics. PMLR, 2019, pp. 1177-1184.
[26] Rafael Oliveira, Lionel Ott, and Fabio Ramos. Bayesian optimisation under uncertain inputs.
2019. arXiv: 1902. 07908
[27]
Nicholas D. Sanders al. Bayesian Search for Robust Optima. 2021. arXiv: 1904. 11416
[28] Bobak Shahriari et al. "Taking the Human Out of the Loop: A Review of Bayesian Optimiza-
tion". In: Proceedings of the IEEE 104.1 (Jan. 2016), pp. 148-175. DOT: 10. JPRCC
2494218
[29] Jasper Snoek, Hugo Larochelle, and Ryan P Adams. "Practical bayesian optimization of machine learning algorithms". In: Advances in neural information processing systems 25 (2012).
[30] Niranjan Srinivas et al. "Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design". In: International Conference on Machine Learning. 2009. URL: https://api. semanticschos
.rg. Ceg//CrrpuIID. 59031327
[31] Zi Wang and Stefanie Jegelka. "Max-value entropy search for efficient Bayesian optimization". In: International Conference on Machine Learning. PMLR. 2017, pp. 3627-3635.
[32] Colin White, Willie Neiswanger, and Yash Savani. "Bananas: Bayesian optimization with neural architectures for neural architecture search".. In: Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 35. 12. 2021, pp. 10293-10301.
[33] Christopher Williams and Matthias Seeger. "Using the Nyström Method to Speed Up Kernel Machines". In: Advances in Neural Information Processing Systems. Vol. 13. MIT Press, 2000.
[34] Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning. Vol. 2. 3. MIT press Cambridge, MA, 2006.
12

---

A Nyström Estimator Error Bound
Nyström estimator can easily approximate the kernel mean embedding Up pr, as well the MMD distance between two distribution density P1 and p2. We need first assume the boundedness of the feature map to the kernel k:
Assumption 2. There exists positive constant K S oo such that suPzex [|o(x)] - K
The true MMD distance between p1 and pa denoted as MMD(pI.P2) The estimated MMD distance when using Nyström sample size mi, sub-sample size h; for p3 respectively, is denoted as MMD Ak) Then the error
(p..m)) MMD(pI, p2) MMD (MMD(p..m.h)))
and now we have the lemma from Theorem 5.1 in 8
Lemma Let Assumption D hold. Furthermore, assume that for i 1,2, the data points "i,... ".,.., are drawn i.i.d. from the distribution ps and that h2 m, sub-samples Xi, ,X, are drawn uniformly with replacement from the dataset {Xi,'.. X Xmm... Then, for any 8 E (0, 1),
holds with probability at least 28
Vlog(hs/0)
12K log(hs/0) V.
Σ 1=1,2
02
+
vm
h4
h4
provided that, for {1,2},
h4 2 (67, 12K log(h//))
where 61 2K (o 4/33Kogg12/// log (128)) and C4 6K Vlo((12/5). The notation Ne
denotes the effective dimension associated to the distribution pk-
Specifically, when the effective dimension N satisfies, for some 0,
either No (o²) I co2-1 for some y E (0, 1),
or No S log(1 3/32))/,,0 for some
Then, choosing the subsample size m to be
h; m, 1/(2 1og(m./8) in the first case
or h4 mm log( mm max(1/6, m in the second case,
we get Err(pi = 0(1
B Proofs of Section
B.1 Exact Kernel Uncertainty GP Formulating
Following the same notation in Section 4 now can construct Gaussian process GP(O, k) modelling functions over P. This GP model then be applied to learn f from a given set of observations D,, {(P:, bi) Under zero mean condition, the value of f(P.) for given P. P follows Gaussian posterior distribution with
(P.) (P.) + +021) yy
(18)
(P.) = -k(P., P.) - k (P.) + =21 -
(19)
where yn == [...... V. k., =
k(P,, PT) and Kany k(P.,
13

---

Now we restrict our Gaussian process in the subspace Px [P,. E 2 C P. We assume the observation y/ = fz,) + ç with the noise Ç. The input-induced noise is defined as Afpe, := f(x;) Ep, = (PPa)) Then the total noise is y/ Ep., [f]= A Ap.,- To formulate the regret bounds, we introduce the information gain and estimated information gain given any {P:) C:
f(y.; (.PPP In IT + o ²K,),
(20)
Fly.; fn ((P del - ²kn),
(21)
and the maximum information gain is defined 9n SUPREPM|R|=n f(y:: (R). Here f, (f())).... ["
We define the sub- nub-Gaussian condition as follows:
Definition 1. For a given og >> O, real-valued random variable & is said to be oe -e- Gaussian if:
YA E R. E[c S eAPo2/2
Now state the lemma for bounding the uncertain-inpe -inputs regret of exact kernel evaluations, which is originally stated in Theorem 5 in (26).
Lemma 2. Let E 1), E Hk. and the corresponding II S b. Suppose the observation noise y: [xx,) is conditionally ac- e-subbGaussian Assume that both i and Pg satisfy the conditions for Afp. to be Gaass-Gasssan,, for a given > O. Then, we have the following results:
The following holds for all x E X and 1:
(P.) b+ Vo ( ( +1+ (n(//)))))))))))))))))))
(n.(P))
(22)
Rumning with upper confidence bound (UCB) coquisitionnfunction a(x[Dn)) = (P,) + (Pa) where
Bn bb + +1
and set a² 11 2/n, the uncertain-inputs cumulative regret satisfies:
R, E OlVnn V1n + In(1/8))
(23)
with probability least 8
Note that although the original theorem restricted to the case when k(p, (wbp, tho) ks the results can be easily generated to other kernels over P, as long as its universal w.r.t C(P) given that X' is compact and the mean map injective 22
B.2 Error Estimates for Inexact Kernel Approximation
Now let us derivative the inference under the introduce of inexact kernel estimations.
Theorem 3. Under the Assumption 1 for
O, let Jin, on, i(y.; f, P.
defined in
B 14 (1) respectively, and fin. In.I(yn : (. as defined in .
19
20
Assume
max f(x) = M, and assume the observation error S; = y, - fz(z) satisfies
A for all i.
Then have the following error bound holds with probability least 1 Ne:
n²
(P.) - (P.)
+
(M + A]ee - O
(24)
-2
(o.( - -P. < (1+
a
ye + 0
(25)
02
3//2
&(P,1. fly,: Pa (al((1)711)
202 e, + O
(26)
14

---

Proof. Denote e(P.. Q) Q) - k(P..Q) Q). (P.) = (i), (.PPPP))....e(PPPPP.)'' and
Enli, e(P, P)). Now according to the matrix inverse perturbation expansion.
(X - +010x)1)
we have
(Kn + + (Kn + (Kn + 1211--'B (K, + + O(J)E,
thus
Pn
(P.) (K) yn (P,, (K,,, + +
yn
=,, +((ee, |E.ll)
(P.)) (h + (e.
P(P.. 2e,, (P)()(()) (K,, + 0011 - kn
+k.(P)" (K, + E (K, + - (P.) +O([ + O(Ile,)l [x,I) + O(1e,, -IE,,]))
Notic that the following holds with a probability least ne, according to the Assumption |
len (P.) + o²1) "
ynl s (22) 12 112lly 1
(M A)es,
ynl (P.) + d'I)-''E, (K, - S |I,, (P.) 1221(1 - 112 E,, 1/2lllylll 1212
n2
S Vno - nee V + +A)
- (M + A),
here the fact that K,, semi-definite (which means [I + 1123 , P.) $ lui - M + A. Combining these results, we have that
n (P.) - (P.) < 02 (M +)c,+
n +
holds with probability at least - NE.
Similarly, we can conduct the same estimation to e, (P) (Kn, + and ka (P) + o²I)- (K, + - and get
(P.)) 6.. oP)) (1+ n ) +
holds with probability at least nE.
It remains estimate the enror for estimating the information gain. Notice that, with probability least 1 - ne,
det (I + -
-22K,Kk)) 2k₂)
[yy.: Pxlppe f(y.; (pt))t)))))
log
det(I + 0 -KK,)
log det (02I + Kn,) -E,)))
Tr(log (0²I K,,)K,) k, -1E,))
1 Tr(- ( ( - -'E,)) + 0 E, 11
3/2
202 + O,
15

---

here the second equation uses the fact that det(AB-)) = (A) det(B)) and the third and fourth
equations use log de A) Trlog(1 A) Tr(A 4 ....). The last inequality follows from the fact
Tr(o'I - & II R,, ) ||= EnllE & 0 e
and K, is semi-definite.
With the uncertainty bound given by Lemma 3 let us prove that under inexact kernel estimations, the posterior mean concentrated around the unknown reward function f
Theorem 4. Under the former setting as in Theorem 3 with probability at least 1 6 -ns, let
dv oo + or laking a =1+ the following holds for all x E X:
(Px) f(P+)) S3,0 Bu + n) + (n + n²) (M + A)ec:
(27)
where Bn
ln) + 1)
Proof. According to Lemma
equation
22
we have
(Px) f(Pa) s (non(Pa)
with
3n = -bbboy/2
(1 S An.
Notice that
f(Pa)..
(28)
We also have
25 which means
( + -1++ n)e
(29)
combining
24
and 129
we finally get the result in
B.3
Proofs for Theorem
Now we can prove our main theorem
T
Proof of Theorem 7 Let maximize f(P)) over X. Observing that at each round 1, by the
choice of H., to maximize the aquisition function a(x[2n,1) = -a- (PP) + 9n-1ō (Pa). we have
F, -
S un-1 (P..) + 3n (P- - (n-1((a.)) + (n-1ōn--(Pa.) + 2Err(n 1,
< 2
2Err(n 1,
Here denote Err(n, e = (3n ān (Pu 12 + (n n²) (M + A)ez. The
second inequality follows from 27
i(Pa-) Nn (Pa-) 5 (Ps.) + Err(n - 1, 1,6,)
(P_.(P2))-
Bn-1ôn-1 (P...) Err(n 1, ee)c
and the third inequality follows from the choice of vn:
(n-1(Pa-)) + 3n-1ōn-1(P2) < pn (Pa..) Bn-1ōn-1 (Pa.).
16

---

Thus we have
1
Rn Σ
re S 23n
Σ
8 01-1
(P,,) +
Σ t=1
Err(t 1,),).
1=1
t=1
From Lemma 4 in (9], we have that
1
n
Σ oP.. s V (( o- - I +
2 ec).
and thus
σ -0 Vngn(mn lnd) + V Vnāce + vna(tn In dyoee)
2,3n
(30)
On the other hand, notice that
L,e (Gr In Ino)e
(fn
+ (n²² + + (n²2
(31)
we find that the eg term in 30 can be controlled by in 31,. thus we immediately get the result.
B.4 Proofs for Theorem
Proof. Define the square of the MMD distance between PA, Pa. as dlr (11, have
das(#1,
Mx, (z) (r')drdz' + f k(x, (t) (x')dddd
-2 { ((x,x')P (3:) Pr., (2' (x')drda''
na
[ (k(2 - x1, x1) + k(x 22, - x2) - -2k(x - 21,2'
not hard verify that dM shift invariant: dw (21, - dM 22 0), and dM has r-th bounded derivatives, thus (*1 522) = hi Pz,: Pa,) exp(- (#1, x²2) is shift invariant with r-th bounded derivatives. Then take f(z)) as the Lebesgue measure over X, according to Theorem 4, symmetric compact operator
[77, the integral operator Tk.p fx) = fx K (f,
(X.,), and the spectrum of T satisfies
An (Te, O(n -1-/dd.
((di1)
Then according to Theorem in [30]. have 9/2 O(n ,+)()+1] log(n)), which finish the proof.
C Evaluation Details
C.1 Implementation
In our implementation of AIRBO, we design the kernel & used for MMD estimation be linear combination of multiple Rational Quadratic kernels as its long tail behavior circumvents the fast decay issue of kernel 6):
(a x
k(z,
Σ
(1+
(32)
2a,12
0:E(0.220.5.12253)
where is a learnable lengthscale and a; determines the relative weighting of large-scale and small-scale variations.
17

---

Depending on the form of input distributions, the sampling and sub-sampling sizes for Nyström MMD estimator are empirically selected via experiments. Moreover, as the input uncertainty already modeled in the surrogate, we employ classic UCB-based acquisition as Eq. with B = 2.0 and maximize via an L-BBFGSB optimizer.
D
More Experiments
D.1
Comparing with the Other Models
To compare the modeling performances with the other models, we design the input uncertainty to follow a step-changing Chi-squared distribution: Ps = x²(g(x), a 0.01), where g(æ) = 0.5 if x [0.0, 0.6) and g(z) 7.0 when E [0.6,1.0j. Due to this sudden parameter change, the uncertainty at point z = 0.6 is expected to be asymmetrie: 1) on its left-hand side, as the Chi-squared distribution becomes quite Iean and sharp with small value of 9(x)
0.5, the distance from z = 0.6 to its LHS points, Liis E [0.0.0.6). are relatively large, thus their covariances are small, resulting fast-growing uncertainty. 2)eenne,,,,, when E [0.6, 1.0), the g(x) suddenly increases to 7.0, rendering the input distribution a quite flat one with a long tail. Therefore, the distances between 0.6 and its RHS points become relatively small, which leads large covariances and small uncertainties for points in [0.6, 1.0]. As a result, we expect to observe an asymmetric posterior uncertainty at 0.6.
Several surrogate models are employed in this comparison, including:
MMDGP-nystrom
is our method with a sampling size m 160 and sub-sampling size h = 10. Its complexity is O(MNmh) where M and N are the sizes of training and test samples (Note: all the models in this experiment use the same training and testing samples for fair comparison).
uGP(40) the surrogate from 25], which employs an integral kernel with sampling size = 40. Due to its O(M (MMNm) complexity, we set the sampling size = 40 to ensure
similar complexity as ours.
uGP(160) is also the surrogate from (25 but uses a much larger sampling size (m 160). Given the same training and testing samples, its complexity is 16 times higher than MMMDGP-aystrom((I0)))
skl is a robust GP surrogate equipped with a symmetric KL-based kernel, which described in 220
ERBF 13 assumes the input uncertainty to be Gaussians and employs a close-form expected RBF kernel.
GP utilizes a noisy Gaussian Process model with a learnable output noise level.
According Figure 7 our method, MMDGP- (orooomm66611000 160/10), can comprehensively quantify the sudden change of the input uncertainty, evidenced by its abrupt posterior change at = 0.6. However, Figure 7b shows that uGP(40) with the same complexity fails model the uncertainty correctly. We suspect this is because uGP requires much larger samples to stabilize its estimation of the integral kernel and thus can perform poorly with insufficient sample size, so we further evaluate the uGP(160) with a much larger sampling size m 160 in Figure 7c It does successfully alleviate the issue but also results in 16 times higher complexity. Apart from this, Figure 7 suggests the noisy GP model with a learnable output noise level is not aware of this uncertainty change at all it treats the inputs as the exact values instead random variables. Moreover, Figure 77 and show that both the skl and ERBF fail in this case, this may be due to their misassumption of Gaussian input uncertainty.
D.2 Ablation Test for Nyström Approximation
In this experiment, we aim to examine the effect of Nyström approximation for optimization. To this end, we choose to optimize an RKHS function (Figure 4a) under a beta input distribution: P. beta(a 0.4, B 0.2, a 0.1). Several amortized candidates include:
18

---

#MUGPinystromi LeOOIOI
/s=0,215. XE OS 0.723 no.se00.035
XC 1s=0.054, XC 0570.707 noise= 0.035
XC
Is.00.599..... 0S-0.392 noice= 0.063
-0--
Nedd teoce pratetton
Boode model predecdo
40
-00
-20
-20
(a) MMD-GP with Nystrom es- timator, in which the sampling size 160 and sub-sampling size 10.
(b) uGP model that uses an inte- gral kernel (26] and a sampling size of m 40.
(c) uGP with an integral ker-
nel
26 and uses much larger
samping size (m 160).
GP
ski
ER3F KC. 3.747 xc oo 1.201, noisa 0.127
Nc. (0..223,, xc_ oun 1,000, noise-0.055
xc. 1.....1..7. 782, ac_ 0s-0111 noise-0.0008
profition
mican
2500
campeee
get:
Hoddl producion
40
moon
230d
sampaee
MoOd profition
90
25ti mein
20-
20
20
20
-20
(d) Conventional noisy GP model.
(e) GP model with a symmetric KL-divergence kernel (20)
(I) Robust GP model with an ex- pected RBF keree [13
Figure Modeling performance with step-changing Chi-squared distribution.
MHDDG
HHDGT
MMDGR S
201
121
Regrettmeen)
Rogrensed
10P
100
20-1
30"
10:1
20
-2
60
100
20
60
Figure 8: Ablation test for the Nyström approximation.
MMDGP- -nynrrom is our method with Nystrom approximation, in which the sampling size m = 16 and sub-sampling size = 9. Its complexity is ((M Nmh), where M and N are the sizes of training and test samples respectively, is the sampling size for MMD estimation, and indicates the sub-sampling size during the Nystrom approximation.
MMDGP-raw-S does not use the Nystrom approximation but employs an empirical MMD estimator. Due to (MM)²) complexity, we set the sampling size m 12 to ensure similar complexity the MMDGP-nystrom.
MMDGP-rant raw-L also uses empirical MMD estimator, but with larger sampling size (m = -16).
GP utilizes a vanilla GP with learnable output noise level and optimizes with the upper- confidence-boume acquisitior
According Figure B we find that 1) with sufficient computation power, the MMDGP-raw-Le can obtain the best performance by using a large sample size. 2)However, with limited complexity, the performance MMDGP-raw-S degrades obviously while the MMDGP-nystrrm performs much better. This suggests that the Nyström approximation can significantly improve the efficiency with a mild cost of performance degradation. 3) All the MMDGP-based methods are better than the vanilla GP-UCB.
Forr a fair comparison, all the methods in this test use UCB acquisition.
19

---

300
BRBF
systrom
ERBF
DOP
USP
03
Regretimecn))
HegresGoot
10-
20
-00
00
100
20
40
00
80
309
stop
stop
The input GMM distribution.
Figure 10:
Figure 9: Optimization regret on 10D bumped-bowl problem.
MGP
WERBE
62
40
303
125
61
4.0
Figure 11: Simulation results of the push configurations found by different algorithms.
D.3 Optimization on 10D Bumped-Bon -Bowl Problem
To further evaluate AIRBO's optimization performance the high-dimensional problem, employ a 10-dimensional bumped bowl function from 27 19
{
g(x) log (0.8||z/|2 + -100|=11) + 2.54
f(x) = (( where
(33)
h())) Σ: 5x +
Here, is the i-th dimension of æ, 21.2 represents the first dimensions for the variable, and az indicates the rest dimensions. The input uncertainty is designed to follow a concatenated distribution of 2D circular distribution(r 0.5) and multivariate normal distribution with zero mean and diagonal covariance of 0.01.
Figure shows the mean and std values of the optimization regrets. We note that 1 w)we it comes to a -igh-dimessional problem and complex input distribution, the misassumption of Gaussian input uncertainty renders the sk/ and ERBF fail to locate the robust optimum and get stuck at local optimums. 2)rur method outperforms the others and can find the robust optimum efficiently and stably, while the uGP with similar inference cost suffers the instability caused by insufficient sampling and stumbles over iterations, which can evidenced by its high std values of optimization regret.
20

---

D.4 Robust Robot Pushing
This benchmark is based on Box2D simulator from [1] where our objective identify a robust push configuration, enabling robot to push a ball to predetermined targets under input randomness. In our experiment, we simplify the task by setting the push angle to ra arctan ix ensuring the robot is always facing the ball. Also, we intentionally define the input distribution two-component Gaussian Mixture Model as follows:
0.12
-0.2
le
O -1 1 E
0.5 0.5
(T)) GMM
0.3
0.12
le 6
) (34)
le
1
1.02
where the covariance matrix E is shared among components and w is the weights of mixture components. Meanwhile, as the SKL-UCB and ERBF-UCB surrogates can only accept Gaussian input distributions, we choose to approximate the true input distribution with Gaussian. As shown in Figure [0] the approximation error obvious, which explains the performance gap among these algorithms in Figure 5b
Apart from the statistics of the found pre-images in Figure 6 also simulate the robot pushes according to the found configurations and visualize the results in Figure []] In this figure, each black hollow square represents instance the robot's initial location, the grey arrow indicates the push direction and duration, and the blue circle marks the ball's ending position after the push. We can find that, as the GP-UCB ignores the input uncertainty, it randomly pushes these targets and the ball ending positions fluctuate. Also, due to the incorrect assumption of the input distribution, the SKL-UCB and ERBF-UCB fail control the ball'ss 's ending position under input randomness. On the contrary, AIRBO successfully recognizes the twin targets in quadrant as an optimal choice and frequently pushes to this area. Moreover, all the ball's ending positions are well controlled and centralized around the targets under input randomness.
21


Reviews:
[{"review_id": "X38K2tfuws", "reviewer": "Reviewer_jZyz", "summary": "The paper tackles the problem of robust Bayesian Optimization (BO) with uncertain inputs, i.e., the input values are deviated from the intended value before evaluation. The paper proposes a new technique, namely AIRBO (Arbitrary Input uncertainty Robust Bayesian Optimization), that can model the uncertain input and incorporate this uncertainty into the surrogate model, and thus can be used to guide the search of the objective function global optimum. The paper further proposes to use Nystrom approximation to reduce the computational cost. Theoretical analysis is conducted to guarantee the performance of the proposed technique. Experiments are conducted on some synthetic and one real-world problem to evaluate the performance of the proposed method.", "questions": "Please answer my comments and questions in the Weaknesses section.", "limitations": "I can’t find any dedicated section that describes about the limitations of the proposed technique. There are some future work mentioned in the Conclusion, but there are no limitations mentioned there.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "+ The paper’s writing is generally clear and easy to understand. Illustration and figures are plotted nicely to explain the problem setting, as well as the property of the proposed method.\n+ The proposed method of using MMD to construct a kernel that can incorporate uncertainty of the inputs seems to be interesting to me.\n+ Theoretical analysis is conducted to guarantee the convergence property of the proposed method (although note that I have some concerns regarding the theoretical analysis).\n+ The experiments are described very detailed. Various experiments are conducted to understand the performance of each component within the proposed method.\n", "weaknesses": "+ I think the application of this problem setting should be motivated much better as it’s unclear on the significance of the problem tackled in this paper. In the Introduction, the paper only mentions that this problem is quite common for robotics and process controls, but no references or further explanations are provided.\n+ Some notations are quite confusing, making it hard for me to follow the proposed method. In Eq. (8), the paper shows how to compute MMD(P,Q) from the m samples {x_i}_1^m, {y_i}_1^m – but y is defined as the objective function value. Then how can we compute the kernel \\hat{k}(P_{x_i}, P_{x_j}) in Eqs. (10) and (11)? \n+ I also have some concerns with the theoretical analysis. In the theoretical analysis, there is no formal proof to prove the proposed kernel is actually a valid kernel. There is just one paragraph in Lines 170-173 briefly mentioning about this but no formal proof is given to know if the proposed kernel is valid, and in which scenarios it will be a valid kernel. Assumption 1 is too strong, it assumes that we already know that the error function e(P,Q) can be uniformly upper-bounded. In practice, can this assumption be true? How can we know it can be true? Theorem 1 also doesn’t have much meaning to me. In Theorem 1, \\tilde{\\beta}_n is defined with the maximum information gain \\tilde{I} within the formula. This is not really a standard analysis in BO’s theory, normally \\beta is defined by some constants. Is this maximum information gain \\tilde{I} bounded? \n+ Experiments are mostly conducted with synthetic functions. There is only one real-world problem being used in the experiment, and it’s quite simple to me. It’s just a 3-dim problem.\n+ What is the time cost of the proposed method?\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for the constructive comments.  **The reference is listed in the global rebuttal**.\n\n### 1. Motivation\n> I think the application of this problem setting should be motivated much better as it’s unclear on the significance of the problem tackled in this paper. \n\nThanks for bringing it up. This paper aims to tackle a robust optimization problem in which the design parameters are perturbed randomly before the evaluation. Such perturbations are the manifestation of uncertainties in the design process and can arise for several reasons: *e.g*., execution noise during the control process or machining error of manufacturing. The drone measurement in Section 1 gives an example of execution noise, yet some other motivating applications include:  robot exploration [7], robot grasping [3], semiconductor design [4], and more applications can be found in [6].\n\n### 2. Notations\n> In Eq. (8), the paper shows how to compute MMD(P, Q) from the m samples ${x_i}^m, {y_i}^m$ – but y is defined as the objective function value. Then how can we compute the kernel $\\hat{k}(P{x_i}, P{x_j}$) in Eqs. (10) and (11)?\n\nSorry for the abused notation. In Eq.8, the $y$ is not the objective function value, we use $x$ and $y$ to represent the samples from the input distributions $P$ and $Q$ respectively.  The updated Eq.8 should be: \n\nSuppose we can get $m$ samples, $\\\\{ x_i^{(u)} \\\\}\\_{u=1}^m$ and $ \\\\{x_j^{(v)} \\\\}\\_{v=1}^m$, from input distribution $P_{x_i}$ and $P_{x_j}$, the MMD can be empirically estimated via:\n$$\n\\text{MMD}^2(P_{x_i}, P_{x_j}) \\approx \\frac{1}{m(m-1)} \\sum_{1\\leq u, v\\leq m, u\\neq v} \\big( k(x_i^{(u)}, x_i^{(v)}) + k(x_j^{(u)}, x_j^{(v)}) \\big) - \\frac{2}{m^2} \\sum_{1\\leq u,v \\leq m} k(x_i^{(u)}, x_j^{(v)}),\n$$\nwhere $u$ and $v$ are the sample indices and $x_i^{(u)}$ represents the $u$-th sample from $P_{x_i}$, which are consistent with the notations used in Eqs 9, 10, and 11.\n\n### 3. Theoretical analysis\n> In the theoretical analysis, there is no formal proof to prove the proposed kernel is actually a valid kernel. \n\nFor the theoretical analysis, we complement the evidence that our proposed kernel is a valid kernel: by Theorem 2.2 in [11], any type of bi-functions defined on $\\mathcal{P}\\times\\mathcal{P}$ as $\\hat k(P,Q) = \\sum_{i=1}^{\\infty} a_i  \\langle \\psi_P, \\psi_Q \\rangle_k$ with $a_i \\ge 0$, $\\psi_P,\\psi_Q$ is the kernel mean map from $\\mathcal{P}$ to $\\mathcal{H}_k$, then $\\hat k(P,Q)$ is a valid kernel; In addition if $a_i > 0$, this kernel is universal over $C(\\mathcal{P})$, given the mean map $\\psi$ is injective and the space of $P, Q$ defined on is compact. Thus, the Gaussian RBF kernel $\\hat k(P,Q) = \\exp(-\\alpha \\Vert \\psi_P - \\psi_Q \\Vert_k^2)$ is ensured to be valid and universal.\n\n> Assumption 1 is too strong, it assumes...can this assumption be true? How can we know it can be true? \n\nFor the concern on assumption 1, we admit that the original assumption is too strong, and we weaken the uniform condition, only assuming that each evaluation of $e(P,Q)$ has a bound with probability $1-\\varepsilon$. The modified assumption 1 is stated in the global rebuttal.\n\nNote that this assumption is standard in our case: we may assume $\\max_{x \\in \\mathcal{X}} \\Vert \\phi \\Vert_k \\le K $, where $\\phi$ is the feature map corresponding to the $k$. Then when using an empirical estimator, the error between $\\text{MMD}_{\\text{empirical}}$ and $\\text{MMD}$ is controlled by $4K\\sqrt{2\\log(6/\\varepsilon)m^{-1}}$ with probability at least $1- \\varepsilon$ according to Lemma E.1, in [5]. When using the Nystrom estimator, the error has a similar form as the empirical one, and under mild conditions, when $h = O(\\sqrt{m}\\log(m))$, we get the error of the order $O(m^{-{1/2}}\\log(1/\\varepsilon))$ with probability at least $1-\\varepsilon$. \n\nUnder this modified assumption, the result for our theorem 1 is also slightly modified, and we state it in the general Author Rebuttal page. We may see that the final regret is bounded in a probability $1 - \\delta - n\\varepsilon$ compared to the original $1 - \\delta - \\varepsilon$, but it matters a little: as the empirical estimator and Nystrom estimator have errors that are log-scale on $\\varepsilon$, we can take $n\\varepsilon$ small enough without hurting the error order.\n\n> Theorem 1 also doesn’t have much meaning to me. In Theorem 1, \\tilde{\\beta}_n is defined with the maximum information gain \\tilde{I} within the formula. This is not really a standard analysis in BO’s theory, normally \\beta is defined by some constants. \n\nFor theorem 1, the $\\tilde{\\beta}_n$ is defined with the maximum information gain $\\hat\\gamma_n$ because we consider in a general case: for arbitrary input f with bounded RKHS norm, not assuming $f$ is sampled from GP process. One may checkin Theorem 2, [10] for $f$ sampled from GP process, and Theorem 3, [10] for arbitrary $f$ in RKHS space.  for $f$ assume to be sampled from the GP process, we do not need to define $\\beta_n$ with $\\gamma_n$, while for arbitrary $f$ in RKHS space, we need it. \n\n> Is this maximum information gain \\tilde{I} bounded?\n\nWe add the theorem 2 in the global rebuttal to discuss the order of the maximum information gain $\\gamma_n$, showing the sublinear order of the accumulative regret bound. \n\n### 4. Experiments\n> Experiments are mostly conducted with synthetic functions. There is only one real-world problem being used in the experiment, and it’s quite simple to me. It’s just a 3-dim problem.\n\nWe appreciate the reviewers' suggestion, we have added an HD optimization problem with a 10-dimensional bumped bowl problem from [8] and set the input uncertainty as a circular distribution. Figure 7b in the supplementary results shows our method yields the best performance and can locate the robust optimum in high dimension efficiently and stably.\n\n### 5. Time\n> What is the time cost of the proposed method?\n\nTable 1 in the paper reports the time cost for our method.\n"}, {"review_id": "akFyHNCoak", "reviewer": "Reviewer_Yh2J", "summary": "The paper proposes a novel Bayesian Optimization (BO) algorithm that explicitly tackles input uncertainty by introducing a new integral probabilistic metric (IPM)-based kernel. The algorithm also utilizes an efficient and stable Nystrom estimator to approximate the Maximum Mean Discrepancy (MMD), which serves as the adopted IPM. Furthermore, the paper extends the GP-UCB framework by incorporating the proposed kernel and derives the corresponding upper bound for the cumulative regret. The empirical study illustrates the effectiveness of the proposed approximation and the BO algorithm.", "questions": "1. What is SKL-UCB mentioned in line 262? Is it referring to the algorithm discussed in line 45 from [19]?\n2. What is the corresponding reference for Oliveira et al. mentioned in line 50?\n3. Despite the explicit treatment of input uncertainty, it is unclear to me whether tolerating the input uncertainty by fixing a larger output noise level could provide a simpler solution. The author claims that the GP model overfits (line 222) and justifies the use of MMD-GP. However, one would expect that fixing a larger noise level for the GP model would help mitigate it as well. Additionally, the existing method addressing heteroscedastic noise proposed in [1] potentially aids in mitigating the input uncertainty purely on the output side. It would be valuable to hear the authors' comments on this.\n4. Could the author provide more intermediate quantitative evidence showing the effectiveness of the proposed MMD-based kernel? For example, any evidence showing that the uncertainty quantification is better than the naïve GP model that doesn't explicitly deal with the input uncertainty?", "limitations": "LImitations are mentioned in the comments above.", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. The paper is generally well-organized and illustrative with figures and tabular results.\n2. The theoretical result extends the existing UCB work with reasonable approximation to tackle the input uncertainty.\n3. The limitation of the proposed AIRBO on dealing with high-dimensional input and lack of discussion over other IPMs is clearly stated.", "weaknesses": "1. The theoretical bound is arguably sound. Due to the numerical approximation of the Maximum Mean Discrepancy (MMD), it may potentially result in pseudo-metric or even worse outcomes in practice. Therefore, it remains unclear whether the upper bound of the maximum information gain $\\gamma_T$ still holds and guarantees that the UCB cumulative regret is sublinear.\n\n2. The essential justification for the proposed Nystrom estimator lies in assumption 1. While the author claims in lines 210 and 211 that the approximation error can be fairly small, the paper lacks sufficient discussion or substantiation to support this claim.\n\n3. The empirical results are limited. The regret curve lacks statistical significance, and more extensive empirical studies on real-world applications would be desirable.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for the constructive comments.\n\n### Theoretical bound\n> The theoretical bound is arguably sound. Due to the numerical approximation of the Maximum Mean Discrepancy (MMD), it may potentially result in pseudo-metric or even worse outcomes in practice. Therefore, it remains unclear whether the upper bound of the maximum information gain $\\gamma_T$ still holds and guarantees that the UCB cumulative regret is sublinear.\n\nWe modify the expression of our key result, showing that the cumulate regret is bounded by the approximating error $e_\\varepsilon$, and the maximum information gain $\\hat\\gamma_T$ of the **exact kernel $\\hat k$**, thus no need to calculate the approximate information gain. One can check Theorem 1 in our general rebuttal part. We also discuss the order of $\\hat\\gamma_T$ in Theorem 2, showing that in mild assumption, $\\hat\\gamma_T$ satisfies an order depending on its differentiability and dimension, and thus the UCB cumulative regret is sublinear.\n\nThe modified theorem 1 can be checked in the global rebuttal. \n\n\n### Nystrom errors\n> The essential justification for the proposed Nystrom estimator lies in assumption 1. While the author claims in lines 210 and 211 that the approximation error can be fairly small, the paper lacks sufficient discussion or substantiation to support this claim.\n\nWe have updated assumption 1 and added a remark after assumption 1 to discuss the general approximation error of the Nystrom estimator. Please refer to the global rebuttal.\n\n### Experiments\n> The empirical results are limited. The regret curve lacks statistical significance, and more extensive empirical studies on real-world applications would be desirable.\n\nThanks for pointing out this issue. We have rerun all the experiments with more trials to provide a statistically-meaningful result: Figure 9 in the attached supplementary_results.pdf reports all the updated experiment results. In particular, a performance comparison on the RKHS function with Gaussian input uncertainty is reported in Figure 9a, the result on the 1D double-peak function with beta inputs is presented in Figure 9b, while 9c and 9d compare the performance on a real-world robot pushing problem.\n\nIn addition, we also add an HD optimization problem with a 10-dimensional bumped bowl problem from [2] and set the input uncertainty as a circular distribution. Figure 7b in the supplementary results shows our method yields the best performance and can locate the robust optimum in high dimension efficiently and stably.\n\n\n### Questions\n> 1. What is SKL-UCB mentioned in line 262? Is it referring to the algorithm discussed in line 45 from [19]?\n\nYes, the SKL-UCB employs a GP surrogate equipped with a symmetric KL-based kernel, as described in[19], to model the uncertain inputs. We have updated the citation in Iine 262.\n\n\n> 2. What is the corresponding reference for Oliveira et al. mentioned in line 50?\n\nSorry for the mistake. We have added the correct reference to [1] in the updated paper.\n\n\n> 3. Despite the explicit treatment of input uncertainty, it is unclear to me whether tolerating the input uncertainty by fixing a larger output noise level could provide a simpler solution. The author claims that the GP model overfits (line 222) and justifies the use of MMD-GP. However, one would expect that fixing a larger noise level for the GP model would help mitigate it as well. Additionally, the existing method addressing heteroscedastic noise proposed in [1] potentially aids in mitigating the input uncertainty purely on the output side. It would be valuable to hear the authors' comments on this.\n\nThis is an interesting viewpoint. In our understanding, explicit modeling of the input uncertainty not only helps to distinguish the input perturbation from other sources of randomness (e.g., measure noise) but also enables us to predict the **expected** function value under such an input uncertainty, which then can be used to guide the search for the robust optimum. This is quite reasonable and intuitive in our problem setting. On the contrary, simply setting a larger output noise to absorb the input uncertainty may mix the different sources of randomness together, rendering the output noise level hard to learn.\n\n\n> 4. Could the author provide more intermediate quantitative evidence showing the effectiveness of the proposed MMD-based kernel? For example, any evidence showing that the uncertainty quantification is better than the naïve GP model that doesn't explicitly deal with the input uncertainty?\n\nThanks for this constructive comments. We have conducted a new experiment for modeling the uncertainty with diff.\nIn this test, we intentionally design the input uncertainty to be a \"step-changing\" Chi-squared distribution, whose degrees of freedom parameter $df$ is 0.5 when $x \\in [0.0, 0.6)$ and suddenly changes to 7.0 if $x \\in [0.6, 1]$. Figure 8 in the supplementary results visualizes the uncertainty quantification of different surrogate models (The numbers followed by the surrogate name are the sampling sizes, *e.g.*, *MMDGP-nystrom(160/10)* is our method with a sampling size of 160 and sub-sampling size of 10).\n\nWe can clearly observe that i) *MMDGP-nystrom* can comprehensively model the input uncertainty, evidenced by the abrupt change of its posterior distribution at the location $x=0.6$. ii) uGP with a small sampling size fails to quantify the uncertainty but a large sampling size with a higher computation cost help alleviate this problem. iii) A naïve GP model that does not explicitly deal with the input uncertainty cannot be aware of this uncertainty change at all.\n\n### Reference\n[1] Oliveira, Rafael, Lionel Ott, and Fabio Ramos. “Bayesian Optimisation under Uncertain Inputs.” PMLR, 2019.\n\n[2] Sanders, Nicholas D., Richard M. Everson, Jonathan E. Fieldsend, and Alma A. M. Rahat. “Bayesian Search for Robust Optima.” arXiv, 2021.\n"}, {"review_id": "mJTgUJw5dZ", "reviewer": "Reviewer_JQUf", "summary": "This work focuses on the situations where input uncertainty arises and the input values are unobservable, and introduces to measure the distance of uncertain inputs through MMD when training the Gaussian process surrogate. The authors theoretically and empirically demonstrate the effectiveness of the proposed method.", "questions": "Please see the weakness.", "limitations": "No broader societal impacts are provided.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1.\tThe proposed method is sound and MMD can be a metric to measure the distances of uncertain inputs.\n\n2.\tThe experiments are extensive to show the effectiveness of the proposed method.\n", "weaknesses": "1.\tMMD-GP needs to query m times more samples than GP, which can restrict its application since BO is usually applied to the tasks where evaluating a query can be time-consuming. What about the time cost of this work compared to other methods for input uncertainty?\n\n2.\tAfter sampling m times for one query x, we can also get m rewards/performance. I wonder how to deal with the m rewards/performance? Line 139 shows that $D_n=${$(x_i, y_i)|x_i \\sim P_i$}. However, it seems that $y_i$ also should be a distribution.\n\n3.\tOther typos:\n\n3.1\tEq. 2 use $\\xi$ to denote the noise while line 92 utilize $\\epsilon$.\n\n3.2\tLine 151: “the computation and space complexities of the empirical MMD estimator scale quadratically with the sampling size N” Should N be m?\n\n3.3\tWhat does the training and testing samples stand for in line 162? Do you mean the number of observations and the number of queries?\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for the insightful comments.\n\n### 1. Sampling issue\n> MMD-GP needs to query m times more samples than GP, which can restrict its application since BO is usually applied to tasks where evaluating a query can be time-consuming. What about the time cost of this work compared to other methods for input uncertainty?\n\nThanks for pointing this out and sorry for our obscure expression! In our problem setting, we assume our input determinant point $x_i$ is blurred by some noise $\\delta(x_i)$, and then evaluate a time-consuming $f(x_i + \\delta(x_i))$ to get $y_i = f(x_i + \\delta(x_i)) + \\xi_i$ with additional noise $\\xi_i$. However, when doing MMD-based kernel evaluation, we do not have to query the time-consuming function $f$ but only assume $\\delta(x_i)$ can be easily sampled, and thus we have enough samples $\\\\{ \\delta^j(x_i) \\\\}_{1\\le j\\le m}$ for each $1 \\le i \\le n$ to calculate the MMD distance between $x_i + \\delta(x_i)$ and $x_l + \\delta(x_l)$. When calling the evaluation of $f$, we only input one sample $x_i + \\delta(x_i)$, and return one evaluation $y_i$.\n\n### 2. Misleading notations\n> After sampling m times for one query x, we can also get m rewards/performance. I wonder how to deal with the m rewards/performance. Line 139 shows that $D_n = \\\\{ (x_i,y_i) \\vert x_i \\sim P_i\\\\}$. However, it seems that $y_i$\n also should be a distribution.\n\nThanks for pointing this out and sorry for our obscure expression! This question can also be solved by our last answer, as in each query, we only input one $x_i$ to the function $f$ and thus generate one $y_i$ . Also there are some typo in the definition of $D_n$ in our origina paper: $D_n = \\\\{ (\\\\hat x_i,y_i) \\vert \\\\hat x_i \\sim P_{x_i}\\\\}$, here $\\hat x_i$ is one random evaluation of $P_{x_i}$ when we intend to input $x_i$ into the function $f$.\n\n### 3. Typos\n> Other typos\n\nThanks for pointing this out, we have modified all the mentioned typos in the updated paper."}, {"review_id": "KPzYfDKfOA", "reviewer": "Reviewer_UjwV", "summary": "The paper proposes a new method for Bayesian optimization under uncertain inputs, where the distribution of an input can be complex and unknown but can be sampled from. The paper proposes an MMD based kernel between probability distributions and the use of the Nystrom approximation to make the computations tractable. A UCB-type algorithm is used. The paper provides regret bounds accounting for the approximation error due to sampling / Nystrom approximation, and compares the algorithm to some baselines.", "questions": "No additional questions.", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 1, "contribution": 3, "strengths": "1. The motivation of requiring to do BO with uncertain inputs with complex distributions that may not be known in closed form is compelling.\n2. The regret bound incorporating approximation error and the corresponding implications on the sampling size $m$ are useful and interesting.\n\nIn general I believe the work has potential given that the technical issues raised in the Weaknesses section are sufficiently addressed.", "weaknesses": "1. **Validity of kernel**. The paper designs an MMD-based kernel between probability distributions $\\hat k = \\eta(\\text{MMD}(P, Q))$. When designing a new kernel, it is important to prove that it is a valid kernel. However, no such proof is given (or references to previous literature that does something similar). $\\eta$ is said to be a \"normalizing function with range $[0, 1]$\". What conditions does $\\eta$ need to fulfill in order for $\\hat k$ to be a valid kernel?\n\n2. **Empirical evaluation**. The empirical evaluation is unconvincing for the following reasons:\n* The Nystrom approximation is adopted for efficiency. It is important to measure how much using this approximation affects performance beyond simply a qualitative comparison by looking at the posteriors in Figure 3. The empirical evaluation should do an ablation study that includes MMDGP without the Nystrom approximation, at a sampling size that results in the same inference time as that with the Nystrom approximation. What if removing the Nystrom approximation results in better performance with no decrease in efficiency?\n* The closest work seems to be Oliveira et. al. (2019). Why isn't their algorithm one of the baselines?\n* The error bars are simply too large and overlap too much. Any performance improvement could be due to randomness; please run with more trials to decrease the error bars so that the performances can be meaningfully compared.\n\n3. **Writing**. The writing has much to be improved:\n* Line 174 \"One important theoretical guarantee to conduct GP model...\" and in Theorem 1 \"...running Gaussian Process with acquisition function...\" I believe you mean \"Bayesian optimization\" in these contexts instead of \"Gaussian process\". The GP is the model, the algorithm is BO which relies on a GP. The thing being run is the algorithm, not the model.\n* Line 151 \"sampling size $N$\" do you mean $m$?\n* Many grammatical and typo errors: the last word in the title is not capitalized; line 178 \"For $\\hat k$ be radial kernels\" and \"For $\\hat k$ be linear kernel\", should be, for example, \"If $\\hat k$ is a radial kernel\"; Lines 201 and 205 are missing periods; Line 256 \"can be more pronounced impact\" should be \"can have\".\n* [24] and [25] point to the same reference.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for the insightful comments, especially regarding the constructive suggestions on the kernel validity, ablation test, and experiment.\n\n### 1. Validity of kernel\n\n>  When designing a new kernel, it is important to prove that it is a valid kernel. \n>  $\\eta$ is said to be a \"normalizing function with range [0, 1]\". What conditions does $\\eta $ need to fulfill for $\\hat{k}$ to be a valid kernel?\n\nThanks for pointing this out! In our paper, we focus on the case that $\\eta(x) = \\exp{(-\\alpha x)}$ and $\\hat{k} = \\exp(-\\alpha \\text{MMD}(P,Q; k))$. It can be proved to be a valid kernel by using the Theorem 2.2 from [1]: for any type of bi-functions defined on $\\mathcal{P}\\times\\mathcal{P}$ as $\\hat k(P, Q) = \\sum_{i=1}^{\\infty} a_i  \\langle \\psi_P, \\psi_Q \\rangle_k$ with $a_i \\ge 0$, $\\psi$ is the kernel mean map from $\\mathcal{P}$ to $\\mathcal{H}_k$, then $\\hat k(P, Q)$ is a valid kernel; In addition if $a_i > 0$, this kernel is universal over $C(\\mathcal{P})$, given the mean map $\\psi$ is injective and the space of $P, Q$ defined on is compact. In this case, an MMD kernel with RBF mapping $\\hat k(P, Q) = \\exp(-\\alpha \\Vert \\psi_P - \\psi_Q \\Vert_k^2)$ is ensured to be valid and universal.\n\n### 2. Ablation test for Nystrom approximation\n\n> The Nystrom approximation is adopted for efficiency. It is important to measure how much using this approximation affects performance. \n> The empirical evaluation should do an ablation study that includes MMDGP without the Nystrom approximation, at a sampling size that results in the same inference time as that with the Nystrom approximation. \n\nThanks for this constructive suggestion! We have conducted an ablation test for the Nystrom approximation. In this experiment, we employ an RKHS function and set the input uncertainty to follow a beta distribution (see Sec 5.2 in the paper). Several candidates are to study the effect of Nystrom approximation: \n\n+ *MMDGP-nystrom* is our method with Nystrom approximation, in which the sampling size $m=16$ and sub-sampling size $h=9$. Its complexity is $O(MNmh)$, where $M$ and $N$ are the sizes of training and test samples respectively, $m$ is the sampling size for MMD estimation, and $h$ indicates the sub-sampling size during the Nystrom approximation.\n+ *MMDGP-raw-S* does not use the Nystrom approximation but employs an empirical MMD estimator. Due to its $O(MNm^2)$ complexity, we set the sampling size $m=12$ to ensure a similar complexity as the *MMDGP-nystrom*.\n+ *MMDGP-raw-L* also uses an empirical MMD estimator, but with a larger sampling size ($m = 16$).\n+ *GP* utilizes a vanilla GP with a learnable output noise level and optimizes with the upper-confidence-bound acquisition [^a]. \n\nFigure 7a in the supplementary_results.pdf shows that i) with sufficient computation power, the *MMDGP-raw-L* can obtain the best performance by using a large sample size.  However, ii) with the same inference complexity, the *MMDGP-nystrom* performs much better than the *MMDGP-raw-S*, suggesting the Nystrom approximation can significantly improve the efficiency with a mild cost of performance degradation. iii) All the MMDGP-based methods are better than the vanilla GP-UCB.\n\n### 3. Experiments\n\n> The closest work seems to be Oliveira et. al. (2019). Why isn't their algorithm one of the baselines?\n> The error bars are simply too large and overlap too much. please run with more trials to decrease the error bars so that the performances can be meaningfully compared.\n\nThanks for pointing this out! We have cited this paper and included the uGP-UCB method in our baselines and rerun all the experiments with more trials to provide a  statistically-meaningful comparison.\n\nThe uGP-UCB method in [2] employs an integral kernel over probability measures $\\mathcal{P}$: $\\hat{k}(P, P^\\prime):= \\int_\\mathcal{X} \\int_\\mathcal{X} k(x, x^\\prime) dP(x) dP^\\prime(x^\\prime)$. Since it does not mention how to compute this kernel and no public code is available, we compute this integral kernel via sampling, resulting in an inference complexity of $O(MNm^2)$. Here $M$ and $N$ are the sizes of training and test samples respectively, and $m$ is the number of samples used for estimating the integral. \n\nIn the supplementary results, Figure 9a compares our method with the baselines on an RKHS function (Figure 4a in the paper) with a Gaussian input uncertainty. All the robust methods perform well except the vanilla GP-UCB stacks to a local optimum. Also, we notice that the *uGP* method performs slightly better than the others in this case, but we will see later that our method outperforms in more complex distributions and high-dimensional cases.\n\nFigure 9b further compares these models on a double-peak function (Figure 4c in the paper) and a beta distribution. In this case, we observe that the *MMDGP-nystrom* quickly converges to the optimum while the *uGP* hobbles with a larger variance. Also, as mentioned in the paper, the *skl* and *ERBF* fail to locate the robust optimum due to the mismatching of their assumptions and true input uncertainty.\n\nWe also re-evaluate on the robot pushing problem with a multi-modal Gaussian mixture distribution. Each method is tested for 48 trials and the robust regrets are summarized in Figures 9c and 9d. We can observe our algorithm outperforms the others in terms of both optimization efficiency and stability.\n\n### 4. Writing\n\n> The writing has much to be improved...\n\nThanks for pointing them out and sorry for such a rough presentation. We have revised them accordingly and performed careful proofreading for the updated version.\n\n\n[1] Andreas Christmann and Ingo Steinwart. “Universal Kernels on Non-Standard Input Spaces”. NIPS. 2010\n\n[2] Oliveira, Rafael, Lionel Ott, and Fabio Ramos. “Bayesian Optimisation under Uncertain Inputs.” PMLR, 2019. \n\n[^a] For simplicity, all the methods in this work use an acquisition of upper confidence bound.\n "}]

Author rebuttal:
We sincerely thank all the reviewers for their insightful comments and constructive suggestions. In this global rebuttal, we mainly provide the updated version of our theoretical results and a discussion on the limitations, followed by an attached pdf for the supplementary experiments (The other detailed feedback can be found under each reviewer's comments).

### Theoretical Results
For the theoretical part, we summarize the amended theorems/assumptions/descriptions here.

**Assumption 1** For any $\varepsilon > 0$, $P,Q \in \mathcal{P_{\mathcal{X}}}$, we may choose an estimated $\tilde k(P,Q)$ such that the error function $e(P,Q)$ can be upper-bounded by $e_\varepsilon$ with probability at least $1-\varepsilon$, that is, $\mathbb{P}\left(|e(P,Q)| \le e_\varepsilon\right) > 1-\varepsilon.$ 

**Discuss on Assumption 1** Note that this assumption is standard in our case: we may assume $\max_{x \in \mathcal{X}} \Vert \phi \Vert_k \le K $, where $\phi$ is the feature map corresponding to the $k$. Then when using empirical estimator, the error between $\text{MMD}_{\text{empirical}}$ and $\text{MMD}$ is controlled by $4K\sqrt{2\log(6/\varepsilon)m^{-1}}$ with probability at least $1- \varepsilon$ according to Lemma E.1, in [5]. When using the Nystrom estimator, the error has a similar form as the empirical one, and under mild conditions, when $h = O(\sqrt{m}\log(m))$, we get the error of the order $O(m^{-{1/2}}\log(1/\varepsilon))$ with probability at least $1-\varepsilon$. 


**Theorem 1** Let $\delta >0, f\in \mathcal{H}\_k,$ and the corresponding $ \Vert \hat f \Vert_{\hat k}  \le b, \max_{x\in \mathcal{X}} |f(x)| = M$. Suppose the observation noise $\zeta_i = y_i -f(x_i) $ is $\sigma_\zeta$-sub-Gaussian, and thus with high probability $|\zeta_i|< A$ for some $A>0$. Assume that both $k$ and $P_x$ satisfy the conditions for $\Delta f_{P_x}$ to be $\sigma_E$-sub-Gaussian, for a given $\sigma_E > 0$. Then, under Assumption1 with $\varepsilon >0$ and corresponding $e_{\varepsilon}$, setting $\sigma^2 = 1+\frac{2}{n}$, running Gaussian Process with acquisition function 
$\\tilde \\alpha(x|\\mathcal{D}_n) = \\tilde\\mu_n(P_x)  + \\beta_n  \\tilde\\sigma_n(P_x),$ where
$\\beta_n =b + \\sqrt{\\sigma_E^2 + \\sigma\_\\zeta^2}\\sqrt{2(\\hat \\gamma\_n + 1 - \\ln\\delta)  },$
  we have that the uncertain-inputs cumulative regret satisfies:

$\\tilde R_n \\in O ( \sqrt{n\hat\gamma_n(\hat\gamma_n - \ln \delta)}  + n^2\sqrt{(\hat\gamma_n - \ln \delta) e_\varepsilon}  + n^3 e_\varepsilon )$
 with probability at least $1-\delta - n\varepsilon$. Here $\tilde R_n = \sum_{t=1}^n \tilde r_t$, and $\\tilde r\_t = \\max\_{x\\in \\mathcal{X}} \\mathbb{E}\_{P\_x}[f] - \\mathbb{E}\_{P\_{x\_t}} [f]$


**Theorem 2** Suppose $k$ is $r$-th differentiable with bounded derivatives and translation invariant, i.e., $k(x,y) = k(x-y,0)$.  Suppose the input uncertainty is i.i.d., that is, the noised input density satisfies $P_{x_i}(x) = P_0(x - x_i), \forall x_i \in \mathcal{X}$. Then if the space $\mathcal{X}$ is compact in $\mathbb{R}^d$, the maximum information gain $\hat\gamma_n$ satisfies 
\begin{equation}
    \hat\gamma_n = O(n^{(d^2 +d)/(r+d^2+d)} \log(n)).
\end{equation}
Thus, when $r > d(d+1)$, the accumulated regret is sub-linear with respect to $n$ assuming $e_\varepsilon$ is small enough.

### Limitations
Here we provide a discussion of the limitation & further work:
1. The choices of sampling and subsampling size are highly related to the distribution type, dimension, and rank of the kernel matrix. For now, we set these parameters heuristically via experiments, it is worth devising a learnable way to determine their best values automatically.
2. As the power of MMD-based two-sample tests drops polynomially with problem dimension [9], our current method can only alleviate this issue with a quadratic computation cost, which hinders a scale to very high-dimensional problems. We leave the exploration to HD problems for future work.

### Reference

[1] Andreas Christmann and Ingo Steinwart. Universal Kernels on Non-Standard Input Spaces. In:NeurIPS, 2010

[2] Srinivas N, Krause A, Kakade S M, et al. Gaussian process optimization in the bandit setting: No regret and experimental design[J]. arXiv, 2009.

[3] Nogueira, Jose, Ruben Martinez-Cantin, et al. “Unscented Bayesian Optimization for Safe Robot Grasping.” In IROS, 2016.

[4] Tsan Sheng Ng, Yang Sun, et al.,Semiconductor lot allocation using robust optimization, European Journal of Operational Research, Volume 205, Issue 3, 2010.

[5] Chatalic, Antoine et al. Nyström Kernel Mean Embeddings. ICML 2022.

[6] Gabrel, Virginie, Cécile Murat, and et al. “Recent Advances in Robust Optimization: An Overview.” *European Journal of Operational Research* 235, no. 3 (June 2014): 471–83.

[7] Oliveira, Rafael, Lionel Ott, and et al., “Bayesian Optimisation under Uncertain Inputs.” AISTATS, 2019.

[8] Sanders, Nicholas D., Richard M. Everson, et al., Bayesian Search for Robust Optima. arXiv, 2021.

[9] Ramdas, Aaditya et al. “On the Decreasing Power of Kernel and Distance Based Nonparametric Hypothesis Tests in High Dimensions.” AAAI Conference on Artificial Intelligence (2014).

[10] Srinivas N, Krause A, Kakade S M, et al. Gaussian process optimization in the bandit setting: No regret and experimental design[J]. arXiv preprint arXiv:0912.3995, 2009.

[11] Christmann A, Steinwart I. Universal kernels on non-standard input spaces[J]. Advances in neural information processing systems, 2010, 23.

Meta-review: There were initially some non-minor concerns, particularly around the kernel validity and restrictiveness of the main assumption.  The concerns mostly seem to have been alleviated following the rebuttal.  However, the paper might be approaching the boundary of what amount of required editing is acceptable for camera-ready.  The authors should very carefully incorporate the required changes, including some of the writing clarity as mentioned by Reviewer UjwV.

In addition to the reviewer comments, I would mention some weaknesses in the main theoretical result.  The weakness in the sample size (e.g., n^6 scaling) is already discussed following the theorem, which is appreciated.  Another weakness is having gamma_n * sqrt(n) dependence – in the vanilla setting, several algorithms are now known getting the much better sqrt(n * gamma_n) dependence instead: 
- Finite-time analysis of kernelised contextual bandits
- A domain-shrinking based Bayesian optimization algorithm with order-optimal regret performance
- High-dimensional experimental design and kernel bandits
- Gaussian process bandit optimization with few batches
Since the authors are tackling a more challenging setting and have a relatively practical focus, I wouldn’t reject based on this limitation, but I do suggest discussing it.
Decision: Accept (poster)