Maximally Robust Satisficing Bayesian Optimization
Samuli Kinnunen1 1
Petrus Mikkola1
arXiv:2607.13652v1 [cs.LG] 15 Jul 2026
Arto Klami1
Department of Computer Science, University of Helsinki, Finland 2 ASM International N.V.
Abstract Many design tasks can be cast as black-box function optimization, enabling use of Bayesian optimization to find an ideal design with minimal number of trials. However, often we do not actually need the optimum but instead a sufficiently good solution is enough, for instance a material that is durable enough for its intended use. In most cases there are multiple satisfactory solutions, forming a superlevel set of the function, raising a key question of which one to prefer. We answer this by explaining why robustness to input perturbations that may occur when the solution is deployed is a good criterion and by introduce a Bayesian optimization method that efficiently finds satisficing solutions that are robust to maximally large perturbations. In contrast to previous works, we assume the inputs can be accurately controlled during optimization, but will be perturbed after the deployment.
1
Antti Niskanen2
INTRODUCTION
Bayesian optimization (BO) is a sample-efficient framework for optimizing expensive black-box functions, applied in diverse domains such as hyperparameter tuning Cho et al. [2020], materials discovery Chitturi et al. [2024], chemical synthesis Shields et al. [2021], and robotic control LechuzSierra et al. [2024]. For the standard problem of finding the global optimum, x∗ = arg maxx f (x), there are highly polished tools. However, for many design tasks where the evaluations are extremely costly, e.g. due to involving production and testing of a new material, even a sample-efficient search may still be difficult in practice. Often we do not, however, need the global optimum x∗ , but already a sufficiently good solution that meets a given minimum quality, f (x∗ ) ≥ t, is enough. For instance, a subcontractor may be asked to produce a component that
achieves a certain level of durability, reducing the difficulty of the problem. We call such a solution satisficing, following the nomenclature introduced by Herbert Simon [Simon, 1955] and subsequently used e.g. in operations research and decision-making [Reverdy et al., 2016, Long et al., 2023, Saday et al., 2025]. The set of all satisficing solutions forms a superlevel set St = {x : f (x) ≥ t} of the function. In principle any solution within St would work, but outside of idealized conditions they are not alike. One would intuitively steer away from solutions close to the threshold of f (x) = t and prefer a solution with a higher value, but we need a more formal criterion. We argue that robustness to input noise or variation is a good one: A solution is good if f (x∗ ) remains satisficing when x∗ is perturbed. We draw motivation from applications where the inputs can be accurately controlled in the design (optimization) phase but are unreliable after deployment. As a concrete example, consider a cake recipe designed by professionals that measure the ingredients accurately and use consistent kitchen equipment. When commoners bake the cake following the recipe, they will not measure their tablespoons exactly and their ovens do not maintain the temperature perfectly. What ultimately matters is how good their cake tastes, and hence a recipe robust to this kind of deviations is a better one. The same principle is seen in many scientific applications: manual procedures and environmental factors like humidity introduce variability in laboratory chemistry [Aldeghi et al., 2021], sim-to-real transfer [Doersch and Zisserman, 2019] is a whole subfield studying specifically how to train reliable real-world models using clean simulation environments, and in design and manufacturing [Daulton et al., 2022] fabrication tolerances constraint the realizability of designs. Often the solution that is optimal under noiseless inputs can perform poorly or even lead to unsafe behavior under perturbations Oliveira et al. [2019]. To address this need, we introduce a new problem: Given a threshold t that determines the minimum acceptable value
(a) Maximally Robust Satisficing
0.8
x2
0.6 0.4 0.2 0.25
0.50 x1
0.75
10−1
10−2
10−3
1.00
(c) MRSBO acquisition values
MRSBO AdveRS2
0
50 Iteration
100
1.0
(d) AdveRS2 acquisition values
0.8
0.8
0.6
0.6
1.0
0.75 0.60
x2
x2
0.90
0.4
0.4
0.2
0.2
0.0 0.00
0.0 0.00
0.45 0.30
0.25
0.50 x1
0.75
1.00
0.15 0.25
0.50 x1
0.75
1.00
Normalized AF value
0.0 0.00
2020, Shapiro et al., 2023] accounting also for uncertainty in the perturbation distribution. These methods ensure the expected performance to be good, but there are no guarantees for the solution to satisfy f (x) > t for any x.
(b) Optimization performance Log Inference Robustness regret
1.0
0.00
Figure 1: (a) The maximally robust satisficing solution (red) is a point that can be modified the most while still staying within the superlevel set defined by the threshold that needs to be satisfied (blue). Note that the global maximum (black cross) can be arbitrarily far. (b) Our MRSBO method drastically outperforms the best previous method (AdveRS2) adapted for our setup; it also optimizes for the robustness radius, but assumes input perturbations already during optimization. (c-d) Acquisition functions for the two methods during late stages of optimization. MRSBO explores the parts of the level set that are informative of the robustness radius, whereas AdveRS2 wastes evaluations for learning the target within the superlevel set.
for the function, find a solution x∗ that is robust to maximally large input deviations. During optimization we assume f (x) can be evaluated for arbitrary x, but the solution is eventually deployed in an environment where the inputs are perturbed by an unknown noise. That is, instead of f (x∗ ) we observe f (x∗ + δ). The task is illustrated in Figure 1, showing also how our method maximizing the information gain about the maximal robustness radius solves the problem efficiently by primarily evaluating the function around selected parts of the level set, not wasting evaluations to learn the function accurately within the superlevel set. While our problem setup is new, it is closely related to the concept of Robust Bayesian Optimization (RBO), generally referring to optimization methods finding solutions that remain effective under input perturbations. Previous RBO methods can be broadly categorized into two types. Expectation Robustness [Yang et al., 2023, Fröhlich et al., 2020] optimizes expected performance under a known perturba tion distribution p(δ), x∗ = arg maxx Eδ∼p(δ) f (x + δ) , with Distributionally Robust Optimization [Kirschner et al.,
Worst-case Robustness [Bertsimas et al., 2010, Bogunovic et al., 2018] is closer to our setup, searching for a solution that performs optimally under the most adverse perturbation as x∗ = arg maxx minδ f (x + δ). Most methods in this family, like StableOpt by Bogunovic et al. [2018], directly optimize for the worst-scare performance, but recently Saday et al. [2025] connected the worst-case robustness to satisficing solutions. They also maximize the magnitude of the perturbations that can be allowed while ensuring the solution remains above the threshold t. However, there is a fundamental difference in the learning setup: Similar to other worst-case robustness methods, Saday et al. [2025] assume already the optimization is done in an adversarial environment. This is in stark contrast to our setup where the inputs are clean during optimization. Their method can be used to solve our problem by artificially inducing perturbations during optimization, but this requires assuming a maximum level of perturbation. Figure 1 shows how this approach in practice fails even when given oracle information about the maximum perturbation, with substantially worse convergence, poor final regret, and wasteful evaluations. Contributions. Our main contributions are: (1) We formalize Maximally Robust Satisficing Bayesian Optimization problem setup for learning solutions that remain sufficiently good under maximal input perturbation. (2) We provide practical algorithm MRSBO for solving it, based on the information gain about the maximally robust solution with an efficient construction for the acquisition function. (3) We demonstrate that MRSBO vastly outperforms robust and non-robust alternatives on a variety of synthetic and real-world robust optimization problems.
2
BACKGROUND
2.1
BAYESIAN OPTIMIZATION FRAMEWORK
Bayesian Optimization, in its standard form, is a global optimization method for solving x∗ = arg max f (x), x∈X
(1)
where X ⊆ Rd is a d-dimensional design space and f : X → R is an unknown function that is expensive to evaluate. To optimize the function, we can only evaluate inputs x to produce an evaluation y = f (x). BO algorithm learns cheap surrogate model, typically a Gaussian process, that approximates the real function and selects the next evaluation by maximizing an acquisition function to identify the evaluation that maximally benefits the optimization process. Surrogate model approximates the true objective
function f . A Gaussian process (GP) Williams and Rasmussen [1995] is the most commonly used surrogate model, providing a distribution over function values f (x) characterized by a mean function µ(x) and a covariance (kernel) function k(x, x′ ). Given observations Dt = {(xi , yi )}ti=1 , the posterior predictive distribution at a test point x is Gaussian with mean µt (x) = kt (x)⊤ (K + σϵ2 I)−1 y and variance σt2 (x) = k(x, x) − kt (x)⊤ (K + σϵ2 I)−1 kt (x),, where kt (x) = [k(xi , x)]ti=1 , K = [k(xi , xj )]ti,j=1 is the kernel matrix of the training inputs, y = [y1 , . . . , yt ]⊤ , and σϵ2 denotes the observation noise variance.
Algorithm 1: Maximally Robust Satisficing Bayesian Optimization (MRSBO) Input: Design space X , GP prior (µ(x), k(x, x′ )), acquisition function α(x), budget T Output: Recommended point x̃∗ 1 Initialize data D0 2 for t = 1, . . . , T do 3 {(x∗,(s) , r∗,(s) )}ns=1 ← S AMPLE C ONTINUOUS 4 xt ← arg max α x′ , {(x∗,(s) , r∗,(s) )}ns=1 , Dt−1 x′ ∈X
yt ← f (xt ) + ϵt 6 D t ← Dt−1 ∪ {(xt , yt )} Acquisition function (AF) α(x) scores the value of eval7 end uating an input point x. At each Bayesian optimization ∗ iteration, the next query point is selected by solving 8 x̃ ← arg maxx∈X ρ̃(x) xt+1 = arg maxx∈X α(x). Popular single-objective acquisition functions include Upper Confidence Bound (UCB), the motivational example cases, whereas they assume the Expected Improvement (EI) Jones et al. [1998] and Maxinputs are perturbed also during optimization (and that the value Entropy Search (MES) Wang and Jegelka [2017]. perturbed inputs can, perhaps unrealistically, be observed for learning the proxy).
3
5
PROBLEM DEFINITION
We consider design tasks where a solution x∗ is determined in controlled settings, but after fixing the design it is deployed in an environment where every input x ∈ X may experience an unknown perturbation δ whose magnitude is bounded by an unknown robustness radius r > 0. For any r, the admissible perturbations are collected in Br = δ ∈ Rd : ∥ δ ∥ < r , and we define as satisficing the x ∈ X that remain satisfactory under all admissible perturbations: the worst–case function value stays above a given threshold t ∈ R as inf f (x + δ) ≥ t.
δ∈Br
To find the maximally robust solution, we define the objective function ρ(x) as the maximum robustness radius r such that all perturbed inputs yield satisficing values: ρ(x) = sup r > 0 : inf f (x + δ) ≥ t and x + δ ∈ X . δ∈Br
(2) This immediately gives the final learning task that provides both the optimal solution x∗ and the corresponding robustness radius r∗ as x∗ = arg max ρ(x),
(3)
r = ρ(x∗ ).
(4)
x∈X ∗
Our objective is identical to the one recently proposed by Saday et al. [2025], but the learning setup is different. We assume exact evaluations f (x) during optimization, matching
4
METHOD
This section introduces the MRSBO method for solving the above problem, summarized in Algorithm 1. It follows a standard Bayesian optimization protocol: At each iteration, candidate maximizers are sampled (Line 3), and the acquisition function is optimized to select the next evaluation point (Line 4). The proxy is updated with the new observation, and we continue until exhausting an evaluation budget T . We here design an information-theoretic acquisition function and explain how it can be computed efficiently. The acquisition maximization step can be carried out using standard gradient-based optimization methods, e.g., as discussed by Wilson et al. [2018]; we leave description of this standard procedure in Appendix B.1. 4.1
INFORMATION GAIN ABOUT (x∗ , r∗ )
We propose an acquisition strategy that maximizes the information gain about the maximum robustness radius r∗ and its corresponding robust maximizer x∗ . Following the symmetric mutual information formulation [Hernández-Lobato et al., 2014, Wang and Jegelka, 2017], this approach avoids direct entropy computation of the maximizer and robustness radius, which are difficult quantities even to sample. Formally, the acquisition function for candidate x is α(x) = I({y, x}; {x∗ , r∗ }) = I({x∗ , r∗ }; {y, x})
= H(p(y | x, Dt )) − Ep(x∗ ,r∗ |Dt ) H(p(y | x, Dt , x∗ , r∗ ))
(5)
The first term of Eq. 5, H(p(y | x, Dt )), is the entropy of a Gaussian random variable and admits the analytic form 1 H p(y | x, Dt ) = log 2πeσt2 , 2
Z(y) = (6)
where σt2 denotes the posterior predictive variance at x. The expectation in the second term of Eq. 5 is taken with respect to the joint posterior of the robust maximizer and robustness radius, and is estimated using Monte Carlo sampling. We next discuss the approximation of the conditional entropy, H(p(y | x, Dt , x∗ , r∗ )), starting with a noiseless case for communicating the core idea even though in practice we always assume noisy observations. Noiseless conditional entropy. In the absence of observation noise (σϵ2 = 0), the computation can be separated into two cases based on the distance d(x, x∗ ). If the candidate x is within the robustness radius (d(x, x∗ ) ≤ r∗ ) then, by definition of robustness, we are guaranteed to have f (x) ≥ t. When d(x, x∗ ) > r∗ , we do not have any information about the values of f (x) unless making additional assumptions about the function. f (x) can still be arbitrarily high since the superlevel set can take any form. Hence, the conditional entropy can be expressed as: H(p(y | x, Di , x∗ , r∗ )) ( √ ∗ ∗ log( 2πe σZ) + αφ(α) 2Z , if d(x, x ) ≤ r , ≈ 1 2 otherwise, 2 log(2πeσt ),
(7)
where α = (a − µ)/σt is the standardized lower truncation point, with a being the lower bound. Z and φ are defined as 1 1 2 Z = 1 − Φ(α) φ(ξ) = √ exp − ξ , 2 2π where Z(x) = 1−Φ α can be computed using the standard normal cumulative distribution function Φ. For smooth functions we could attempt lower-bounding the function also outside of the robustness radius: For x just outside the robustness radius, we necessarily have f (x) close to t. Exploiting this in practice, however, would require additional assumptions about the smoothness of the function, which are technically challenging to specify correctly [Lederer et al., 2019] and yield little empirical benefit, as shown in Appendix A.3. We therefore recommend using the above conditional entropy approximation. Noisy conditional entropy. With observation noise, the guarantee that f (x) > t for d(x, x∗ ) ≤ r∗ no longer holds, but the distance still provides information. For computing the entropies, we adapt the approach of Takeno et al. [2020] proposed for noisy truncation in multi-fidelity BO context. The conditional predictive density can be expressed as p(y | x, Dt , f > t) = Z(y) p(y | x, Dt ),
where p(y | x, Dt ) is the posterior predictive density and
(8)
p(f > t | y, x, Dt ) . p(f > t | x, Dt )
The denominator corresponds to a Gaussian tail probability, while the numerator is computed from the conditional distribution f | y, x, Dt ∼ N (u, s2 ), with u=µ+
σf2 (y − µ), 2 σf + σϵ2
s2 = σf2 −
σf4 . 2 σf + σϵ2
Finally, the conditional entropy H(p(y | x, Dt , x∗ , r∗ )) is approximated via Monte Carlo samples {yi }N i=1 as N
−
1 X log p(yi | x, Dt , f > t) if d(x, x∗ ) ≤ r∗ , N i=1 1 log(2πeσt2 ), 2
4.2
otherwise.
SAMPLING FROM p(x∗ , r∗ |Dt )
Evaluation of the acquisition function requires computing expectations with respect to the posterior p(x∗ , r∗ | Dt ). Since this distribution is generally intractable, we approximate the expectation using Monte Carlo integration. This reduces the problem to efficiently sampling pairs (x∗ , r∗ ). Following Hernández-Lobato et al. [2014, 2016], we approximate samples from the GP posterior using Random Fourier Features (RFF). This yields differentiable sample paths for the sampling procedure in Algorithm 2. RFF representation By the Bochner’s theorem, the spectral density of the GP kernel k exists and we can approximate k as an empirical expectation of the feature inner product k(x, x′ ) ≈ ϕT ϕ. For the radial basis function (RBF) kernel with the length scales ℓ and the output scale σf , ϕ is an m-dimensional random feature function defined as r 2σf ϕ(x) = cos(W x + b), m where the rows of W ∈ Rm×d are sampled independently as w ∼ N (0, I/ℓ) and the phases are sampled as b ∼ U(0, 2π). We approximate the GP prior using a Bayesian linear model f˜(x) = ϕ(x)⊤ θ, where θ ∼ N (0, I) is a Gaussian-distributed weight vector. The posterior of θ|Dt is Gaussian, θ ∼ N A−1 Φ⊤ y, σ 2 A−1 , where Φ = [ϕ(x1 ), ..., ϕ(xt )]⊤ is the feature matrix and A = Φ⊤ Φ + σ 2 I. Robust maxima To generate a single approximate posterior sample fs , we draw both the random feature parameters (ws , bs ) and the corresponding posterior weight vector θs
0
0.8
x2
0.6 0.4
0.8
2
0.6
4
0.2
5
0.25
0.50
x1
0.75
6
1.00
2 3
0.4
4
0.2 0.0 0.00
5
0.25
0
0.4
0.8
2
0.6
4
0.2
5
x1
0.75
6
1.00
0.75
1.00
6
approx. level set points x* True level set
1.0
1 3
0.50
x1
1 0
1 2
x2
x2
0.6
1
Normalized function value
approx. level set points robust maximizer candidates True level set
0.25
0.50
(b)
0.8
0.0 0.00
0
1
(a) 1.0
1
3
0.4
4
0.2 0.0 0.00
Normalized function value
0.0 0.00
1 3
approx. level set points True level set
1.0
Normalized function value
1
x2
approx. level set points True level set
Normalized function value
1.0
5
0.25
(c)
0.50
x1
0.75
1.00
Algorithm 2: S AMPLE C ONTINUOUS Input: Surrogate GP model GP, number of level set samples S, number of basis functions m, tolerance τ , threshold t Output: {(x∗s , rs∗ )}Ss=1 for s = 1, . . . , S do Construct differentiable sample function f˜s (x) = ϕs (x)⊤ θs Approximate level set samples Cs ← L EVEL S ETA PPROX(f˜s , t, τ ) Initialize robust maximizer candidates Ss ← I NITIALIZE M AXIMIZATION(f˜s , Cs ) Compute robust maximizer x∗s = arg maxx∈Ss − β1 LogSumExp(−d(x, Cs )) Compute robustness radius rs∗ = minc∈Cs d(x∗s , c) end
6
(d)
Figure 2: Sampling of {(x∗s , rs∗ )}Ss=1 . For each s, spacefilling initial candidates (a) are refined using Gauss-Newton to find the level set (b). A robust maximizer (d) is found by gradient-based optimization of LogSumExp distance to the level set samples, multi-started from a diverse subset of the level set samples pushed to the superlevel set (c). according to the generative process above. This yields a differentiable sample function f˜s (x) = Φs (x)⊤ θs . This enables efficient Algorithm 2 for finding the robust maximizer (x∗s , rs∗ ) for each f˜s . It is visually illustrated in Figure 2, consisting of four key steps explained next. First, we construct an estimate of the (inflated) level set Cs = {x ∈ X : |fs (x) − t| < τ }, for a tolerance τ . The set Cs is approximated by first sampling k space filling samples (x1 , . . . , xk ) (Fig. 2 (a)), and using the Gauss-Newton algorithm starting from each of these initial samples to minimize the squared residual (fs (xi ) − t)2 until the tolerance is met; see Algorithm 14 in Appendix. This gives k samples covering the level set (Fig. 2 (b)). The robust maximizer is then obtained by solving " # X 1 ∗ xs = arg max − log exp −β d(x, c) , x∈Ss β
(9)
c∈Cs
where the differentiable objective corresponds to a softmin approximation of the distance to the level set. We first
generate a set of n ≪ k superlevel set points Ss = {x ∈ X : f˜s (x) ≥ t} by running a few iterations of gradient ascent for fs starting from the level set points Cs (Fig. 2 (c)). The set of n samples are selected to be maximally diverse via a max–min distance criterion; see Algorithm 19 in Appendix. In the final step, x∗s is found as the maximizer of the objective 12 (i.e. LogSumExp(−d(x, Cs ))) by running multistart gradient ascent starting from the points Ss and picking the best one x∗s . This approximates the robust maximizer of fs (Fig. 2 (d)), and the associated robustness radius is rs∗ = minc∈Cs d(x∗s , c).
Computational cost Even though this algorithm requires an explicit construction of the level set and involves two separate optimization stages, we can use gradient-based methods for all parts and the sampling can be done in parallel for all of the MC samples. For one posterior sample, sampling the approximate level set scales linearly with the number of approximate level set points k, basis functions m, and the dimensionality d, with complexity O(kmd). The computation of the robust maximizer is dominated by evaluating the distances between the n maximizer candidates and the k level set points, which scales as O(knd). To accelerate this phase when k is large, we restrict the LogSumExp distance computation in each iteration to the 100 nearest points, since distant points have a negligible contribution. To maintain accuracy, the set of nearest points is re-estimated every 10 iterations. In practice the computational cost is not an issue in typical applications where the function evaluation cost dominates. For example, in a three-dimensional problem the cost for selecting the next candidate was on average 2.4 seconds, only three times the cost of standard BO acquisition (2).
4.3
ROBUSTNESS REGRET
We propose two robustness-based regret metrics for evaluation. The expected robustness regret at iteration i is Ri = |ρ(x∗ ) − Ep(x̃∗ |Dt ) [ρ(x̃∗i )]|,
(10)
where x̃∗i denotes the estimated maximizer of Eq. (3) at iteration i. The expectation is taken over GP sample paths, and the true objective value is evaluated for each sample path. This can be computed for any method, but requires GP approximations similar to what we use in sampling. As a more tractable alternative, we define the inference robustness regret, relying only on the GP posterior mean: R̃i = |ρ(x∗ ) − ρ(x̃∗i )|,
(11)
where ρ̃i is the posterior mean of ρ at iteration i, and x̃∗i is the estimate of x∗ at iteration i, that is x̃∗i = arg maxx∈X ρ̃i (x). This surrogate regret is inexpensive to compute and enables efficient evaluation.
5
RELATED WORK
As explained in Introduction, our setup is close to the broad literature on robust BO, especially to worst-case robust optimization, but differs by assuming perturbation-free inputs during optimization. Ignoring this distinction for the moment, the algorithm itself can be contrasted with prior work. AdveRS-2 [Saday et al., 2025] and StableOpt [Bogunovic et al., 2018] both rely on the Upper Confidence Bound (UCB) acquisition function, which scales poorly in higher dimensions. StableOpt requires solving an inner optimization problem inf δ∈Br f (x + δ) ≥ t for each candidate point during acquisition optimization. AdveRS-2, in turn, discretizes the design space to identify the point that maximizes the distance to the unsatisficing sub-level set, with exponential cost. We avoid both limitations. The sampling of the maximizers (x∗ , r∗ ) is independent of the candidate and hence we can first sample the maximizers and only then optimize the acquisition function efficiently using gradients. Information-theoretic acquisition functions have been used in Expectation Robustness setting Fröhlich et al. [2020] and in a specific worst-case contextual variable robustness usecase Weichert et al. [2024], and have long-standing roots in standard Bayesian optimization Hernández-Lobato et al. [2014]. Our formulation is closely related to the symmetric mutual-information acquisition functions [HernándezLobato et al., 2014, Wang and Jegelka, 2017]. It also aligns with approaches that jointly consider information about both the maximum value and its location Hvarfner et al. [2022], which has been shown to yield more informative acquisitions. Unlike their method, however, ours does not require conditioning the surrogate model on a fantasized maximum, avoiding an additional computational step.
Satisficing BO is intrinsically related to estimation of the level set where f (x) = t. Finding the level set of the true (black-box) function is an extremely hard problem, especially in higher dimensions: Bachoc et al. [2021] showed that the sample complexity depends on the packing number of the (inflated) level set, while Gotovos et al. [2013] introduced a GP-UCB-style algorithm with probabilistic convergence guarantees under the assumption that f is drawn from a GP on a discrete domain. While we need to estimate the level set within the sampling algorithm, it is critically important to note that we do not need the level set of the true f (x) but only the level sets of GP samples fs (x). These samples are smooth, we are not constrained in sample evaluations, and we can use gradients. This makes learning the level set relatively easy for reasonably smooth targets. Moreover, we do not need perfect coverage of the level set for finding the robust maximizer, as will be empirically demonstrated in Section 6.2. One relevant line of research relates to the nature of the inputs. Fröhlich et al. [2020] and Daulton et al. [2022] consider robustness specifically for controllable inputs, whereas Weichert et al. [2024], Toscano-Palmerin and Frazier [2022] explicitly assume there are two distinct sets of inputs: design parameters to be optimized and uncontrollable environmental parameters with respect to which the solution should remain robust. Our presentation assumes all inputs will be noisy after deployment, but we can easily extend the approach to accommodate this setting by introducing a binary mask m ∈ {0, 1}d that specifies which coordinates will be perturbed. This changes the admissible permutation set to Br (m) = δ ∈ Rd : ∥m ⊙ δ∥ < r and (1 − m) ⊙ δ = 0 but requires no other modifications.
6
EXPERIMENTS
We compare MRSBO with two worst-case robust optimization algorithms and one standard maximization method, explained below, using identical GP surrogate with RBF kernel implemented using BoTorch [Balandat et al., 2020] for all methods. The proxy hyperparameters are optimized after each evaluation. Each method is evaluated over 30 independent random initializations, using identical initial designs across all baselines, and the performance is evaluated using the inference robustness regret of Eq. 11, computed with respect to the maximally robust estimate x̃∗ obtained based on the surrogate model. See Appendix B for complete experimental setup. The code for reproducing all experiments is available at https://github.com/KinnunenSamuli/MRSBO/. Max-value Entropy Search (MES) Wang and Jegelka [2017] is included as a strong baseline to show how standard optimization algorithms looking for arg maxx f (x) ignoring robustness of the satisfiability threshold t work in our setting. MES selects points that maximize the mutual
information about the unknown maximum value y ∗ using αM ES (x) = H(y | x, Dt ) − Ep(y∗ |Dt ) [H(y | x, Dt , y ∗ )] . StableOpt Bogunovic et al. [2018] optimizes worst-case performance under adversarial perturbations with known maximum level ϵ with UCB acquisition αSO (x) = min UCB(x + δ), δ∈∆(x)
where ∆(x) = {x′ − x : x′ ∈ X , d(x, x′ ) < ϵ}. AdveRS-2 Saday et al. [2025] maximizes our objective but in an adversarial environment. They maximize the optimistic robustness radius under worst-case input perturbations, subject to the constraint f (x) > t, with acquisition function αA-RS-2 (x) = d(x, x′ )
s.t.
UCB(x + δ) ≥ t.
For both robust baselines, we follow the authors’ recommendations and use β = 2 to control the UCB exploration level. To apply these methods in our setup where the optimization is carried out without adversarial perturbations, we simulate the perturbations by drawing them uniformly within a ball of ϵ-radius. In absence of prior knowledge the ϵ parameter could not be determined prior to running the algorithm in realistic settings, but we consider an idealized setup where we have access to the maximally robust solution r∗ and set ϵ = r∗ . This can be considered an optimal choice, but we will also illustrate the effect ϵ has on the methods. 6.1
BENCHMARKING
Illustration Figure 3 illustrates the task and the different ways the algorithms solve it. The target is a two-dimensional Gaussian mixture with multiple modes, showcasing a typical scenario where the global maximizer of the function (at [0.95, 0.95]) is vastly different from the true robust maximizer at x∗ = [0.51, 0.51]. MES obviously finds the global maximizer and hence has substantially worse final regret due to the small robustness radius for that solution, and StableOpt that ignores the satisficing constraint converges to the widest mode – it is indeed a robust solution, but not optimal for our objective. Both MRSBO and AdveRS-2 converge (approximately) to the robust satiscifing maximizer, but MRSBO achieves it an order of magnitude faster. To provide insight on how the methods work, we show histograms of the first 20 evaluations, accumulated over the 30 repeated runs. MRSBO spends most of the evaluations exploring the parts of the level set that are informative of the robustness radius, and interestingly often does not evaluate the function at all around the eventual solution. This is in contrast with AdveRS2 that spends most acquisitions learning the function within the superlevel set. MES naturally explores the modes to find the global optimum. StableOpt explores the perimeter of the largest mode, demonstrating somewhat similar evaluation preference as ours.
Table 1: Optimization performance (mean and standard error of R̃) measured with the budget T = 20d evaluations, with bold indicating the best method for each task. See Fig. 11 in Appendix for regret curves and exact threshold values. Method
MRSBO
MES
A DVE RS2
S TABLE O PT
low threshold Branin 0.14 ± 0.04 0.44 ± 0.08 0.09 ± 0.02 0.19 ± 0.04 GM 2D 0.00 ± 0.00 0.03 ± 0.00 0.01 ± 0.00 0.04 ± 0.01 GM 3D 0.00 ± 0.00 0.05 ± 0.01 0.02 ± 0.00 0.04 ± 0.01 RP 3D 0.13 ± 0.03 0.32 ± 0.06 0.14 ± 0.03 0.27 ± 0.07 Hart. 3D 0.00 ± 0.00 0.04 ± 0.00 0.01 ± 0.00 0.05 ± 0.01 high threshold Branin 0.06 ± 0.01 0.09 ± 0.02 0.05 ± 0.03 0.13 ± 0.03 GM 2D 0.00 ± 0.00 0.01 ± 0.00 0.00 ± 0.00 0.03 ± 0.01 GM 3D 0.00 ± 0.00 0.03 ± 0.00 0.00 ± 0.00 0.04 ± 0.01 RP 3D 0.13 ± 0.03 0.26 ± 0.04 0.11 ± 0.03 0.22 ± 0.04 Hart. 3D 0.00 ± 0.00 0.01 ± 0.00 0.01 ± 0.00 0.02 ± 0.00
Low-dimensional targets Optimization performance is evaluated on four synthetic benchmarks (Hartmann, Branin, and two Gaussian mixtures) and one real-world task (Robot Pushing Kaelbling and Lozano-Pérez [2017], Wang and Jegelka [2017]). For each function we consider three satisticing thresholds t and the optimization is initialized with 5 random evaluations. Table 1 shows the regrets for two choices of the threshold t for each data, with evaluation budget T = 20d; see Appendix C for full regret curves and a third (intermediate) threshold. MRSBO is in general the best solution, reaching zero regret in six scenarios, but AdveRS2 is slightly better in Branin, as well as Robot Pushing for the higher threshold. The quality of MES depends purely on how good the global maximizer happens to be; it achieves good regret when it coincides with the robust maximizer. Higher-dimensional targets Previous robust BO methods are restricted to targets of very low d due to their inherent max-min optimization structure, often constructing a grid for finding the maximum of α(x) and scaling exponentially with d. Our acquisition function is differentiable and we can find the maximum with gradient-based methods as in Wilson et al. [2018], and hence we avoid this scaling. Constructing the level set for sampling the robust maxima for individual samples intuitively gets harder in higher dimensions, but there is no clear computational bottleneck. Figure 4 shows regret curves for a six-dimensional target function (Hartmann 6D) for two thresholds t for MRSBO and MES; StableOpt and AdveRS-2 would no longer scale to this. For sufficiently low threshold MRSBO retains strong performance, converging fast to a good solution. When t approaches the maximum, here y ∗ = 3.3, the superlevel set shrinks and all satisficing solutions are here around
1.0
MRSBO
AdveRS2 ( = 0.15)
AdveRS2 ( = 0.25)
AdveRS2 ( = 0.35)
MES
StableOpt ( = 0.15)
StableOpt ( = 0.25)
StableOpt ( = 0.35)
True solution
0.8 0.6 0.4 0.2
1.0 0.8 0.6 0.4 0.2 0.0 0.0
0.5
1.0
0.0
0.5
1.0
MRSBO AdveRS2 ( = 0.25)
0.0
0.5
1.0
AdveRS2 ( = 0.35) MES
0.0
0.5
StableOpt ( = 0.15)
1.0
Log Inference Robustness regret
0.0
Optimization performance 10−1
10−2
10−3 0
50
100
StableOpt ( = 0.35) Iteration
StableOpt ( = 0.25)
0
80 160 240 320 400 480 560 Iteration
MRSBO
100
10−2 0
80 160 240 320 400 480 560 Iteration
the global maximizer. MES is then more efficient. See Appendix C for results on other problems of 4-6 dimensions. 6.2
80 60 40
LEVEL SET QUALITY
The quality of the robust maximizer and robustness radius samples {x∗ , r∗ }Ss=1 depends on well the level set samples Cs approximate the true level set. Since Algorithm 14 rejects samples that fail to converge, the accuracy is determined by how well Cs covers the full level set. For low-dimensional targets we can evaluate the coverage using a dense grid over the space: We identify the cells that intersect the true level set and then measure which fraction of these cells is covered, in the sense that at least one sample in Cs is in its immediate proximity, up to the grid resolution. Figure 5 illustrates the level set coverage and its effect on the robust maximizer identification for a bi-modal Gaussian mixture of varying dimensionality, showing both the level set coverage and the fraction of runs that converge to the correct robust maximizer as a function of k, the number of level
d=2 d=3 d=4
20 0
MES
Figure 4: Performance on 6-dimensional Hartmann target.
Level Set Coverage Fraction of grid points within 2h of s
102
103 kbase (k = kbase × d)
Success rate (%)
10−2
Hartmann 6D, t=2.5 10−1
Coverage (%)
Hartmann 6D, t=1 10−1
Log Inference Robustness Regret
Log Inference Robustness Regret
Figure 3: The first four columns show histograms of the first 20 evaluations aggregated over 30 runs for each acquisition strategy, with the right column showing the solution and regrets. MRSBO directly focuses on evaluations informative of the robustness radius and MES explores the modes to find the global optimum. StableOpt and AdveRS2 correctly explore roughly the area where the robust maximizer is, but AdveRS2 is highly inefficient and StableOpt that ignores the satisficing threshold converges to a wrong solution. For both methods, higher adversarial perturbation ϵ increases exploration.
100 80 60 40 20 0
x * Recovery % seeds: x * x * < 0.05
d=2 d=4 d=6 d=8 d = 10 103 kbase (k = kbase × d)
Figure 5: Left: Level set coverage. Fraction of the true level set of a bi-modal Gaussian mixture covered by the approximation set Cs . Since the coverage is estimated using a dense grid, the coverage is shown only for d ≤ 4. Right: Robust maximizer identification. Fraction of robust maximizer estimates x∗s that converge to the analytic true solution of the bi-modal Gaussian mixture.
set samples. The coverage improves with k and degrades with d, as expected. Importantly, partial coverage is often enough to find the correct maximizer. For example, with d = 4 any k is sufficient for finding the correct maximizer even though the coverage is always below 80%, and even for d = 10 (for which coverage cannot even be estimated) we can reliably find the true robust maximizer with large k. See Appendix A.2 for analysis of the coverage for other targets. 6.3
PARAMETER SENSITIVITY
Sampling of robust maximizers (x∗ , r∗ ) depends on a few technical hyperparameters. Figure 6 shows how the robust
Log inference Robustness regret
GM 3D
Hartmann 3D
Acknowledgements
10 2 10 3 1000 2000 3000 # level set samples
0
50 100 # (x * , r * ) samples
1
5 25 100 500 # basis functions
Figure 6: Sensitivity to technical choices: number of level set samples, candidates for robust maximizer, and the basis functions.
inference regret at budget T = 20d depends on these choices for two example targets. We observe that rather small values can be used, speeding up the computation: hundreds of level set samples k and tens of MC samples S are enough in practice, even though we used more in all previous experiments. However, it is important to use sufficiently many RRF basis functions m.
7
CONCLUSION
This work was a part of Finland’s Ministry of Education and Culture’s Doctoral Education Pilot under Decision No. VN/3137/2024-OKM-6 (The Finnish Doctoral Program Network in Artificial Intelligence, AI-DOC). This work was supported by ASM International N.V. The authors acknowledge support from CSC – IT Center for Science, Finland, for computational resources. SK, PM and AK additionally acknowledge the research environment provided by ELLIS Institute Finland and support by the Research Council of Finland Flagship programme: Finnish Center for Artificial Intelligence FCAI. References Matteo Aldeghi, Florian Häse, Riley J Hickman, Isaac Tamblyn, and Alán Aspuru-Guzik. Golem: an algorithm for robust experiment and process optimization. Chemical Science, 12(44):14792–14807, 2021. François Bachoc, Tommaso Cesari, and Sébastien Gerchinovitz. The sample complexity of level set approximation. In International Conference on Artificial Intelligence and Statistics, pages 424–432, 2021.
We introduced a new learning setup for reliable optimization of black box functions when the optimization is carried out in a reliable controlled environment but the eventual solution is deployed in an environment where the inputs may be perturbed. The goal is to find a solution that remains satisfactory, above a given threshold, under maximally large perturbations. We introduced the problem and connected it to the robust Bayesian optimization, explaining how the learning environment is different from previous works. We provided an algorithm that solves the problem using a justified information-gain acquisition, showed how it allocates the evaluations efficiently to explore the level set, and demonstrated the method outperforms standard and robust BO methods in empirical comparisons.
Maximilian Balandat, Brian Karrer, Daniel Jiang, Samuel Daulton, Ben Letham, Andrew G Wilson, and Eytan Bakshy. BoTorch: A framework for efficient Monte-Carlo Bayesian optimization. Advances in neural information processing systems, 33:21524–21538, 2020.
The key technical contribution is the algorithm for sampling robust maximizer candidates for individual GP samples, involving explicit sample-based approximation for the level set. We showed that already rather crude approximations are enough even in dimensionalities higher than typically considered in Robust BO literature, but ultimately any method relying on explicit representation of the level set will fail when the level set is highly complex, i.e. for highly nonsmooth target functions in high dimensions. Despite the strong empirical performance, our analysis provides no convergence guarantees or regret bounds; establishing these is an important direction for future work. However, for such problems there are no other good solutions either. Our method is inefficient when the satisficing threshold is very close to the global function maximizer, but for this regime standard BO is a highly practical solution.
Ilija Bogunovic, Jonathan Scarlett, Stefanie Jegelka, and Volkan Cevher. Adversarially robust optimization with Gaussian processes. Advances in neural information processing systems, 31, 2018.
Felix Berkenkamp, Andreas Krause, and Angela P Schoellig. Bayesian optimization with safety constraints: safe and automatic parameter tuning in robotics. Machine learning, 112(10):3713–3747, 2023. Dimitris Bertsimas, Omid Nohadani, and Kwong Meng Teo. Robust optimization for unconstrained simulation-based problems. Operations research, 58(1):161–178, 2010.
Sathya R Chitturi, Akash Ramdas, Yue Wu, Brian Rohr, Stefano Ermon, Jennifer Dionne, Felipe H da Jornada, Mike Dunne, Christopher Tassone, Willie Neiswanger, et al. Targeted materials discovery using Bayesian algorithm execution. npj Computational Materials, 10(1):156, 2024. Hyunghun Cho, Yongjin Kim, Eunjung Lee, Daeyoung Choi, Yongjae Lee, and Wonjong Rhee. Basic enhancement strategies when using Bayesian optimization for hyperparameter tuning of deep neural networks. IEEE Access, 8:52588–52608, 2020.
Samuel Daulton, Sait Cakmak, Maximilian Balandat, Michael A Osborne, Enlu Zhou, and Eytan Bakshy. Robust multi-objective Bayesian optimization under input noise. In International Conference on Machine Learning, pages 4831–4866. PMLR, 2022.
Juan G Lechuz-Sierra, Ana Elvira H Martin, Ashok M Sundaram, Ruben Martinez-Cantin, and Máximo A Roa. Bayesian optimization for robust robotic grasping using a sensorized compliant hand. IEEE Robotics and Automation Letters, 9(11):10503–10510, 2024.
Carl Doersch and Andrew Zisserman. Sim2real transfer learning for 3D human pose estimation: motion to the rescue. Advances in Neural Information Processing Systems, 32, 2019.
Armin Lederer, Jonas Umlauft, and Sandra Hirche. Uniform error bounds for Gaussian process regression with application to safe control. Advances in Neural Information Processing Systems, 32, 2019.
Lukas Fröhlich, Edgar Klenske, Julia Vinogradska, Christian Daniel, and Melanie Zeilinger. Noisy-input entropy search for efficient robust Bayesian optimization. In International Conference on Artificial Intelligence and Statistics, pages 2262–2272. PMLR, 2020.
Daniel Zhuoyu Long, Melvyn Sim, and Minglong Zhou. Robust satisficing. Operations Research, 71(1):61–82, 2023.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. Advances in neural information processing systems, 31, 2018. Alkis Gotovos, Nathalie Casati, Gregory Hitz, and Andreas Krause. Active learning for level set estimation. In Proceedings of the Twenty-Third international joint conference on Artificial Intelligence, pages 1344–1350, 2013. Daniel Hernández-Lobato, Jose Hernandez-Lobato, Amar Shah, and Ryan Adams. Predictive entropy search for multi-objective Bayesian optimization. In International Conference on Machine Learning, pages 1492–1501. PMLR, 2016. José M Hernández-Lobato, Matthew W Hoffman, and Zoubin Ghahramani. Predictive entropy search for efficient global optimization of black-box functions. Advances in Neural Information Processing Systems, 27, 2014. Carl Hvarfner, Frank Hutter, and Luigi Nardi. Joint entropy search for maximally-informed Bayesian optimization. Advances in Neural Information Processing Systems, 35: 11494–11506, 2022. Donald R Jones, Matthias Schonlau, and William J Welch. Efficient global optimization of expensive black-box functions. Journal of Global optimization, 13(4):455–492, 1998. Leslie Pack Kaelbling and Tomás Lozano-Pérez. Learning composable models of parameterized skills. In 2017 IEEE International Conference on Robotics and Automation (ICRA), pages 886–893. IEEE, 2017. Johannes Kirschner, Ilija Bogunovic, Stefanie Jegelka, and Andreas Krause. Distributionally robust Bayesian optimization. In International Conference on Artificial Intelligence and Statistics, pages 2174–2184. PMLR, 2020.
Rafael Oliveira, Lionel Ott, Vitor Guizilini, and Fabio Ramos. Bayesian optimisation for safe navigation under localisation uncertainty. In Robotics Research: The 18th International Symposium ISRR, pages 489–504. Springer, 2019. Paul Reverdy, Vaibhav Srivastava, and Naomi Ehrich Leonard. Satisficing in multi-armed bandit problems. IEEE Transactions on Automatic Control, 62(8):3788– 3803, 2016. Artun Saday, Yaşar Cahit Yıldırım, and Cem Tekin. Robust satisficing Gaussian process bandits under adversarial attacks. Advances in Neural Information Processing Systems, 38:50381–50408, 2025. Alexander Shapiro, Enlu Zhou, and Yifan Lin. Bayesian distributionally robust optimization. SIAM Journal on Optimization, 33(2):1279–1304, 2023. Benjamin J Shields, Jason Stevens, Jun Li, Marvin Parasram, Farhan Damani, Jesus I Martinez Alvarado, Jacob M Janey, Ryan P Adams, and Abigail G Doyle. Bayesian reaction optimization as a tool for chemical synthesis. Nature, 590(7844):89–96, 2021. Herbert A Simon. A behavioral model of rational choice. The quarterly journal of economics, pages 99–118, 1955. Shion Takeno, Hitoshi Fukuoka, Yuhki Tsukada, Toshiyuki Koyama, Motoki Shiga, Ichiro Takeuchi, and Masayuki Karasuyama. Multi-fidelity Bayesian optimization with max-value entropy search and its parallelization. In International Conference on Machine Learning, pages 9334– 9345. PMLR, 2020. Saul Toscano-Palmerin and Peter I Frazier. Bayesian optimization with expensive integrands. SIAM Journal on Optimization, 32(2):417–444, 2022. Zi Wang and Stefanie Jegelka. Max-value entropy search for efficient Bayesian optimization. In International Conference on Machine Learning, pages 3627–3635. PMLR, 2017.
Dorina Weichert, Alexander Kister, Sebastian Houben, Patrick Link, and Gunar Ernis. Robust entropy search for safe efficient Bayesian optimization. In Uncertainty in Artificial Intelligence, pages 3711–3729. PMLR, 2024. Christopher Williams and Carl Rasmussen. Gaussian processes for regression. Advances in Neural Information Processing Systems, 8, 1995. James Wilson, Frank Hutter, and Marc Deisenroth. Maximizing acquisition functions for Bayesian optimization. Advances in neural information processing systems, 31, 2018. Lin Yang, Junlong Lyu, Wenlong Lyu, and Zhitang Chen. Efficient robust Bayesian optimization for arbitrary uncertain inputs. Advances in Neural Information Processing Systems, 36:43156–43176, 2023.
Maximally Robust Satisficing Bayesian Optimization (Supplementary Material)
Samuli Kinnunen1 1
Petrus Mikkola1
Antti Niskanen2
Arto Klami1
Department of Computer Science, University of Helsinki, Finland 2 ASM International N.V.
A
METHOD DETAILS
A.1
COMPUTATIONAL DETAILS OF THE SAMPLING FROM p(x∗ , r∗ | Dt )
The computational procedure of the sampling algorithm was outlined in 4.2. The method consists of two main stages. First, we construct a differentiable sample path f˜s and approximate the corresponding level set. Second, we initialize robust maximizer candidates and compute the robust maximizer. To approximate the level set (see Algorithm 14), we construct an estimate of the inflated level set Cs = {x ∈ X : |fs (x) − t| < τ } , where τ is a tolerance parameter. The algorithm begins by initializing candidate level set points using Latin Hypercube Sampling (LHS), which ensures space-filling coverage and promotes a uniform exploration of the domain. For each initialized point, we solve a root-finding problem by minimizing the squared residual (fs (xi ) − t)2 . Each point is optimized independently using the Gauss–Newton method. This choice is motivated by computational efficiency, as it only requires first-order derivatives while providing fast local convergence. After constructing the approximate level set samples, we select a maximally diverse subset of these points. From this subset, gradient ascent is applied to generate superlevel set points Ss , which serve as initialization for the robust maximization step. This initialization procedure is described in 19. The robust maximizer is obtained by solving " x∗s = arg max x∈Ss
# X 1 − log exp −β d(x, c) . β
(12)
c∈Cs
In practice, the dominant computational cost arises from evaluating the distances between the n maximizer candidates and the k level set points, which scales as O(knd). When k is large, this step is accelerated by restricting the LogSumExp distance computation to the 100 nearest level set points at each iteration, since distant points contribute negligibly. To preserve accuracy, the nearest-neighbor set is recomputed every 10 iterations. A.2
ADDITIONAL ANALYSIS OF THE LEVEL SET COVERAGE
The coverage of the approximate level set CS determines the quality of the robust maximizer and robustness radius samples {x∗ , r∗ }Ss=1 . We analyze the level set coverage in Fig. 7 for all differentiable synthetic functions used in our experiments up to four dimensions as the analysis requires constructing a dense grid which is computationally infeasible in higher dimensions. The coverage improves with k (the number of level set samples) and decreases with d, as expected.
Algorithm 3: L EVEL S ETA PPROX Input: Sample path f˜s , threshold t, # of initializations n, tolerance ϵ, damping λ, max iters niter Output: Approximate level set samples Cs 1 Initialize Cs = ∅; 2 for k = 1, . . . , n do 3 Sample initial point x(0) ∼ LHS; 4 x ← x(0) ; 5 for j = 0, . . . , niter − 1 do 6 r ← f˜s (x) − t; 7 if |r| ≤ ϵ then 8 Add x to Cs ; 9 break 10 end 11 g ← ∇f˜s (x); r g; 12 x←x− 2 ∥g∥ + λ 13 end 14 end
Algorithm 4: I NITIALIZE M AXIMIZATION Input: Level set samples Cs , sample path f˜s , number of seeds n, learning rate η, ascent steps T , bounds (ℓ, u) Output: Superlevel set points Ss (1) 1 Randomly select c ∈ Cs ; (1) 2 Initialize I ← {c }; 3 for k = 2, . . . , n do 4 For each c ∈ Cs , compute d(c) = minx∈I ∥c − x∥; 5 c⋆ = arg maxc∈Cs d(c); 6 I ← I ∪ {c⋆ }; 7 end 8 Ss ← ∅; (0) 9 for each x ∈ I do (0) 10 x←x ; 11 for t = 1, . . . , T do 12 x ← x + η∇f˜s (x); 13 x ← min max(x, ℓ), u ; 14 end 15 if f˜s (x) ≥ t then 16 Ss ← Ss ∪ {x}; 17 end 18 end 19 return Ss ;
Gaussian Mixture (t=0.8)
100
100
80
80
60 40
d=3 d=4
20 0 102
103
kbase (total level-set seeds = kbase × d) (a) Hartmann
Coverage (%)
Coverage (%)
Hartmann (t=2.0)
60 40
d=2 d=3 d=4
20 0 102
103
kbase (total level-set seeds = kbase × d) (b) Gaussian Mixture
Branin (t=-5.0)
Coverage (%)
100 80 60 40 20 0
d=2 102
103
kbase (total level-set seeds = kbase × d) (c) Branin
Figure 7: Level set coverage. Fraction of the true level set, identified as cells of a dense grid that intersect the level set, covered by the approximation set Cs . A cell is considered covered if at least one sample is in its immediate proximity, up to the grid resolution. Covering level sets in higher dimensions is more difficult, but increasing the sample size improves coverage monotonically.
0.775
0.8
0.550
0.4
0.475
0.2 0.0 0.0
0.400 0.2
0.4
Observations
x1
0.6
0.8
1.0
0.6
0.48
0.4
0.36
x2
GP Mean
x2
0.625
0.8
0.60
0.60
0.700 0.6
0.72
0.72
0.8
Lipschitz (L = )
1.0
0.325
Acquisition Function
0.850
0.6
0.48
x2
1.0
The Effect of Lipschitz Constant (GM 2D, t = 0.8, S = 1) Lipschitz (L = 2 f/ min) 1.0 0.84
Acquisition Function
GP Mean Prediction
0.36
0.4
0.24
0.24 0.2 0.0 0.0
Predicted level set
0.12 0.2
0.4
x1
0.6
0.8
1.0
0.00
0.2 0.0 0.0
0.12 0.2
0.4
x1
0.6
0.8
1.0
0.00
Figure 8: Left: GP mean prediction. The GP posterior mean fitted using 5 observations. Middle: α(x) with small L. Acquisition function with L = 2σf /ℓmin , where σl is prior variance and ℓmin is the smallest length scale over the two dimensions. The queries expand notably outside of the robustness radius. Right: α(x) with conservative L = ∞. Without additional smoothness assumption, the acquisitions are more concentrated and more commonly inside the robustness radius.
A.3
CONDITIONAL ENTROPY WITH LIPSCHITZ CONTINUITY ASSUMPTIONS
Computation of the information gain in Eq. (5) requires approximating the conditional entropy H p(y | x, Dt , x∗ , r∗ ) . We approximate the conditional entropy by leveraging the property that if the distance to the robust maximizer satisfies d(x, x∗ ) ≤ r∗ , then by the definition of robustness, function values are constrained to f (x) > t. On the contrary, when d(x, x∗ ) > r∗ , we do not have such bounds since the function could decrease arbitrarily fast. For sufficiently smooth functions, however, it is possible to construct a lower bound also outside the robustness radius. In the following we discuss one possible approach for doing this, pointing out challenges in obtaining practical gains from such a bound. If we assume the function is Lipschitz continuous with Lipschitz constant L ≥ 0, we can write a common lower bound that holds both within and outside of the robustness radius as H p(y | x, Dt , x∗ , r∗ ) ≈ H(p(y | y ≥ κ(x))) √ κ̃(x) φ κ̃(x) = log 2πe σ(x) Z(x) + , 2 Z(x)
(13) (14)
where φ is the probability density function of the standard normal distribution. Here Z(x) = 1 − Φ κ̃(x) can be computed using the standard normal cumulative distribution function Φ. Now the truncation level κ(x) depends on L: ( κ(x) =
t, ∗
t − L d(x, x∗ ) − r ,
if d(x, x∗ ) ≤ r∗ ,
if d(x, x∗ ) > r∗ ,
used here in a standardized form κ̃(x) = κ(x)−µ(x) . Figure 8 shows how the acquisition function would change when using σ(x) Eq. 14 for computing the information gain; it encourages the algorithm to sample candidates further outside the current robustness radius. Assuming known L, however, is far from realistic and in practice we would need to estimate it. It is possible, but not easy, to estimate the Lipschitz constant L for a Gaussian process proxy [Lederer et al., 2019, Berkenkamp et al., 2023]. Since we only observe f through a finite set of noisy samples, we do not observe L directly but must infer it. Moreover, since f ∼ GP(µt , kt (x, x′ )) is a random function, the constant L is itself a random quantity: different sample paths of the GP admit different rates of change. Consequently, L can only be estimated up to a probabilistic
Log Robustness Regret
GM 3D (t = 0.8, 30 seeds) 10−1
10−2
10−3 0
6
12
18
24
30
36
42
48
Iteration
L=∞
L = 1.6L ̂
L = 3.2L ̂
L = 4.8L ̂
L = 6.4L ̂
L = 8L ̂
Figure 9: Optimization performance with different Lipschitz constant estimates. Typical estimators set L = cL̂ for σ L̂ = ℓmin , where ℓmin is the smallest length scale across input dimensions and c depends on several details of the setup. We intentionally consider also overly optimistic estimators with too small c, to show the method will then fail. The most pessimistic estimator L = ∞ used in all other experiments is highly competitive, but with careful finite estimate (here c = 6.4) we may gain marginal improvement in early stages of the optimization. bound that holds with probability 1 − δ. Lederer et al. [2019] provides a closed-form estimate Lf of a GP sample f s q q q √ 2d ∂ 2 log max k ∂1 (x, x) + 12 6d max max k ∂1 (x, x), r Lk1 x∈X x∈X δ . . Lf = . s q q q √ 2d 2 log max k ∂d (x, x) + 12 6d max max k ∂d (x, x), r L∂kd x∈X x∈X δ
(15)
2
with each element bounding the rate of change along one input dimension. Here k ∂i (x, x′ ) = ∂x∂i ∂x′ k(x, x′ ) is partial i
derivative kernel, L∂ki its Lipschitz constant, r the diameter of the domain X, d the input dimension, and δ ∈ (0, 1) the failure probability. All terms are computable analytically for differentiable kernels, but require kernel-specific derivations and the overall estimate depends also on the optimization domain. Instead of constructing a specific estimator, we demonstrate the effect different estimates would have for the actual optimization quality. For the RBF kernel used in our experiments, we have L ∝ σf /l and hence we plot in Figure 9 the regret as a function of the constant multiplier for this base scale. We see that clearly too optimistic estimators (small L) fail completely and for larger multipliers the convergence curves appear highly similar to L = ∞ that corresponds to arbitrarily non-smooth function and the computation collapses back to theoriginal formulation in Eq. 7. With well-chosen L we see very slight improvement over L = ∞ during the first 20-25 iterations. Given the difficulty of estimating L and the marginal gains that can be achieved, we recommend not attempting to bound the function values outside the robustness radius. A.4
CONDITIONAL ENTROPY WITH NOISY OBSERVATIONS
In the noiseless-observation case in Section 4.1, we used the approximation H p(y | x, Dt , x∗ , r∗ ) ≈ H(p(y | y ≥ κ(x))) ,
so the conditional distribution becomes a truncated normal with an analytical entropy. In the noisy case this implication no longer holds for the observation y; instead it only applies to the latent function value f . Following the derivation of the truncation with noisy observations by Takeno et al. [2020], we derive p(y | x, Dt , f > κ(x)) using Bayes’ theorem: p(f > κ(x) | y, x, Dt ) p(y | x, Dt ) p(y | x, Dt , f > κ(x)) = . p(f > κ(x) | x, Dt ) Here p(y | x, Dt ) is Gaussian under the GP predictive model with observation noise, and p(f > κ(x) | x, Dt ) is the corresponding tail probability for the latent f .
Assume the observation model y = f + ϵ with ϵ ∼ N (0, σϵ2 ), and the GP predictive for the latent f given (x, Dt ) is f ∼ N (µ, σf2 ). Then the joint marginal is 2 σ + σ2 y µ ∼N , f 2 ϵ f µ σf
σf2 σf2
.
Therefore, f | y, x, Dt ∼ N (u, s2 ), with u=µ+
σf2 (y − µ), σf2 + σϵ2
s2 = σf2 −
σf4 . σf2 + σϵ2
Hence, u − κ(x) p(f > t | y, x, Dt ) = Φ , s
µ − κ(x) p(f > t | x, Dt ) = Φ , σf
and the Bayes factor is Z(y) =
p(f > κ(x) | y, x, Dt ) . p(f > κ(x) | x, Dt )
Finally, the conditional entropy is H p(y | x, Dt , f > κ(x)) = −
Z
p(y | x, Dt , f > κ(x)) log p(y | x, Dt , f > κ(x)) dy.
Using p(y | x, Dt , f > κ(x)) = Z(y) p(y | x, Dt ), this becomes Z H = − Z(y) p(y | x, Dt ) log Z(y) p(y | x, Dt ) dy, which can be approximated numerically, e.g. by Monte Carlo with samples yi ∼ p(y | x, Dt , f > κ(x)): N
H≈− A.5
1 X log p(yi | x, Dt , f > κ(x)). N i=1
TECHNICAL HYPERPARAMETERS
Sampling the robust maximizers {(x∗s , rs∗ )}Ss=1 is the most computationally demanding component of our method. Three main technical hyperparameters must be specified. First, the number of Monte Carlo samples S used to approximate the Information Gain in Eq. 5. Second, to obtain differentiable sample paths, the Gaussian process (GP) surrogate is approximated using Random Fourier Features. The number of basis functions m controls the accuracy of this approximation. Third, for each sampled path, we approximate the corresponding level set. The number of the level set samples k affects the accuracy of the estimated robustness radius: denser samples yield more precise estimates. A sensitivity analysis of these hyperparameters is presented in Fig. 10. In this analysis, we vary one hyperparameter at a time while fixing the others to S = 100d, m = 100, k = 100. Overall, MRSBO is relatively robust to the choice of hyperparameters. The main exception is the number of basis functions. In the Gaussian Mixture (3D) problem, using too few basis functions leads to substantially higher regret.
Table 2: Average per-iteration runtime for selecting the next candidate in the Gaussian Mixture 3D benchmark. Method A DVE RS2 MRSBO MES S TABLE O PT
A.6
Mean [s]
Std [s]
2.4 2.4 0.8 7.8
0.4 0.4 0.1 1.5
RUN TIME
In Table 2, we show average per-iteration runtime for selecting the next candidate in the Gaussian Mixture 3D benchmark. Only the computational cost of acquisition optimization (i.e., selecting the next evaluation point) is included. The time required for GP hyperparameter fitting is excluded, as it is identical across all methods. The runtimes are computed using Intel Core Ultra 5 125U × 14 with 16 GiB of RAM memory. The technical hyperparameters are the same as in the experiments specified in Table 4.
B
EXPERIMENT DETAILS
B.1
IMPLEMENTATION
Surrogate Model
We use the same GP surrogate model for all methods: f (x) ∼ GP µ(x), k(x, x′ ) ,
where the prior mean is constant, µ(x) = m, and the radial basis function (RBF) kernel is defined as ! D 1 X ′ 2 ′ 2 (xd − xd ) , k(x, x ) = σf exp − 2 2ℓ d=1
where σf2 > 0 is the output scale, determining the overall magnitude of the covariance, and ℓ > 0 is the length-scale PD parameter, controlling the smoothness of the function with respect to the inputs. The term d=1 (xd − x′d )2 denotes the squared Euclidean distance between x and x′ . The GP hyperparameters include the kernel parameters m, ℓ, σf , and the observation noise standard deviation σ. After each evaluation, these hyperparameters are fitted by maximizing the log marginal likelihood: Z log p(y | X) = log p(y | f (X)) p f (X) | X df. For the length-scale and output-scale parameter, we impose a Gamma prior, ℓ ∼ Gamma(α, λ), with shape parameter α and rate parameter λ as defined in Table 3. The GP model is implemented using SingleTaskGP from BoTorch Balandat et al. [2020] and GPyTorch Gardner et al. [2018]. Optimization of the Acquisition Function For MES and MRSBO, we optimize the acquisition function using multi-start gradient descent Wilson et al. [2018]using the implementation of BoTorch. For AdveRS-2 and StableOpt, due to the discretization-based implementation of AdveRS-2 and the nested max–min optimization in StableOpt, we evaluate the acquisition function on an equidistant grid of size 50d in 2D and 30d in 3D, selecting the input that maximizes the acquisition value. These algorithms were originally developed for adversarial optimization settings, where the maximizer of the acquisition function, xt = arg max α(x), x
10 2 10 3 0
20
40
60
Iteration
80
Log Inference Robustness Regret
Log Inference Robustness Regret
level set samples=150 level set samples=300 level set samples=600 level set samples=1500 level set samples=3000
10 1
level set samples=150 level set samples=300 level set samples=600 level set samples=1500 level set samples=3000
10 1
10 2
0
(a) GM 3D t = 0.8
20
40
60
Iteration
80
(b) Hartmann 3D t = 2
# of (x * , r * ) samples=1 # of (x * , r * ) samples=5 # of (x * , r * ) samples=25 # of (x * , r * ) samples=100
10 1 10 2 10 3 0
20
40
60
Iteration
80
Log Inference Robustness Regret
Log Inference Robustness Regret
(a) Sensitivity to the number of approximate level set samples
# of (x * , r * ) samples=1 # of (x * , r * ) samples=5 # of (x * , r * ) samples=25 # of (x * , r * ) samples=100
10 1
10 2 0
(c) GM 3D t = 0.8
20
40
60
Iteration
80
(d) Hartmann 3D t = 2
10 1 number of basis functions=1 number of basis functions=5 number of basis functions=25 number of basis functions=100 number of basis functions=500
10 2
10 3 0
20
40
60
Iteration
(e) GM 3D t = 0.8
80
Log Inference Robustness Regret
Log Inference Robustness Regret
(b) Sensitivity to the number of robust maximizer samples (r∗ , x∗ )
number of basis functions=1 number of basis functions=5 number of basis functions=25 number of basis functions=100 number of basis functions=500
10 1
10 2 0
20
40
60
Iteration
80
(f) Hartmann 3D t = 2
(c) Sensitivity to the number of basis functions Figure 10: Sensitivity analysis across technical hyperparameters for Gaussian Mixture (GM) 3D (left column) and Hartmann 3D (right column).
Table 3: Gamma priors for the length-scale ℓ and output-scale σf2 for different benchmark functions. The Gamma distribution is parameterized by concentration (shape) α and rate λ.
Function Branin Hartmann 3D Hartmann 4D Hartmann 6D GM 3D GM 4D GM 6D Robot Push 3d
Length-scale ℓ α λ 1.0 4.0 6.0 17.0 6.0 17.0 6.0 17.0 6.0 17.0 6.0 17.0 6.0 17.0 3.0 19.0
Output-scale σf2 α λ 3.0 2.0 11.0 10.0 11.0 10.0 11.0 10.0 11.0 10.0 11.0 10.0 11.0 10.0 3.0 2.0
is adversarially perturbed before evaluation: y = f (xt + δt ) + ϵ. In StableOpt, the perturbation is defined as the worst-case perturbation computed using the GP lower confidence bound LCB = µ(x) − βσ: δt = arg min LCB(xt + δ). δ∈∆ϵ (xt )
AdveRS-2 considers several perturbation schemes. In our experiments, following the experimental setup of Saday et al. [2025], we use random attack perturbations defined as δt ∼ Uniform ∆ϵ (xt ) . The perturbation δϵ ∈ ∆(x) is controlled by a stability parameter ϵ. In our experiments, we set ϵ equal to the admissible perturbation of the maximally robust solution, i.e., ϵ = r∗ . However, in real-world settings, determining an appropriate value for ϵ is challenging, since it requires prior knowledge of the input perturbation levels, which is typically unavailable. Technical hyperparameters. As shown in Section A.5, the proposed method is relatively insensitive to the choice of technical hyperparameters. Therefore, we selected the values primarily according to our computational budget as shown in Tables 4 and 5. Table 4: Technical hyperparameters for problems with dimension d ≤ 3. Hyperparameter Number of Approximate level set base samples Number of basis functions Number of (x∗ , r∗ ) samples
Value 500 100 100
Table 5: Technical hyperparameters for problems with dimension d ∈ {4, 6}. Hyperparameter Number of Approximate level set base samples Number of basis functions Number of (x∗ , r∗ ) samples
B.2
Value 1000 100 100
ROBOT PUSHING BENCHMARK
The robot pushing optimization problem Kaelbling and Lozano-Pérez [2017] takes a pushing action as an input and outputs the object’s distance to the target location. The objective is to minimize the distance to the target location. However, the objective is converted to a maximization problem that takes values f ∈ [0, 5]. The input parameters consist of: rx ∈ [−5, 5]
ry ∈ [−5, 5] rθ ∈ [0, 30], where (rx , ry ) denotes the robot location and rθ denotes pushing angle. We use three different Satisficing thresholds t ∈ (3, 3.5, 4). We adapt the open source implementation by Wang and Jegelka [2017] available at https://github.com/zi-w/ Max-value-Entropy-Search/tree/master/test_functions. B.3
HARTMANN TEST FUNCTION
The d-dimensional Hartmann function (with d ∈ {3, 4, 6}) is defined on the hypercube [0, 1]d . We negate the function to be used in maximization defined as 4 d X X 2 H(x) = αi exp− Aij xj − 10−4 Pij , x = (x1 , . . . , xd )⊤ ∈ [0, 1]d , (16) i=1
j=1
where αi ∈ R, Aij ∈ R, and Pij ∈ R are fixed coefficients. For the six-dimensional case (d = 6), α = (1.0, 1.2, 3.0, 3.2),
10 3 17 3.5 1.7 8 0.05 10 17 0.1 8 14 , A= 3 3.5 1.7 10 17 8 17 8 0.05 10 0.1 14 1312 1696 5569 124 8283 5886 2329 4135 8307 3736 1004 9991 P = 2348 1451 3522 2883 3047 6650 . 4047 8828 8732 5743 1091 381 The function is multi-modal. In the six-dimensional setting, it possesses six local minima and a global maximum at z⋆ = (0.20169, 0.150011, 0.476874, 0.275332, 0.311652, 0.6573), with optimal value H(z⋆ ) = 3.32237. We use the open source implementation of Botorch Balandat et al. [2020]. Details about the adaptation to the lower dimensional problems available at https://botorch.readthedocs.io/en/latest/_modules/botorch/ test_functions/synthetic.html#Hartmann We use three different Satisficing thresholds t ∈ (1, 2, 3) if d ≤ 4 else t ∈ (1, 2, 2.5). B.4
BRANIN TEST FUNCTION.
The Branin function is a two-dimensional test function defined on the domain (x1 , x2 ) ∈ [−5, 10] × [0, 15]. It is given by B(x1 , x2 ) = x2 − bx21 + cx1 − r where the constants are b=
5.1 , 4π 2
c=
5 , π
2
+ 10(1 − t) cos(x1 ) + 10,
r = 6,
t=
1 . 8π
(17)
(18)
The function is multimodal and possesses three global minimizers, z1 = (−π, 12.275),
z2 = (π, 2.275),
z3 = (9.42478, 2.475),
(19)
at which the global minimum value is B(zi ) = 0.397887,
i = 1, 2, 3.
(20)
We use three different Satisficing thresholds t ∈ (−20, −10, −5). We use the open source implementation of Botorch Balandat et al. [2020] available at https://botorch. readthedocs.io/en/latest/_modules/botorch/test_functions/synthetic.html#Branin . B.5
GAUSSIAN MIXTURE TEST FUNCTION.
The Gaussian mixture test function is defined on the unit hypercube x ∈ [0, 1]d , with d ∈ N. It consists of a mixture of Gaussian components, a single mode for each corner and additional mode in the center. Modes are defined as: M X 1 ⊤ −1 (21) f (x) = ai exp − (x − µi ) Σi (x − µi ) , 2 i=1 where M is the number of modes, ai > 0 are amplitudes, µi ∈ Rd are the component means, and Σi ∈ Rd×d are diagonal covariance matrices. The mean vectors of the Gaussians are given by the 2d shifted corner points and one shifted center point. Let v(k) ∈ {0.2, 0.8}d for k = 1, . . . , 2d enumerate all corner vectors, and let 1 ∈ Rd denote a vector of ones. Define the corner means µk = v(k) + 0.15 1,
k = 1, . . . , 2d ,
(22)
and define the center mean µ2d +1 = 0.5 1 + 0.15 1.
(23)
Thus, the total number of modes is M = 2d + 1. Base variances. All modes are initialized with σi2 = 0.02 and then modified for selected modes: the first corner mode is 2 set to σ12 = 1.2 · 0.02, the corner mode with the highest amplitude uses σm = 0.5 · 0.02 = 0.01, and the center mode uses 2 σc = 0.01. Dimension-dependent scaling. The covariance matrices are isotropic and dimension-scaled as Σi = σi2 dα Id , where Id is the d × d identity matrix and
1.0, 0.8, α= 0.85, 0.6,
(24)
d = 2, d = 3, d = 4, d = 6.
Amplitudes. All amplitudes are ai = 1 except for the corner mode with the highest amplitude, where am = 1.1, and the center mode, where ac = 0.4. We use three different Satisficing thresholds t ∈ (0.7, 0.8, 0.9) if d ≤ 4 else t ∈ (1.0, 1.1, 1.2).
C
ADDITIONAL RESULTS
The complete experimental results are shown in Fig. 11 and Fig. 12. For each benchmark and method, we report Robust Inference Regret curves over the full optimization horizon. These results complement the findings presented in Section 6 by providing a more comprehensive view of the empirical performance across all settings. In addition, we include a development view of the optimization process in Figure 13, illustrating how the proposed acquisition criterion encourages exploration during the early stages and gradually shifts toward exploitation by reducing uncertainty in the relevant regions of the level set.
10−1
0
15
30
45
60
75
Log Inference Robustness Regret
Branin 2D, t=-10 Log Inference Robustness Regret
Log Inference Robustness Regret
Branin 2D, t=-20 100
100
10−1
10−2
90
0
15
30
Iteration
10−3 45
60
75
10−3
90
0
15
30
75
0
15
30
75
15
30
45
60
75
45
60
75
15
30
0
15
30
45
60
MES
AdveRS2
60
75
90
75
90
100
10−1
0
75
90
Iteration
MRSBO
45
15
30
45
60
Iteration
10−2
90
90
10−3
90
10−1
Iteration
75
10−2
0
Log Inference Robustness Regret
Hartmann 3D, t=2 Log Inference Robustness Regret
Hartmann 3D, t=1
10−3
60
RP 3D, t=4
10−1
0
45
Iteration
Iteration
30
30
10−1
90
100
90
10−2
15
75
Iteration
10−1
0
15
Gaussian Mixture 3D, t=0.9
Log Inference Robustness Regret
Log Inference Robustness Regret
Log Inference Robustness Regret Log Inference Robustness Regret
60
0
RP 3D, t=3.5
10−1 45
60
90
10−3
Iteration
100
30
45
75
Iteration
10−3
90
60
10−2
90
10−2
RP 3D, t=3
15
75
10−1
Iteration
0
60
Log Inference Robustness Regret
Log Inference Robustness Regret
Log Inference Robustness Regret
10−3 60
45
45
10−1
Gaussian Mixture 3D, t=0.8
10−2
45
30
Iteration
10−1
30
15
Gaussian Mixture 2D, t=0.9
10−2
Gaussian Mixture 3D, t=0.7
15
0
Iteration
10−1
Iteration
0
90
Log Inference Robustness Regret
Log Inference Robustness Regret
Log Inference Robustness Regret
10−2
30
75
10−2
Gaussian Mixture 2D, t=0.8
10−1
15
60
10−1
Iteration
Gaussian Mixture 2D, t=0.7
0
45
Branin 2D, t=-5
100
Hartmann 3D, t=3
10−1
10−2
0
15
30
45
60
Iteration
StableOpt
Figure 11: Comparison of AFs d ≤ 3
75
90
240
320
400
480
0
240
320
400
240
320
400
480
160
3 × 10−2 2 × 10−2
10−2
560
0
80
160
240
320
160
240
320
400
400
480
10−2
560
0
80
160
240
320
400
GM 4D, t=0.7
GM 4D, t=0.8
GM 4D, t=0.9
400
480
Log Inference Robustness Regret
Iteration
320
10−2
10−3
560
0
80
160
240
320
560
10−1
Iteration
240
480
Hartmann 6D, t=2.5
10−2
0
80
Iteration
10−1
560
4 × 10−2
Iteration
Log Inference Robustness Regret
160
480
6 × 10−2
Log Inference Robustness Regret
Hartmann 6D, t=2
400
480
480
560
480
560
10−2
560
0
80
160
240
320
400
Iteration
GaussianMixtureDense6D, t=1.0
GaussianMixtureDense6D, t=1.1
GaussianMixtureDense6D, t=1.2
10−2
80
160
240
320
400
480
560
Log Inference Robustness Regret
Iteration
Log Inference Robustness Regret
Iteration
10−1
0
160
Hartmann 6D, t=1
10−3
80
80
Iteration
10−2
0
10−2
Iteration
10−2
80
10−1
560
Log Inference Robustness Regret
Log Inference Robustness Regret
160
10−1
0
Log Inference Robustness Regret
80
Hartmann 4D, t=3 Log Inference Robustness Regret
10−2
0
Log Inference Robustness Regret
Hartmann 4D, t=2 Log Inference Robustness Regret
Log Inference Robustness Regret
Hartmann 4D, t=1
10−1
10−2
0
80
160
240
Iteration
320
400
480
560
Iteration MRSBO
10−1
10−2 0
80
160
240
320
Iteration MES
Figure 12: Comparison of AFs 4D and 6D.
400
480
560
0.850 0.775
0.8
0.700
GP Mean
0.6
0.99 0.90
0.8
0.81 0.6
0.72
0.625 0.550
0.4
0.475 0.2 0.0 0.0
Mean Prediction
1.0
GP Mean
Mean Prediction
1.0
0.63 0.4
0.54 0.45
0.2
0.400 0.2
0.4
Observations Predicted level set True level set
0.6
0.8
1.0
0.24 0.21
0.8
0.0 0.0
0.2
0.4
Observations Predicted level set True level set
Predicted robust maximum Inference robustness radius Candidate
Acquisition Function
1.0
0.325
0.36 0.6
0.8
1.0
Predicted robust maximum Inference robustness radius Candidate
Acquisition Function
1.0
0.21 0.18
0.8
0.15
0.15
Acquisition Function
Acquisition Function
0.18 0.6
0.6
0.12
0.12 0.4
0.09 0.06
0.2
0.09
0.4
0.06 0.2
0.03
0.03 0.0 0.0
0.2
0.4
0.6
(a) Iteration 1
0.8
1.0
0.00
0.27
0.0 0.0
0.2
0.4
0.6
0.8
1.0
0.00
(b) Iteration 9
Figure 13: Progress of MRSBO algorithm in 2-dimensional Gaussian Mixture with t = 0.7 and random initialization of 5 evaluations. The top panel shows the GP mean with predicted level set (red dashed), true level (blue dashed), inference robust maximizer (green marker), the selected candidate (black triangle), and observations (red points). The bottom panel shows the acquisition function with the selected candidate (black triangle). The model increasingly refines the level set approximation and concentrates exploration near the borders of the estimated maximally robust satisficing region.