Bayesian Optimization with Informative Covariance
Authors: Afonso Eduardo
Michael U. Gutmann
Venue: Accepted by TMLR ()

Abstract: Bayesian optimization is a methodology for global optimization of unknown and expensive objectives. It combines a surrogate Bayesian regression model with an acquisition function to decide where to evaluate the objective. Typical regression models are given by Gaussian processes with stationary covariance functions. However, these functions are unable to express prior input-dependent information, including possible locations of the optimum. The ubiquity of stationary models has led to the common practice of exploiting prior information via informative mean functions. In this paper, we highlight that these models can perform poorly, especially in high dimensions. We propose novel informative covariance functions for optimization, leveraging nonstationarity to encode preferences for certain regions of the search space and adaptively promote local exploration during optimization. We demonstrate that the proposed functions can increase the sample efficiency of Bayesian optimization in high dimensions, even under weak prior information.

Bayesian Optimization with Informative Covariance
Afonso Eduardo
afonso. www..eewwwwwddddd....o
School of Informatics University of Edinbungh
Michael U. Gutmann
michael.gutmannnBeee
School Informatics
University of Edinburgh
Abstract
Bayesian optimization is methodology for global optimization of unknown and expensive objectives. It combines surrogate Bayesian regression model with acquisition function to decide where to evaluate the objective. Typical regression models are given by Gaussian processes with stationary covariance functions. However, these functions are unable to ex- press prior input-dependent information, including possible locations of the optimum. The ubiquity of stationary models has led to the common practice of exploiting prior information informative mcan functions. In this paper, we highlight that these models can perform poorly, especially in high dimensions. We propose novel informative covariance functions for optimization, leveraging nonstationarity to encode preferences for certain regions of the search space and adaptively promote local exploration during optimization. We demonstrate that the proposed functions can increase the sample efficiency of Bayesian optimization in high dimensions, even under weak prior information.
1 Introduction
Bayesian optimization (BO) is methodology for global function optimization that has become popular since the work of Jones et al. (1998). Other influential works date back to Kushner (1962) and Močkus (1975). BO has been applied to broad range of problems (sce c.g. Shahriari al., 2016), including chemical design (Gómez- (Gómaz-Bombarelli et al., 2018), neural architecture search (White et al., 2021) and simulation-based inference (Gutmann et al., 2016; Cranmer et al., 2020).
As a Bayesian approach, the idca is to place a prior over the objective function, typically a Gaussian process (GP) prior characterized by a mean and covariance function. The prior, combined with an observation model and past. function evaluations (observations), allows the computation of a posterior predictive distribution. The next acquisition is determined by the posterior predictive and a utility function, treating the optimization decision problem.
Standard BO is generally effective in relatively low-dimensional problems. However, as dimensionality in- creases, statistical and computational problems become more noticeable. In particular, the curse of dimen- sionality manifests itself while learning the probabilistic surrogate of the objective function and during the acquisition process. Indeed, as training data become sparser, the estimation of an accurate surrogate becomes more difficult and the reliance on extrapolation incrcases. In this context, previous work has relied on certain structural assumptions (Binois & Wyeoff, 2022), namely low elfective dimensionality (Garnetl al., 2014; Li et al., 2016; Wang al., 2016; Letham et al, 2020; Moriconi et al., 2020; Raponi et al., 2020; Eriksson & Jankowiak, 2021; Grosnit et al., 2021) and additivity (Kandasamy et al., 2015; Gardner et al., 2017; Mutný & Krause, 2018; Rolland et al., 2018; Wang et al., 2018). Local methods based on space-partitioning schemes and trust. regions have also been proposed (Assael et al., 2014; McLeod et al., 2018; Eriksson et al., 2019; Wang et al., 2020b; Diouane et al., 2021). Notably, these methods employ stationary covariance functions.
By construction, stationary covariance functions are translation invariant and thus unable express spatially-varying information. As shown in Figure 1, stationary models can lead to poor optimization of

---

Stationary
Stationary Quadratic Mean
Monstationary IteX0
Ground Truth
Creditle
Initial Chaurvations
Acouiskions
Figure Overconfident models are unable to capture the ground truth, causing BO to yield suboptimal results. The stationary model with informative quadratic mean fits the trend too well, resulting in small residuals. These residuals lead to smaller prior variance and longer lengthscale that, in turn, decrease epistemic uncertainty. Stationarity posits that these hyperparameters are constant, but this assumption does not hold for functions with irregular bumps. The proposed informative priors, being nonstationary, allow for the specification of spatially-varying properties.
relatively simple onc-dimensional objectives, suggesting issues in higher- sigher-dimensional domains. A common practice is to include input-dependent information via informative mean functions (see e.g. Snoek et al., 2015; Acerbi, 2019; Järvenpää et al., 2021). However, as demonstrated in Figure 1, this practice does not necessarily alleviate the limitations of stationary covariance functions. In contrast, nonstationarity allows for the incorporation of input-dependent information directly in covariance functions. The crux lies in the design of flexible and scalable informative GP priors
In the remainder of this paper, after we introduce the background in Section 2, we lay the foundation for our methodology by discussing the advantages of second-order nonstationary for optimization in Section 3. focusing on spatially-varying prior variance, lengthscales and high-dimensional domains. This includes a discussion regret bounds for covariance functions with spatially-varying parameters, which we have not previously found in the literature. We then make the following contributions:
In Section 4, we design informative covariance functions that leverage nonstationarity to incorpo- rate information about promising points. Our GP priors encode preferences for certain regions and allow for multiscale exploration during optimization. To our knowledge, we are the first to consider informative GP priors for optimization with spatially-varying prior variances and lengthscales, both induced by shaping function that captures information about promising points.
In Section 5, we empirically demonstrate that the proposed methodology can increase the sample efficiency of BO in high-dimensional domains, even under weak prior information. Additionally, we show that it can complement existing methodologies, including GP models with informative mean functions, trust-region optimization and belief- angiee-e acquisition functions.
Section 6 discusses related work and Section 7 concludes the paper.
2
Background
2.1
Bayesian Optimization
In global optimization, the goal is to find the global minimizer of an objective function f: X + Y R,
x arg mig f(x),
#EX
where the domain X C R,,, for objectives of D variables, and f is function whose expression and gradient are lypically unavailable. Evaluation is possible al any point x; E X, yielding the observation yı = -(:;;, but its cost is high enough to justify the use of sample-efficien methodologies that aim to minimize the number of evaluations by careful selection. In this context, Bayesian optimization (BO) emerges as competitive approach. The idea is train a Bayesian regression model A that explains the evidence collected up step n, D,, {( yr) igno and provides a distribution over possible functions. Given this probabilistic
2

---

Algorithm Bayesian Optimization (BO)
Input: objective lunction S acquisition function a, statistical model M, initial evidence Dno
repeat
#n+1 arg vegmaxale aax M)
Find best candidate
Un+1 f2n+11
Evaluate candidate
Dn+1 D, U 1, Un+1)
Update evidence set
until stopping condition is met
surrogate, an acquisition function a assigns a utility each candidate point, guiding the selection of the next acquisition, i.e, the next point to evaluate. The procedure is outlined in Algorithm 1, where the stopping condition may be based on maximum evaluation budget and tolerance. For more details on BO, see e.g. Shahriari et al. (2016), Frazier (2018) and Greenhill et al. (2020).
2.2 Gaussian Process Regression
BO is traditionally paired with Gaussian process (GP) models. For given GP prior f GP( Ce) with mean function mg: X Y, covariance function Co: X x X + Rt, and an additive Gaussian observation model - f()) +,, , ~ ( the univariate posterior predictive distribution of the latent f at test input æ is Gaussian with mean m, and variance of
mn (a)) mo. (a)) Ca (x) + +
m.).
(2)
(xr) = C,(
C₂ (2,, x. = Co. (x;,, x.)) - Cn (x,) + +
(4)
C, (x) = (Co. (x, x. ..... Co. (T))
(5)
where the entries in the Gram matrix Cn are [Cr]us = (x, Rz), and the mean vector mp has en- E D,,. Hyperparameters 0" can, c.g., be fixod a priori or learned via empirieal Bayes, by introducing a hyperprior p(0) and maximizing the unnormalized posterior p(0)D..) p(y Xn, ())( where the marginal likelihood is p(y Xx, 0) my Cn off. Xn {x.)isn For comprchensive introduction to GPs, see Rasmussen & Williams (2005).
tries [m,,]: me- (x;). Y(z,,y)
2.3 Stationary Covariance Functions
Most popular covariance functions are translation invariant, depending only the relative position of the input points. These functions compute (fxff based on translation-invariant distance between and æz. Valid covariance functions must be symmetric and positive definite. In this family, we find the squared exponential or Gaussian covariance,
Ca(o,, = of exp( (z),)) (x,, V(: IT
xz)),
(6)
where of is the prior variance and dx is the Mahalanobis distance function with matrix A. For simplicity and scalability, it is customary to constrain the matrix to be diagonal with entries [A]ad = X e R d E (1, D), known as (squared) lengthscales. In this case, the Mahalanobis distance turns into the weighted Euclidean distance, which we denote by dwg. These lengthscales control the exponential decay of correlation. For instance, a short. lengthscale leads to sample paths, realizations, that vary more rapidly.
The Gaussian covariance, being infinitely differentiable, may generate unrealistically smooth sample paths. A more general class of functions is generated, or reproduced, by the Matérn covariance function which controls the differentiability of the process via smoothness parameter v, recovering CG when V + oc. For optimization purposes, good empirical results have been found with v = 5/2 (Snoek et al., 2012), for which
Cma(;; 1.) ob(1 1+ V5draa2s my) ( (2,, mzy)
(-V5ddn
(7)
3

---

An attractive property of the Gaussian and mixtures thereof, such Matérn, is their proven universality (Micchelli et al., 2006). Asymptotically, these covariance functions, or kernels, are capable of approximat- ing any contimuuus function on any compact subset of the domain. In practice, however, other covariance functions may be preferable in the small sample regime, espocially if the clicited prior includes relevant information, or structure, that accelerates the learning process.
2.4 Acquisition Functions
While many rules to determine promising candidate solutions have been proposed (see e.g. Wilson et al., 2018; Neiswanger et al., 2022), those with an analytical expression remain the most popular. Such acquisition functions can be evaluated without Monte Carlo or quadrature approximations. Given a Bayesian model, myopic acquisition functions suggest one acquisition at time, depending only on the univariate posterior predictive distribution with mean mn and variance Un: i.e., a(x Dn, = a(x may m,) In this context, popular choice is the Lower Confidence Bound (LCB) critcrion (Srinivas et al., 2010),
LCB(x) = m,, - (x),
(8)
where the factor 3n is related to a confidence level, balancing exploitation and exploration. Another popular alternative is the Expected Improvement (EI), which has been found to be more effective than LCB when properties of the objective function, e.g. norm bounds, are unknown (Snoek et al, 2012; Chowdhury & Gopalan, 2017; Merrill et al., 2021). This acquisition function is defined as
EI(x)) = ()((())) ([)) [xx]] = (,,())((2())) <<<(x)) = (xx)) + +((=()):0, 1), x(x) (f(bbttt) - mn (x)))/o(xx).
(9)
(10)
where I(x:) = max(0, f f( is the improvement over the incumbent (* and FNr is the standard Normal cumulative distribution function (CDF). For noisy functions, (Xbeet) is replaced by the minimum posterior predictive mean, mn = min (x))
2.5 Regret
In the context of optimization, there are several useful metrics to measure performance. The instantaneous regret is the loss incurred at step n, Tn = f(æ,) - ff("), where H, is the acquisition and xt is the global minimizer of f. The simple regret is the minimum instantaneous regret incurred up to step or, equivalently, the regret incurred by the incumbent solution. The cumulative regret is defined as the sum of instantancous regret over N steps, Ry EnSN Tn- If the cumulative regret grows at sublinear rate, then the average regret goes zero, Rw/N implies that simple regret goes to zero, ensuring that the global optimum is asymptotically found. For more details regret bounds, see e.g. Vakili et al. (2021) and Gupta et al. (2022).
= minegn Tt,
O, which is also known as no regret. This property
3 Benefits of Second-Order Nonstationarity for Optimization
We now turn our attention to GP priors with spatially-varying properties, namely prior variance and length- scale, which provide the foundations for the informative covariance functions in Section 4. While an analysis in terms of cumulative regret. is customary, we focus on the instantaneous regret because we are interested in the relatively small sample regime, for which the asymptotics are not relevant. Furthermore, by demon- strating that the instantancous regrets In are smaller (or larger), follows that the cumulative regret Ry is also smaller (or larger).
Spatially-varying prior variance The prior variance of quantifies the (epistemic) uncertainty of the surrogate model of f before observing data. Relatedly, posterior predictive variances (x),, defined in Equation 3, determine the mutual information between observations yn and the objective function f. MI (yx:f) = 2 ΣN log (1 + 2 (0_=20=_((,,)) (x,)), see e.g. (Srinivas et al., 2010, Lemma 5.3). Then, for sta- tionary covariance functions, any point 2n, where (x,) ~ 0 provides equal information about f because
The L.B needs be minimized and thus corresponds a negative acquisition function a.
4

---

Uninformative Mean
Informative Mean
Ground Truth
(Statiinary) (Norstationay)
Figure 2: Narrower credible intervals (CT) can be obtained with a spatially-varying prior variance. Typically, GP priors BO uninformative, characterized by constant prior mean functions and stationary covariance functions that posit constant prior variance. Then, for popular acquisition functions, all candidates are considered equally good priori. Conversely, spatially-varying prior variances can encode preferences. In this example, points near the center are more informative priori.
Constant Lengthscale
vpatial Lengthscale
Figure 3: Given the evaluation budget (N 20), model with spatially- lengthscales better approximate the function than one with constant lengthscale. While the latter must set shorter global lengthscale, a multiscale approach can locally assign shorter lengthscales.
the predictive variance is approximately constant, ( oî. Conversely, surrogate models of with spatially-varying prior variance spatially informative, in the absence of neighboring observations, in which case (xx) x (æx).
In terms of optimization, spatially-varying prior variance can express a priori preferences for certain regions, as illustrated in Figure 2. Notably, compared to the stationary case, better optimization performance can be achieved: For popular acquisition functions, the worst-case instantaneous regret is proportional to the predictive standard deviation (Lemma A.1), provided that the interval [f(x) m,,(=) s (,),( holds for some factor Bn. Then, recall that the predictive variance is upper bounded by the prior variance, of (z). While the latter is constant for stationary covariance functions, covariance functions with ipatially-arrying prior variances allow for narrower intervals and, in lurn, lighter regret bounds. As caveat, standard regret analysis assumes that 1) intervals include the objective function and 2) global optima of acquisition functions can be found, but in practice these assumptions do not necessarily hold.
Spatially-varying lengthscale lerms sample efficiency, stationary covariance functions are ill-suited for approximating objectives with spatially-varying properties. Intuitively, a function that varies rapidly in one region may be characterized by short lengthscales, but is, otherwise, more efficiently described by longer lengthscales. In order to capture the finest details of such function, stationary covariance would require the shortest possible lengthscale, resulting in additional observations to approximate the function and ultimately find the optimum. In contrast, multiscale approach only requires more neighboring observations where the local lengthscales are shorter (Figure 3).
In more detail, let the objective be a realization from a locally stationary process, (<=)= Im with disjoint X C X, m (1,..., M), M E N+. Each component is member of a class of functions forming the Hilbert space Hc.. (Xm) with reproducing stationary kernel Om and bounded norm || S Bm If f is to be estimated with stationary covariance function C, then it must reproduce the space (x( 2 U Ho (Xm) with E UmSM Xm and norm bound 2 vΣm B2, by (Aronszajn, 1950, Theorem $6). This corresponds to search in the space of more complex functions, for which more
5

---

observations are required because the global lengthscale must be set as the shortest from {Cm JMS Indeed, for shorter lengthscales, X S X, we have larger posterior predictive variances, xax) 2 o (æ), which lead to greater information gains MIn. 2 MI and larger norms, |l 2 |Ilfl by (Bull, 2011, Lemma
4), contributing to a larger BR,x' 2 Bn,A because 3n = Sup llfllltooo ,o,, ((MI, + - log 6), d E (0, 1], by (Chowdhury & Gopalan, 2017, Theorem 2). Then, by Lemma A.1, wider intervals [f(x)) - m. (x)] S f, (2:) lead to larger regret bounds, and hence worse performance.
The argument above provides additional justification for spatially-varyingg prior variances, in which case each Cm is characterized by different but constant &.mm Then, the stationary covariance function must set of maxm o.m: leading again to wider intervals, and hence worse performance by Lemma A.1. Finally, as M grows, it becomes more challenging to estimate and maintain M local models. Nonstationary approaches that do not rely on space partitions can avoid these computational and statistical problems. Notably, in the limit, each subset is a singleton, Xm = {zm]. Vam E X, and the restriction f| is reproduced by (m æ.y) = Jm with norm |lf [f(m)) For an uninformative prior mean, mo(#) = O, the narrowest interval, before observing data, (f(mm) - (2m))
00 is obtained with a spatially- or, equivalently, with Pn 1 and spatially-varying prior standard deviation (xx [f(x)]. Thus, the optimal prior variance depends on the shape of the objective function, as previously suggested by Figure 2.
varying fo(x) [f(x) and constant on
High-dimensional domains In high-dimensional Euclidean spaces, most of the volume is on the boundary, and for stationary GPs, epistemic uncertainty, measured by the posterior predictive variance, increases with the distance to the training data, tilting the acquisition step toward blind exploration of the boundary, typically faces of a hypercube (Binois & WycofT, 2022). This problem is known as the boundary issue (Swersky, 2017; Oh et al., 2018).
Morcover, uninformative priors provide no practical finite-time convergence guarantecs. Begin by recalling that the motivation behind BO is the optimization of expensive functions, which puts a practical constraint on the size of the evaluation budget. However, unless the total budget N increases exponentially, the space cannot be covered in such a way that a global contraction of the posterior predictive variance is ensured, see e.g. (Kanagawa et al., 2018, Section 5.2) and (Wüthrich et al., 2021, Appendix C). Indeed, in the relatively small sample regime, there is at least one point F. whose predictive variance is approximately equal to the (constant) prior variance, (zx) og. This is explained by the absence of neighboring observations, in which case the training data provides a negligible reduction of uncertainty because (x-) R 0. As a result, the worst-case instantaneous regret, being proportional to the predictive standard deviation (Lemmaa A.1), remains approximately constant for fixed Sn, or alternatively is nondecrcasing for nondecrcasing Bn. (Chowdhury & Gopalan, 2017, Theorem 2), and thus there is no guarantee of improvement.
These problems can be alleviated by more informative GP priors. In particular, spatially-varying prior variances allow the surrogate to remain spatially informative, even in the absence of neighboring observations, and ppatiall-vvaryine lengthscales can encourage greedier strategies, in an attempt to improve upon the incumbent solution when given small budgets.
4 Informative Covariance Functions
For the construction of our informative covariance funetions, we start with set of anchors (m) ), i.e., set of promising candidate solutions. This prior belief can be interpreted as point-mass distribution over promising points. A more conservative prior p(x") is obtained by adding an uninformative slab that ensures the optimal solution is included in its support, and forming a kernel density estimate by associating (!)
non-negative weight wi, distance function di and kernel ki with each anchor x
1
1 Σ [ - 1) kn (dile", (0(==)=x()) :0))) L LSL
p(x*) x x(x*),
6(z))=1+ 1+
(11)
Distance functions and kernels characterize the neighborhood of anchors. Larger weights correspond to the belief that better candidates can be found in a given neighborhood. Next, we introduce spatially-varying
6

---

A,= 0.3
1= 0.1
Aa
0.1
0.1
Prior Variance (d)
0.3
0.5
0.5
6
0.8
0.8
1.0
4-
2-
1.00
-3.75
-3.50
0.50
5.75
2.00
0.75
1.00
-1.00
-0.75
-0.50
-233
0.00
0.25
0.50
0.75
1.20
Figure 4: For a Gaussian kar and weighted Euclidean distance d, = dwE, decreasing the ratio Tr increases the prior variance around the anchor 0. Stationarity is recovered when To = 1. The lengthscale Ag in dwE controls the rate at which the prior variance decreases to the stationary baseline.
prior variances by rescaling the stationary covariance function Cs(2,, (g))
Ci(,,, = oz a,) (z; Ry),
x.;) = ob ((1))) o M(x,),
(12)
where op is a scaling constant and Cz a) = 1. The modulating oz(a,, is symmetric and separable in æs and Ej and hence a valid covariance function (Genton, 2001). The informative covariance Ci xz); being the product two covariance functions, is also a valid covariance function.
The prior covariance function of (.), is generalization of the constant prior (co) oy. By recalling that covariance functions compute Cfv(f(x;)) fa;,)) the intuition is follows: For two points æ; and æj that are close to anchors, i.e., large ((x;) and Φ(æ;)) ), this prior posits that their values, f(x;) and f(æ;), are highly correlated because both should be close to the minimal value and hence similarly small. Then, if the probability for a point, e.g. Tj, to be a good candidate decreases, the covariance becomes smaller because fxz)); being less constrained, can take on greater range values.
In Section 5, we focus on the special case of single anchor, for which we obtain
o n.) = ob/ 1+ (1_
To -1) li (do (2:, 200) _- ro
-111 1) ho (do (2,,,
1+
(13)
where ki k1, do di and To A 1/w, E (0, I]. We find this parametrization more convenient because inverse weights are bounded. The ratio ro controls the degree of stationarity, with ro resulting in a stationary covariance and r. + 0 in increased nonstationary effects, as demonstrated in Figure 4.
Next, we introduce spatially-varying lengthscales via input warping. In input warping (Snoek et al., 2014; Risser, 2016). the stationary covariance function Cs is computed in a transformed space, i.e., after applying a (possibly nonlinear) warping function hy to æ. Endowed with such a transformation, the informative covariance function Ci(,,, 2,), in Equation 12, is given by
a,) = af m, hxx)))
(14)
Recall that inputs in (2,, 2,) are already being compared to anchors in the shaping function & (Equa- tion 11). Now, notice that anisotropic stationary covariance function, equipped with the Mahalanobis distance dm, is equivalent to isotropic covariance function, equipped with the standard Euclidean distance and linear transformation hx(x)) = A- Σ... Then, since the warping function can be any real vrctor-value function, we propose nonlinear transformation that shrinks the lengthscales da around anchors, thereby expanding their neighborhoods. In particular, for a single anchor, we have
hx()) = u (.)
"사 (x) 1 + (rx - Id mo));
(15)
with ratio r> 1).2 Note that the matrix A is diagonal with squared lengthscales X and the neighborhood of 2n is characterized by the kernel ly and distance d2. The shortest lengthscale Vrxud-
SISL (1/wn - 1)As (di (zz a The general case with multiple anchors can be given by the kernel mixture ux = 1+1//L EIs an expression similar that of the shaping function & Equation 11. The inverse weights 1/w E (0, 1] shorter Iengthscales.
7

---

5 Experiments
5.1 Setting
Informative covariance functions The shaping function & should capture the existing knowledge about the objective, possibly by being elicited from experts or learned on problems that are known to be similar. This function, based on kernel mixtures, can express arbitrary beliefs, but its optimality crucially depends on its alignment with the shape of the objective function, as discussed in more detail in Scction For black-box optimization, in which little is known in advance, a default kernel and distance function may be used, and their hyperparameters learned by empirical Bayes. To avoid overfitting and reduce computational complexity, regularizing hyperpriors and hyperparameter tying are effective techniques that may also be used. In particular, the parameters and functions that control the spatially-varyingg prior (co)variacce and lengthscales can be tied by setting the ratio rp _ To = rA, kernel ko _ hig = 1x and distance do _ d, = d,- In our experiments, we specify a Gaussian ko, equipped with a weighted Euclidean distance dg = dwE and a shared lengthscale vector, which also parametrizes A. Dirac delta hyperpriors for anchors have also been adopted, but note that none of these choices are strictly required. More implementation details can be found in Appendix B, and we refer to Appendix F for the settings of ro, where we perform sensitivity analysis.
An overview of the main methods using the proposed covariance functions (I) is included in Appendix D. These methods are I+XO and I+XA, where +XO denotes a fixed anchor 2o (placed at the center of the search space) and +XA denotes an adaptive anchor (placed at the incumbent solution). More methods can be developed by incorporating complementary methodologies. In this paper, we test I+XA+TR (trust-region optimization, Appendix D.1), I+XA+QM (informative mean functions, Appendix D.2), I+XA+SAAS (sparsity- inducing hyperpriors, Appendix H) and I+XA+GKEI (belief-augmented acquisition functions, Appendix T).
Baselines Previous work has assumed that optimal solutions are generally closer to the center of the search space rather than its boundary. As dimensionality increases, BO with a stationary covariance (8) may allocate a large portion of the budget toward blind exploration of the boundary, which by assumption does not contain the optimum, resulting in worse performance. To mitigate boundary overexploration, two approaches have been proposed: quadratic mean functions (+QM) (Snoek et al., 2015) and cylindrical covariance functions (C) (Oh et al., 2018). In particular, the latter are nonstationary covariance functions that aim to expand the innermost region, and have been shown to outperform stationary additive models scveral high-dimensional additive objectives, without having to infer their structure. Further details can be found in Appendix C. Additionally, the use of trust regions (+TR), despite not specifically intended to address the boundary issue, has become a popular acquisition strategy for high-dimensional BO. This strategy can promote exploitation by constraining acquisitions to an adaptive region centered at the incumbent solution.
In our experiments, we refrain from direct comparisons with existing BO packages due to the diverse training and acquisition routines. Instead, we reimplemented the BO baselines to study the effect of different GP priors on optimization performance under controlled conditions. We also investigate the interplay between GP priors and trust regions. Our implementation of +TR follows TuRBO-1 (Eriksson et al., 2019). The GP models are global, trained on all observations, and the trust region is given by a box whose side lengths are /hubled/halled after consecutive successes successes/failures failures. The hyperparameter values of +TR are identical to those of TuRBO-1, except for the failure tolerance, which is originally given by the number of dimensions. On a 100-dimensional problem, this mcans that the trust region is only shrunk after 100 consocutive failures, half of our evaluation budget. To increase exploitation, this tolerance is set to 10. Other baselines, such as the popular CMA-ES (Hansen, 2016), are described and shown in Appendix D.
Performance metric We report the normalized improvement over initial conditions (Hoffman et al., 2011),
s(no)) best (aee)-((an f(x)
[(x ((oo(nottt)
NT,
n= (O.....N.. NN:
(16)
where 16 is the number of initial obscrvations, f (mo) is the lowest initial value, (no is the incumbent solution, x= is the global minimizer and N = 200 is the number of acquisitions. 3 Best
Experiments with larger evaluation budget can found in Appendix
8

---

10
QBranin,200
QBranin,50p
QBranin.10D
Herma Red Improvenss
SHOM
SwTR
11X5 Proscsedi
X+ (Prposed?
1-X6-TR (Preacsed)
0.5
60
33
100
320
340
360
360
20
93
530
320
340
1.0
380
300
65
30
100
320
140
180
380
300
S+OM
S+
(Poposed)
(+.A (Prppsedt
I+XA+TR (Propoced)
0--
Noma
Stybiinsk/7ong..00
StyblinskiTang.500
Stybiinski-ang..1000
0.8
(Popases)
0.0
Figure Performance on test functions with 20, 50 and 100 dimensions. Additional results are shown Appendix D (methods and test functions) and Appendix J (larger evaluation budget). Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range. Abbreviations: Stationary (S), Cylindrical (C) and Informative (I) covariances; Quadratic Mean (+QM); Origin (+X0) and Adaptive (+XA) anchors; acquisitions within Trust Region (+TR).
Unlike simple regret, this performance metric is invariant to initial conditions and to the range of f, always starting from 0 and only attaining the maximum value of 1 if global optimum is found (zero simple regret). Additionally, our analysis does not hinge on the final best values because these can be decciving. For instance, a method that rapidly achieves good, though suboptimal, solution should not be considered inferior to another method that requires more evaluations, only to attain a marginally better value. The former is more sample efficient, which a desirable feature when functions are expensive evaluate. Some results are summarized by reporting the mean NI over N acquisitions (or, equivalently, the normalized area under the NI curve), with 0 indicating worst. and optimal performance. Respective optimization trajectories are shown in Appendix D.
5.2 Results Test Functions
We first focus on objectives whose landscapes are well-defined and can be easily controlled, allowing to examine whether more informative priors translate into higher sample efficiency. Our test functions are roughly bowl-shaped, making the quadratic functions, such Branin and Rosenbrock, have been previously used to test C (Oh et al., 2018). However, the Branin function is known for having multiple global optima, some of which are located close to the boundary. To overcome this, we design the QBranin function, which is obtained by adding a quadratic term, resulting in bowl-shaped function with single global optimum. The Styblinsky-Tane function, while multimodal, has unique global optimum. A more detailed characterization of these test functions, as well additional benchmarks (Levy and shifted objectives), can be found in Appendix B.3.
prior appropriate choice for these problems. Test

---

Table 1: Mean normalized improvement (over 200 acquisitions) on shifted (S35 and S50) Rosenbrock func-
tions (mean = standard deviation). Values 0.5 suggest superlincar improvement rates. Additional mothods, lest functions and optimization trajectories can be found in Appendix D. Abbreviations: Stationary (S): Cylindrical (c) and Informative (I) covariances; Quadratic Mean (+QM); Origin (+X0) and Adaptive (+XA)
anchors; acquisitions within Trust Region (+TR).
S35Rosenbrock
S50Rosenbrock
Mothod
50D
100D
50D
100D
0.503 + 0.078
0.283 + 0.052
0.690 0.054
0.411 0.048
S+QM
0.507 + 0.086
0.218 = 0.083
0.548 = 0.193
0.104 = 0.157
S+TR
0.677 + 0.030
0.486 + 0.030
0.752 + 0.020
0.535 0.061
C
0.744 + 0.077
0.487 = 0.067
0.700 0.118
0.377 = 0.142
I+X0
0.803 + 0.031
0.644 0.026
0.831 + 0.038
0.729 + 0.027
I+XA
0.803 + 0.016
0.659 + 0.023
0.857 + 0.030
0.726 + 0.039
I+XA+TR
0.800 + 0.038
0.693 + 0.016
0.844 + 0.027
0.720 + 0.031
Fixed anchor The proposed covariance functions are initially tested with a fixed anchor placed at the center (I+X0). Figure shows the performance curves on two objectives (top and middle) where the optimal solution is relatively close to this anchor. These results indicate that second-order nonstationary approaches, such C and I+xo, can significantly outperform stationary methods, S and S+QM, particularly in higher-dimensional domains. While informative mean (S+QM) is not necessarily better than an uninformative constant mean (8), prior information introduced via the covariance function leads to a consistent improvement over S. The reason may be that, in certain dimensions, S+QM has relatively large quadratic weights that hinder exploration away from the center, and in other dimensions, very small weights are unable to alleviate the boundary issue. Compared to S, S+QM and C, the proposed I+XO not only performs better on average but is also more reliable.
In Appendix E, an ablation study provides empirical evidence that incorporating both spatially-varying prior variances and lengthscales is important for tackling nigher-dimensional problems. In particular, there observe that the shorter hengthscaaess near anchors promote local exploration, which is especially effective when anchors are relatively close to optima. Meanwhile, spatially-varying prior variances can mitigate the boundary issue in high-dimensional domains, thereby improving sample efficiency. Next, we discuss the use of an adaptive anchor and the performance on the Styblinsky-Tang function.
Adaptive anchor The assumption that optimal solutions are close to the center is strong. As optima deviate from the center, the performance of C decreases markedly, eventually becoming comparable to that of S (Table 1). Despite the lower apparent vulnerability of I+XO to anchor misspecifications, this motivates the use of informative covariance functions with an adaptive anchor (I+XA). For simplicity, we opt for a greedy strategy where the anchor is given by the incumbent solution. For noise-free objectives, this is the point in the evidence set with the lowest observed value.
In addition to QBranin and Rosenbrock, Figure 5 includes the performance on the Styblinski- Sannin function, whose optimal solution is relatively distant from the center. Notably, the center is a local maximum, being surrounded by exponentially many local modes. As the assumption that the global minimum is close to the center is no longer valid, we observe that the relative performance of C and I+xo, compared to S and S+QM, deteriorates. In particular, C becomes the worst-performing method and I+X0 is worse than S in the 20-dimensional problem. Conversely, I+XA is consistently among the best-performing methods across all tests shown in the figure, and often outperforms the fixed-anchor variant (I+X0).
Trust region The use of an adaptive greedily-chosen anchor is partly motivated by trust-region opti- mization, wherein regions are centered at the incumbent solution. However, the proposed methodology is complementary to trust regions because it extends a stationary covariance function with spatially-varying parameters that are possibly aligned with prior beliefs about the objective. While S+TR can lead to signili- cant improvements over S, overexploration of boundaries in a trust region can still occur because the model is stationary. The method I+XA+TR shows that the combination of informative covariance functions with trust-region optimization can improve sample efficiency.
10

---

RoverT
RoverT,60D
1.0
Efficient Trajectory
Normalized Improvement
1+20 (Popened) (Propoed)
1+ (roposed
0.6
0.4
0.2
0.0
20 40
100 160
180
Figure Map layout of the rover trajectory problem (left) and respective optimization results (right). Additional results and more detailed analysis can be found in Appendix G. Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range. Abbreviations: Stationary (S), Cylindrical (C) and Informative (I) covariances; Quadratic Mean (+QM); Origin (+X0) and Adaptive (+XA) anchors; acquisitions within Trust Region (+TR).
MNISTW,10OD
MNISTW,100D
0.4
0.4
1+X0-- (Freposed)
(Prooeed) Proposell INt trippssll
I-Xiood Propssedd
0.0
0.0
100
120
140
160
180
200
0
20
40
60
80
100
120
140
160
180
200
Figure Bayesian optimization of 100 paramcters of a noural network trained and tested on MNIST.
Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range. Abbreviations: Stationary (S), Cylindrical (C) and Informative (I) covariances; Quadratic Mean (+QM); Origin (+x0), Good (+XGood) and Adaptive (+XA) anchors; acquisitions within Trust Region (+TR).
5.3 Applications
Rover trajectory planning Recent studies have used the rover trajectory (RoverT) problem to showcase BO (Wang et al., 2018; Eriksson et al., 2019; Eriksson & Jankowiak, 2021). The goal to optimize a trajectory such that the rover satisfies the target endpoints, Estart and Wend; while avoiding collisions with objects, as depicted in Figure 6 (left). The trajectory is given by B-spline that is fitted to 30 2-dimensional points, resulting in 60-dimensional optimization problem. A more detailed description and analysis is provided in Appendix G. As we show in Figure 6 (right), I+X0 and I+XA, which best-performing methods. We also observe that the performance of S+TR and I+XA+TR is significantly worse than that of S and I+XA, respectively. The reason may be that points in a Irust region do not significantly improve upon the incumbent solution. For optimal performance, the endpoints must match the two targets because deviation incurs significant loss. However, this may not be feasible when using small box centered a suboptimal solution, e.g. Øbest Consecutive failures decrease the size of trust regions, potentially resulting in even slower progress.
informative covariance, are the
11

---

Neural network optimization In this application, optimize a neural network layer, proposed by Oh et al. (2018). In particular, we perform BO of 100 parameters of a two-layer fully-connected network trained and tested on the MNIST dataset (TeCun, 1998). The architecture is 784 Wi.by, Nhidden Wa.ba, 10, with weights W, biases b, NNhidden = 10 and ReLU activations. The model is trained for 10 epochs with a batch size of 512 and evaluated on the test set using the negative log-likelihood loss. The weights We are found by BO using the test loss, while the remaining parameters are learned on the training set with Adam (Kingma & Ba, 2014). As pointed by Oh et al. (2018), the goal not to find solutions that generalize well, but rather to evaluate the ability of BO to find solutions that perform well on the test set by minimizing the respective loss. The landscape is multimodal, and points at the boundary can yield relatively good solutions.
Results are shown in Figure 7. While S appears to perform slightly better than I+XO and I+XA, models with an informative covariance function do not preclude the use of an informative mean function. For instance, we observe that I+XA+QM can improve upon I+XA toward the end, achieving similar average performance to that of S+QM. Relatedly, in Appendix D, we show that I+XA+QM outperforms both S+QM and I+XA on Styblinsky-Tang demonstrating that informative mean and covariance functions are complementary. The performefonn method is I+XGood, which assumes strong prior information by using fixed anchor placed at solution found by S+QM. This highlights the value of high-quality information and that it can be incorporated into informative GP priors, without modifying the acquisition process.
5.4 Limitations
Despite the overall superior performance shown by methods with informative covariance functions, there is room for improvement that is left for future work. For instance, as dimensionality increases, the weighted Fuclidean distance cannot. differentiate points that are far only in a few dimensions. This can lead to dif- ficulties in addressing the boundary issue in hypercubes. Potential remedies for this problem may include different parametrizations and distances based on the max-norm. Moreover, there is no mechanism in place that prevents greedily-chosee anchors from being close to the boundary. Since I+XA relies on empirical priors, it can be used without access to reliable prior information. However, as shown in Section 5.3 and Appendix D, methods that include reliable information can significantly outperform I+XA.
In Section 4, we provide formulas for covariance functions with multiple anchor points to demonstrate that the proposed methodology can in thcory capture knowledge about the objective that gocs beyond a single anchor. In our experiments, we focus on methods with a single anchor for several reasons. First, we wanted to show that higher sample efficiency is possible with simple design choices. Furthermore, baselines such as S+QM, S+TR and C implicitly assume the existence of a single anchor. By restricting ourselves to this case, the comparisons are fair, and we can assess the robustness of the different methods to anchor misspecification.
Additionally, in sequential myopic BO, which is the setup considered in this paper, it suffices to focus region that includes optimum. The correct identification of such region may be difficult in multimodal objectives, but as we have shown, the anchor can be dynamically adjusted. One option, without having to specify simultaneous anchors, is to maintain posterior beliefs about promising points and sample them at each step. the anchors are not known, a similar strategy to that proposed by Eriksson et al. (2019) to maintain multiple trust. regions is also a possibility. Other options include setting the anchors after performing cluster analysis and learning them by empirical Bayes. In this context, a different set of baselines would also be more appropriate, including methods with multiple trust regions and quadratic mean functions with multiple centers. Moreover, model selection can be used to automatically determine the number and location of anchors (see e.g. Appendix H). A careful treatment of informative covariance functions with multiple anchors requires substantial analysis and testing, and we believe this is best left for future work.
6 Related Work
Second-order nonstationarity Covariance functions with spatially-varying lengthscales were initially ex- plored by Gibbs (1998) and Paciorek & Schervish (2006), and augmented with spatially-varying prior (and noisc) variance by Heinonen et al. (2016) for GP regression. These models are derived from basis function expansions or process convolution (Risser, 2016), wherein the covariance function is obtained via integral transform of base stochastic process with kernel function, which may also depend on auxiliary mixing
12

---

variables (e.g. Wang et al., 2020a). In contrast, we introduce spatially-warying lengthscales via imput warping, resulting in a simpler expression that does not require extra factors to ensure positive definiteness. Addi- tionally, these other models place GP priors on the spatially-varying parameters and infer the respective latent functions by MCMC, being computationally expensive for high-dimensional parameter spaces. In- stead, we use pointwise specifications of nonstationary effects, governed by semiparametric functions whose shape and parameters can be clicited from an expert or Icarned by empirical Bayes. Importantly, while these works focus on the estimation of the covariance function for regression purposes, we leverage second-order nonstationarity to incorporate information about promising points for optimization.
In the context of optimization, input warping was first studied by Snock et al. (2014), focusing on monotonic input deformations via CDFs, specifically the Beta distribution. The motivation follows from a regression perspective, as the resulting GP is capable of approximating more realistic objective functions. More recently, Oh et al. (2018) proposed cylindrical covariance functions for optimization, outperforming the method by Snock ct al. (2014). In this case, after applying the cylindrical transformation, the radius warped by the Kumaraswamy distribution, further expanding the innermost, region of the search space, which is assumed to contain the optimal solution. A different methodology was explored by Martinez-Cantee used lincar combinations of local and global stationary covariance functions with input-dependent weights to induce spatially-varring lengthscales (piecewise constant). While the motivation is that the local assignment of shorter lengthscales promotes local exploration, the local component is not guaranteed to learn shorter lengthscales because both are initialized with the same hyperpriors. None of these works provide a regret analysis or consider spatially-varying prior (co)variances in their models.
(2015), who
High-dimensional Bayesian optimization To tackle the challenges presented by high-dimensional spaces, previous research has relied on certain structural assumptions (Binois & Wycoff, 2022), namely low effective dimensionality (e.g. Letham et al., 2020; Eriksson & Jankowiak, 2021) and additivity (e.g. Kandasamy et al., 2015; Gardner et al., 2017; Rolland et al., 2018). saauce-prtitioning schemes and trust regions have also been proposed (c.g. McLood et al., 2018; Eriksson et al., 2019; Wang et al., 2020b; Dionane et al., 2021). Notably, these methods employ stationary covariance functions. Our work complements these by introducing more informative priors (see e.g. Appendices D and H). Cylindrical covariance functions (Oh et al., 2018) can also be augmented with spatially-varying (co)variances and anchors can be defined in the transformed space.
Priors over the optimum There have been limited attempls to incorporate beliefs about optimal solu- tions into BO. Souza et al. (2021) combined predefined prior distributions with the tree-structured Parzen estimator approach by Bergstra et al. (2011), which estimates the input given the output, as opposed to the standard GP approach of modeling the output given the input. By design, the former is restricted to the Ex- pected Improvement acquisition function. Ramachandran et al. (2020) explored input warping via predefined CDFs that encode prior beliefs. More recently, Hvarfner al. (2022) proposed belief-augmented acquisition functions. As we show in Appendix I, these functions and informative GP priors are complementary.
7 Conclusion
This work proposes informative covariance functions for Bayesian optimization, leveraging nonstationarity to incorporate input-dependent information. Our methodology relies on priors with spatially-varying paramcters to express a priori preferences for certain regions and to promote multiscale exploration. These priors can more efficiently describe wider class of objective functions and result in improved worst-case performance.
Our experiments showed that the proposed covariance functions can significantly increase sample efficiency, even under weak prior information, challenging the prevalent use of stat ionary models for optimization. Ad- ditionally, we showed that our work complements existing methodologies, including models with informative mean functions, trust-region optimization and belief-angiented acquisition functions. Despite our focus on continuous search spaces, we believe that this work can be extended to other types, such as discrete and mixed spaces, by adopting appropriate priors. Our methodology can also accommodate multiple anchors and be combined with other methods, relying for instance on additivity and low effective dimensionality, We recognize that priors over the possible locations of the optimum can be defined in transformed, possi- bly lower-dimensional spaces, but these details are left for future research. Overall, we see this work as a contribution toward the wider adoption of scalable nonstationary and informative models for optimization.
13

---

References
Luigi Acerbi. An Exploration of Acquisition and Mean Functions in Variational Bayesian Monte Carlo. In Symposium on Advances in Approcimale Bayesion Inference, pp. 10. PMLR, 2019.
Nachman Aronszajn. Theory of Reproducing Kernels. Transactions of the American Mathematical Sociely, 68(3):337-404. 1950.
John-Alexander M. Assael, Ziyu Wang, Bobak Shahriari, and Nando de Freitas. Heteroscedastic Treed Bayesian Optimisation. arXiv preprint arXiv:1410. 7172, 2014.
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 Newral Information Processing Systems, 2111244-11335 2020.
James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl. Algorithms for Hyper-Parameter Opti- mization. Advances in Neural Information Processing Systems, 24, 2011.
Mickael Binois and Nathan Wycoff. A survey high-dimensional Gaussian process modeling with applica- tion to Bayesian oplimization. ACM Transactions on Evolulionary Learning and Oplimization, 2(2):1 26, 2022.
Adam D. Bull. Convergence Rates of Efficient Global Optimization Algorithms. Journal of Machine Learning Research, 12()0), 2011.
Sayak Ray Chowdhury and Aditya Gopalan. On Kernelized Multi-armed Bandits. In International Confer- ence on Machine Learning, pp. 844 143-833.77 PMLR, 2017.
Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences, 117(48):30055 -00662,, 2020.
Y. Diouane, V. Picheny, R. L. Riche, and Alexandre Scotto Di Perrotolo. TREGO: a Frust-Region Framework for Efficient Global Optimization. arXiv preprint arXiv:2101.06808 2021.
David Eriksson and Martin Jankowiak. High- liggh-iimesion Bayesian Optimization with Sparse Axis-Aligned Subspaces. In Uncertainly in Artificial Intelligence, pp. 493 503. PMLR, 2021.
David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable Global Optimization via. Local Bayesian Optimization. Advances in Neural Information Processing Systems, 32, 2019.
Peter I. Frazier. A Tutorial on Bayesian Optimization. ar preprint orXic:1807.02811 2018.
Jacob Gardner, Chuan Guo, Kilian Weinberger, Roman Garnett, and Roger Grosse. Discovering and Exploit- ing Additive Structure for Bayesian Optimization. In International Conference Artificial Intelligence and Statistics, pp. 1319. PMLR, 2017.
Jacob Gardner, Geoff Pleiss, Kilian Q. Weinberger, David Bindel, and Andrew G. Wilson, GPyTorch: Black- box Matrix-Matrix Gaussian Process Inference with GPU Accelcration. Advances in Neural Information Processing Systems, 31, 2018.
Roman Garnett, Michael A. Osborne, and Philipp Hennig. Active Learning of Linear Embeddings for Gaussian Processos. In Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, pp. 230 239, 2014.
Marc Genton. Classes of Kernels for Machine Learning: Statistics Perspective. JMLR, 2:299 312, 2001.
Mark N. Gibbs. Bayesian Goussion processes for Regression and Classification. PhD thesis, University of Cambridge, 1998.
14

---

Rafael
Jennifer N. Wei, David Duvenaud, José Miguel Hernández-Lobato, Benjamín Sánchez-Lengee Lengelin, Dennis Sheberla, Jorge Aguilena-Iparraguiee Timothy D. Hirzel, Ryan P. Adams, and Alán Aspuru-Guzik, Automatic Chemical Design Using a Data-Driven Continuous Representation of Molecules. ACS Central Science, 4(2):268-2)) 276, 2018.
Stewart Greenhill, Santu Rana, Sunil Gupta, Pratibha Vellanki, and Svetha Venkatesh. Bayesian Optimiza- tion for Adaplive Experimental Design: A Review. IEEE Access, 8:13937-1333 2020.
Antoine Grosnit, Rasul Tutunov, Alexandre Max Maraval, Ryan-Rhys Griffiths, Alexander I. Cowen-Rivers, Lin Yang, Lin Zhu, Wenlong Lyu, Zhitang Chen, Jun Wang, al. High-Dinensional Bayesian Opti- misation with Variational Autoencoders and Deep Metrie Learning. arXiv preprint arXiv:2106.0800 036600 2021.
Sunil Gupta, Santu Rana, Svetha Venkatcsh, al. Regret Bounds for Expected Improvement Algorithms in Gaussian Process Bandit Optimization. In International Conference on Artificial Intelligence and Slatistics, pp- 8715-8737. PMLR, 2022.
Michael U. Gutmann, Jukka Corander, et al. Bayesian Optimization for Likelihood-Fiee Inference of Simulator-Based Statistical Models. JMLR, 2016.
Nikolaus Hansen. The CMA Evolution Strategy: A Tutorial. arXiv preprint arXiv: arXin:164. 004.00044077 2016.
Nikolaus Hansen, Youhei Akimoto, and Petr Baudis. CMA-ES/pyema, Zenodo, DOI: 10.5281/zen- 2019.
Markus Heinonen, Henrik Mannerström, Juho Rousu, Samuel Kaski, and Harri Lähdesmäki. Non-Stationary Gaussian Process Regression with Hamiltonian Monte Carlo. In Artificial Intelligence and Statistics, pp. 732-740. PMLR, 2016.
James Hensman, Nicolo Fusi, and Neil D. Lawrence. Gaussian Processes for Big Data. Proceedings of the Twendy-Ninth Conference on Uncertainly in Arlificial Intelligence, 2013.
Matthew Hoffman, Eric Brochu, Nando de Freitas, et al. Portfolio Allocation for Bayesian Optimization. In 2011.
UAT, pp.
Carl Hvarfner, Danny Stoll, Artur Souza, Marius Lindauer, Frank Hutter, and Luigi Nardi. wBO: Augment- ing Acquisition Functions with User Beliefs for Bayesian Optimization. In International Conference on Learning Representations, 2022.
Marko Järvenpää, Michael U. Gutmann, Aki Vehtari, and Pekka Marttinen. Parallel Gaussian Process Surrogate Bayesian Inference with Noisy Likelihood Evaluations. Bayesion Analysis, 16(1):147 178, 2021.
Stephen Joe and Frances Y. Kuo. Constructing Sobol Sequences with Better Two-Dimensional Projections.
SIAM Journal on Scientific Computing, 30(5):2635-265) 2008.
Donald R. Jones, Matthias Schonlau, and William J. Welch. Efficient Global Optimization of Expensive Black- Box Functions. Journal of Global optimization, 13(4):455-45) 492, 1998.
M.C. Jones, Kumaraswamy's distribution: beta-type distribution with some tractability advantages. Statistical Methodology, 6(1):70-8) 81, 2009.
Motonobu Kanagawa, Philipp Hemnig, Dino Sejdinovic, and Bharath K. Sriperumbudur Gaussian Processes and Kernel Methods: A Review on Connections and Equivalences. arXiu preprint arXiv:1807.0m 022222 2018.
Kirthevasan Kandasamy, Jeff Schneider, and Barnabás Poczós. High Dimensional Bayesian Optimisation and Bandits via Additive Models. In International Conference Machine Learning, pp. 295 304. PMLR, 2015.
Xiv preprint Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. 142.6612.6998 2014.
15

---

Harold J. Kushner. A Versatile Stochastic Model of Function of Unknown and Time Varying Form. Journal of Mathematical Analysis and Applications, 5(1):150- 1600001661111 1962.
Yann LeCun. The MNIST database of handwritten digits. http:///wnnneeuuncco haamm/erdb/mmsst// 1998.
Ben Letham, Roberto Calandra, Akshara Rai, and Eytan Bakshy. Re- ReeExxmining Linear Embeddings for High-Dimension Bayesian Optimization. Advances in Neural Information Processing Systems, 33, 2020.
Chun-Liang Li, Kirthevasan Kandasamy, Barnabás Poczós, and Jeff Schneider. High Dimensional Bayesian Optimization via. Restricted Projection Pursuit Models. Artificial Intelligence and Statistics, pp. 884- 892. PMLR, 2016.
Ruben Mantine Local Nonstationarity for Efficient Bayesian Optimization. arXiv preprint arXiv:1506. 02080, 2015.
Mark McLeod, Stephen Roberts, and Michael A. Osborne. Optimization, fast and slow: optimally switching between local and Bayesian optimization. In International Conference on Machine Learning, pp. 3443- 3452. PMLR, 2018.
Erich Merrill, Alan Fern, Xiaoli Z. Fern, and Nima Dolatnia. An Empirical Study of Bayesian Optimization: Acquisition Versus Partition. JMLR, 22:4-1, 2021.
Charles A. Micchelli, Yuesheng Xu, and Haizhang Zhang. Universal Kernels. JMLR, 7(Dec):2651-266)) 2006.
Jonas Močkus. On Bayesian Methods for Sceking the Extremum. Optimization Techniques IFIP Technical Conference, pp. 400-404 404.04..55 Springer, 1975.
Riccardo Moriconi, Marc Peter Deisenroth, and K. S. Sesh Kumar. High-dimensional Bayesian optimization using low-dimensional feature spaces. Machine Learning, 109(9):1925 1943, 2020.
Mojmír Mutný and Andreas Krause. Efficient High Dimensional Bayesian Optimization with Additivity and Quadrature Fourier Features. Advances in Neurol Information Processing Systems, 31, 2018.
Willie Neiswanger, Lantao Yu, Shengjia Zhao, Chenlin Meng, and Stefano Ermon. Generalizing Bayesian Optimization with Decision-theorctic Entropics. Advances in Newral Information Processing Systems, 35, 2022.
ChangYong Oh, Efstratios Gavves, and Max Welling. BOCK: Bayesian Optimization with Cylindrical Kernels, In International Conference on Machine Learning, pp. 3868 3877. PMLR, 2018.
Art B. Owen. Scrambling Sobol and Niederreiter- Points. Journal of Complexity, 14(4):466 489, 1998.
Christopher J. Paciorek and Mark J. Schervish. Spatial Modelling Using New Class of Nonstationary Covariance Functions. Ensirommetrics, 17(6):483-506, 2006.
Anil Ramachandran, Sunil Gupta, Santu Rana, Cheng Li, and Svetha Venkatesh. Incorporating Expert Prior in Bayesian Optimisation via Space Warping. Knowledge-Bass Based Systems, 195:105663, 2020.
Elena Raponi, Hao Wang, Mariusz Bujny, Simonetta Boria, and Carola Docrr. High Dimensional Bayesian Optimization Assisted by Principal Component Analysis. In International Conference on Parallel Problem Solving from Nature, pp. 189-183. Springer, 2020.
Carl E. Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. MIT Press, 2005.
Mark D. Risser. Nonstationary Spatial Modeling, with Emphasis on Process Convolution and Covariate- Driven Approaches. Xixiv preprint orXim:1610.02444 2016.
16

---

Paul Rolland, Jonathan Scarlett, Iijja Bogunovic, and Volkan Cevher. high-Dimensional Bayesian Optimiza- tion via Additive Models with Overlapping Groups. In International Conference on Artificial Intelligence and Statistics, pp. 298 307. PMLR, 2018.
Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P. Adams, and Nando De Freitas. Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proceedings of the IEEE, 104(1):148-175, 2016.
Jasper Snoek, Hugo Larochelle, and Ryan P. Adams. Practical Bayesian Optimization of Machine Learning Algorithms.. In Advances in Newral Information Processing Systems, pp. 2951 2959, 2012.
Jasper Snoek, Kevin Swersky, Rich Zemel, and Ryan Adams. Input Warping for Bayesian Optimization of Non-stationary Functions. In International Conference on Machine Learning, pp. 1674-1682. PMLR, 2014.
Jasper Snock, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Mostofa Patwary Mr Prabhat, and Ryan Adams. Scalable Bayesian Optimization Using Deep Neural Networks. PMLR, 2015.
In International Conference on Machine Learning, pp.
Artur Souza, Luigi Nardi, Leonardo B. Oliveira, Kunle Olukotun, Marius Lindauer, and Frank Hutter. Bayesian Optimization with Prior for the Optimum. In Joint European Conference Machine Learning and Knowledge Discovery in Databases, pp. 265 296. Springer, 2021.
Niranjan Srinivas, Andreas Krause, Sham Kakade, and Matthias Seeger. Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design. In Proceedings of the 27th International Conference International Conference Machine Learning, pp. 1015-1022, 2010.
Kevin Jordan Swersky, Improving Bayesian Optimization for Machine Learning using Expert Priors. PhD thesis, University of Toronto, 2017.
Michalis Titsias. Variational Learning of Inducing Variables in Sparse Gaussian Processes. In Artificial Intelligence and Statistics, pp. 574-574. PMLR, 2009.
Sattar Vakili, Kia Khezeli, and Victor Picheny. On Information Gain and Regret Bounds in Gaussian Process Bandits. In International Conference on Artificial Intelligence and Statisties, pp. 82-90. PMT.R, 2021.
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Ev- geni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C. J. Carey, ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Tan Henriksen, F. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fun- damental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272, 2020.
Kangrui Wang, Oliver Hameljjnck, Thcodoros Damoulas, and Mark Steel. Nonstationary Nonseparable Random Fields. In International Conference on Machine Learning, pp. 9887-9897. PMLR, 2020a.
Linnan Wang, Rodrigo Fonsoca, and Yuandong Tian. Learning Search Space Partition for Black-box Op- timization using Monte Carlo Tree Search. Advances in Neural Information Processing Systems, 33: 19511-19522 2020b.
Zi Wang, Clement Gehring, Pushmeet Kohli, and Stefanie Jegelka. Batched Large-scale Bayesian Optimiza- tion in High-dimensional Spaces. In International Conference Artificial Intelligence and Statistics, pp. 745 754. PMLR, 2018
Ziyu Wang and Nando de Freitas. Theoretical Analysis of Bayesian Optimisation with Unknown Gaussian Process Hyper-Parameters. arXie preprint urXiv:1406. 7758, 2014.
Ziyu Wang, Frank Hutter, Masrour Zoghi, David Matheson, and Nando de Feitas. Bayesian Optimization in Billion Dimensions via Random Embeddings. Journal of Artificial Intelligence Research, 55:361-387, 2016.
17

---

Colin White, Willie Neiswanger, and Yash Savani. BANANAS: Bayesian Optimization with Neural Archi- tectures for Neural Architecture Search. In Proceedings of the AAAT Conference on Artificial Intelligence, volume 35, pp. 10293 10301, 2021.
James Wilson, Frank Hutter, and Mare Deisenroth. Maximizing Acquisition Functions for Bayesian Opti- mization. Advances in Neural Information Processing Systems, 31, 2018.
Manuel Wüthrich, Bernhard Schölkopf, and Andreas Krause. Regret Bounds for Gaussian- Process Opti- mization Large Domains. Advances in Neural Information Processing Systems, 34, 2021.
Ciyou Zhu, Richard H. Byrd, Peihuang Lu, and Jorge Nocedal. Algorithm 778: L-BFGS-B: Fortran Subrou- lines for Large-Scale Gound-Constrainee Optimization. ACM Transactions on Mathemalical Software, 23 (4):550-566 1997.
18

---

A Instantaneous Regret Bounds
Lemma A.1. Assume the interval If(x) mn (x)] S (n)n(() holds, Yx E X and that global optima of acquisition functions can be found. Then, for popular acquisition functions, namely LCB and EI, and any coveriance function, possibly noustationary, the worsl-case inslanlaneous regrel depends on the width of the interval, being proportional to the posterior predictive standard deviation, Tn+1 on (n(nn))), where In+1 is the acquisition at step n 1.
Proof. Let x² be the minimizer of f(x). Then, based (Srinivas et al., 2010, Lemma 5.2), we have the following instantaneous regret bound for LCB,
Tn+1 f(xn+1) f(x²)
(17)
$ f(In+1) - mx(x'') + (n)n(("))
(18)
[(xm+1) - TCBB(x*)
(19)
S [(#n+1) - LCB((en11)
(20)
f(xn+1) mm((n++1 + Bn
(21)
S (2
(22)
Regarding EI, we know (x) S ET(x) (x:) (B,, 1) (a)), by (Wang & de Freitas, 2014, Lemma 9),
with B, = 7 S and r(z) =Fx(2> + W[2;; [0. Also, recall that the improvement is given by I, (x)) max(0, mi f(x)) with m, = minl m, (z). Now, we derive a simpler, tighter upper bound on EI,
0, 1). because is nondecrcasing, r'(z) = Fx E
El(x) = (z(2)((2(())
(23)
S =xx(x))(
(24)
2n on (x) nn (x),
(25)
where the first incquality follows from z(x)) (m, mn (x)) /a (z) S 0 and nondccrcasing т. Lot x, arg min mn(x). Then, the instant regret bound can be written
Tn+1
f(xn+1) - f(x²)
(26)
(f(#n+1) m, (m, f(xx))
(27)
$ (f(In+1) m, + (z=)
(28)
1
S (f(#n+1) m. Bn
El(x")
(29)
1
$ (f(#n+1) mm,))
EI [n+11)
(30)
Bn
I (f(In+1) mm,)) +
B, on (I
(31)
(mn( - mm)) +
Bn + Br
on (n((n)1)
(32)
1
Da(wn + +
on (x(4))
(33)
Bn
s Dt + + 11)
1
on (n((n))))
(34)
Bn
where, Equation 32, the difference bounded as Im,(# - m5l S D, (Fa+t))o (Ma+1), with
Da(( = / 2 log(o, we also know that on (x,,) can be lower bounded by a positive value, and that there is a positive constant upper bounding o, (I by construction. Thus, for any covariance function, possibly nonstationary, there exists positive constant D,, such that D 2 D (m,,1)) which yields the result.
by (Wang & de Freitas, 2014, Lemma 10). From the same lemma,
19

---

B Implementation Details
B.1 Gaussian Processes
We used GPy'Torch (Gardner et al., 2018) to implement all GP models in Python, including the proposed informative (I) covariance functions, cylindrical (c) covariance functions (Oh et al., 2018) and axis- alisned quadratic mean (+QM) (Snoek et al., 2015). Unlike the version provided by GPyTorch, our implementation of C follows the original training and prediction routines, which account for the special treatment of the origin, as discussed in more detail in Appendix C.
In order to ensure fair comparisons, all models are optimized by L-BFGS-B (Zhu et al., 1997) with a maximum of 1000 iterations per new acquisition and other options set to default values. In addition, the stationary components Cs are Matérn (v = 5/2), equipped with the weighted Fuclidean distance (lengthscales Xa). Positivity constraints are handled internally by GPyTTrcc using the softplus transform, except for C-specific hyperparameterss which use the log transform (Appendix C). General hyperparameters are given uninfor- mative priors with bounds from (Oh et al., 2018). These include the prior variance of ~ =L(ee 12ee-1,,aa0 e2, and lengthseales Xa ~ U(e-12 12 2/D), de (1...D) where the upper bound corresponds to the maximum Eu- clidean distance in the centered hypercube X = [--, Table 2 summariizes the default settings of I. For noise-free objectives, the noise hyperparameter is set fixed value of =0-3
By default, models characterized by uninformative prior mean function with a uniformlyydistributted constant b (ymmin Imax); where the bounds correspond to the minimum and maximum observed values. Variants +QM use an axis-aligned quadratic mean function m(x) b Es (dd(())) [co]a)²; as suggested by Snoek et al. (2015) and Swersky (2017, Section 4.4). The quadratic weights are independent, distributed as ad ~ 1S * (2). The half-Horseshoe is a anike-and-slab hyperprior with a spike at and a slab on the positive real line, only allowing convex quadratic functions. Quadratic mean functions with adaptive, greedily-chosen centers Cp- = Destt were tested, but no significant or consistent improvement compared to the default center Co 0 was observed. Hence, results using the former (+QM+XA) are only shown in Appendix D.
Table 2: Default settings of informative covariance funct ions.
Quantity
Setting
Stationary corr. Cs
Matérn (v (p= 5/2, oă = 1)
Prior variance of
u(e (2e=12 u(e-12,2v1 2V
Lengthscales Xa
Distances d,, dx
Tied do
Distance do
Weighted Euclidean dwE (shared Aa)
Kernels Ka, kx
Tied kn Gaussian
Kernel ko
Ratios To, TA
Tied ro
Ratio ro
Kumar(3.164, 1000),
Sensitivity analysis (Appendix F)
B.2 Acquisition
Prior to acquiring new data, an initial evidence is formed. The set Dno includes the origin and 15 other pscudo-uniformly distributed points, drawn according to scrambled Sobol sequence (Owen, 1998; Joe & Kuo, 2008). This set provides the initial conditions from which initial models are estimated and their performance measured. For noise-free objectives, the incumbent solution is the point with the lowest observed (mg+n) = arg mingeD., R f(x), nig 16, n = (O...., N), N = 200.
value a hest
The default acqisition function is the Expected Improvement (EI), which is implemented in BoTorch (Ba- landat et al., 2020). Maximization performed on the CPU via gradient- based optimization with multiple restarts. The 20 initial points for the optimization selected from 20010 candidates, including 10 points
20

---

that are chosen heuristically and 20000 points that are randomly generated. The former are obtained by small Gaussian perturbations of the incumbent solution, and the latter are drawn from a scrambled Sobol sequence. The main difference between this procedure and the one implemented by Oh et al. (2018) is the use of L-BFGS-B via BoTorch, as opposed to Adam (Kingma & Ba, 2014).
B.3 Test Functions
Test functions are set up in a way similar to that used by Oh et al. (2018), with their domains adjusted to the centered hypercube x = [| by applying a linear Iranslormation the inpuls. Additionally, original objectives are shifted and scaled satisfy f(x*) = 0 and f(0) 100. Since the resulting functions are non- negative, models are then estimated on log-transformed data, after adding very small offsct. In general, this leads to improved performance because the distribution of function values becomes more Gaussian-like. Performance metries are however computed without the log transform.
QBranin The Branin function originally evaluated on [-5, 10] x [o, 15), and has several global optima. This function is modified by adding a quadratic component (last term), transforming it into a bowl-shaped function with only one global optimum,
(x) [z)2
5.1
+ [ 6)
2
+ 10 (1 10(1-87)
1 8T
cos([2>]] cos( [a>)]) + + 5|0]
(35)
47
Hligher-dimession versions are obtained by additive repetition.
Shifted QBranin The shifted variants SQBranin and SSQBranin are derived from QBranin by shifting the inputs in the original space by and respectively. Additional experiments with these variants are included in Appendix D. Figure 8 provides a 2D illustration of these functions in the transformed space [-1,
Rosenbrock This function is originally evaluated on the hypercube [-5, 10 and is given by
D-1
f(x) = Σ [100[[e]<+1- - [x]3)² + ([x]a 1)²].
(36)
Similar QBranin, the global optimum close the origin in the transformed space, x" -0.2 x However, shown in Figure 9, the Rosenbrock is characterized by narrow banana-shaped valleys, making global optimization more challenging.
Shifted Rosenbrock S35Rosenbrock, S50Rosenbrock, and S65Rosenbrock are three variants of the Rosen- brock function, designed to test the robustness of methods that assume good values near the origin. These functions are obtained by shifting the inputs so that the global minima the transformed space are at 0.35 1, 0.50 x 1, and x² 0.65 1. Figure provides a 2D illustration of these functions.
Levy The original domain is 1010, 10 and the function is defined
((x) = sin sin?(ww))
Σ
(wa- 1)² [1 10 (Twa + 1)] + (w; 1)² [1 + sin 3inw()))
(37)
d=1
Wa [x)
ie DJ.
(38)
4
As shown in Figure 10 (left), the Levy function is characterized by strong sinusoidal components. The global optimum again close to the origin, specifically at x* 0.1 x in the transformed space. Styblinski- Tang The original domain is [-5, 5]" and the function is given by
1 Σ ([=) 16[=]3 +(x]))))
f(x)
(39)
Joint optimization of the Styblinski- Maanisski-Tnngg function is difficult because has exponentially many local modes 2D 1. In contrast to QBranin, Rosenbrock and Levy, the global optimum is relatively far from the origin, -0.88 in the transformed space. Figure 10 (right) provides an illustration of this multimodal function.
21

---

QBranin
SQBranin
SSQBranin
1.00
1.00
1.00
Global Minimum ([-0.333 0.2 1)
Global Minimum (1-0.6 -. 0.467])
Global Minimum (]-0.733 0.6
0.75
Origin
0.75
Origin
0.75
Origin
0.25
0.25
0.25
0.00
0.00
0,00
-0.50
-0.50
-0.50
-0.75
-0.75
-0,75
-1.09
-0.75
-100
-1.0
Figure A 2D slice of QBranin and shifted variants. These test functions are unimodal and bowl-shaped. The purpose of SQBranin and SSQBranin to test the robustness of methods that assume good values near the origin.
Rosenbrock
$35 Rosenbrokk
1.00
1.00
Global Miniuum ([-0.2 -0.21)
Global Minimum ([0.35 0.35] Origin
Origin
0.50
0.50
0.25
0.00
-025
-0.25
0.75
-0.75
-1.09 .000
0.75
0.50
0.25
0.00
0.25
0.50
0.75
1.00
-1.00
1.00
0.75
0.58
0.25
0.00
0.25
0.50
0.75
1.00
S50Rosenbrock
S65Rosenbrock
1.00
1.00
Global Minimum ([0.5 0.5))
Global Minimum ([0.65 0.65)
Origin
Origin
0.50
0.50
0.25
0.00
0.00
0.25
0.25
-075
-0.75
-1.00
0.75 0.50
0.25
0.00
0.25
0.50
0.75
1.00
-1.00
1.00
0.75
0.50
0.25
0.00
0.25
0.50
0.75
1.00
Figure 9: A 2D slice of Rosenbrock and three shifted variants (S35, $50 and S65).
Levy
Styblinski-Tang
Global Minimum (10.1 D.1])
Global Minimum ([-0.581 -0.581)
0.75
Origin
0.75
Origin
0.50
0.25
0.25
0.00
0.00
-0.50
-0.50
-0.75
-0.75
-1.00
.00
-0.75
-0.25
-1.00
-0.5
Figure 10: A 2D slice of Levy and Styblinski-Tangg
22

---

C Cylindrical Covariance Functions
Cylindrical covariance functions (Oh et al., 2018) aim to address the boundary issue that arises in high- dimensional Bayesian optimization. The cylindrical transformation T(x) represents a point x E X in terms of its radius r and angular components a:
(r. - T(x) 1 {{ (11rl12, 2/|[=]1))) if | + 0 (0, Garbitrary). if | =0
(40)
where Garbitrary is a random unit vector. Geometrically, this transformation maps balls of radius R onto the surface of a cylinder of height R. A spherical shell of width dr centered at a point a corresponds to a region in the Euclidean space whose volume increases exponentially with the radius r. As result, algorithms that aim equally covering each volume element in the Euclidean space spend exponentially more time at the boundary of the search space. In contrast, Oh et al. (2018) proposed to spend an equal amount of resources each volume element in the transformed space. Intuitively, the transformation leads to search in the Euclidean space that "expands the region near the center while contracting the regions near the boundary" (Oh et al., 2018, Section 3.2). The assumption is that optimal values are more likely to be found near the origin, so it is beneficial to encourage exploration in this region.
Remarkably, this transformation poses a challenge when comparing any non-origin point to the origin, since the latter is represented by an infinite set of points with radius 0. Oh et al. (2018) proposed using the point in the set that is closest to the point under comparison æ, setting aabiitrary x ////=||2- However, this solution has significant computational implications, requiring custom inference routines. If the origin is included in the training set, the Gram matrix must be recomputed for each test point. This issue is mitigated by block matrix inversion, where the block containing all non-origin points can be procomputed and reused.
Cylindrical covariance functions Ceyl are defined by the product of a 1-dimensional Matérn (v = 5/2), measuring similarity of radii, and function C2 that compares angular components,
3
Cala1, a2) = Σ c((aaa)))
(41)
p=0
where Cp is normalized non-negative weight, given a log-Norma prior log(e) ~ (0,, 2²). To further en- courage the expansion of the innermost region, Oh et al. (2018) used input warping on the radius com- ponent. This involves transforming the radius according to the Kunaraswany cumulative distribution, FKuma(1) (r) = 1- (1 ra)a, with E [0.5, 1] and 3 E [1,2]. Both and 8 are given spike-and-slab hyperprior on the log scale, with spike at log(1). Figure 11 shows possible warpings, including the identily transform obtained with a 1 and B = 1.
0.8
0.6
a-0.25, f-1
0.4
aw0.25, p-2
a=0 25, p-4
2-0.50,
2=0.50, f=2
0.2
==0.50.$=4
a=1 00, pm1
00,
0.0
am1 00, P
0.2
0.6
1.0
Figure 11: Input warping the radius. Concave transformations expand regions of small radii.
23

---

D Additional Results
D.1 Main Methods
S Bayesian Optimization (BO) with an uninformative Gaussian process (GP) prior, characterized by constant mean and a stationary Matérn covariance function.
S+QM s with an axis-aligned quadratic mean function (Snoek et al., 2015; Swersky, 2017).
S+TR S with acquisitions within an adaplive trust. region, i.e., a box centered at the incumbent solution that is farunk/expannee the incumbent solution is the point with the lowest observed value, as described in Appendix B.2.
based on consecutive failures/succosses (Friksson et al., 2019). For frese-free objectives,
C BO with GP characterized by a constant mean and cylindrical covariance function (Oh et al., 2018).
I+X0 (Proposed) BO with a GP featuring a constant mean and an informative covariance function. A single, fixed anchor is placed at the center (origin) of the search space.
I+XA (Proposed) I+XO with an adaptive greedily-chosen anchor, given by the incumbent solution. This is the same point that used in S+TR.
1+XA+TR (Proposed) I+XA with a trust region, using the same adaptation scheme as S+TR.
D.2 Additional Methods
S+QM+XA Quadratic mean function with greedily-chosen center, given by the incumbent solution. In most tests, including the shifted objectives in Figures 12 and 14, S+QM+XA was outperformed by S+QM.
GS+XA Informative covariance functions are based on a search model for the optimum, which is given by a mixture distribution (Equation 11). For single anchor, Gaussian kernel and weighted Euclidean distance, the equivalent distribution a mixture of a uniform and a Gaussian distribution with diagonal covariance matrix. In GS+XA, the mean is set the incumbent solution and the variance is estimated from the data, The standard deviation is the median absolute deviation (MAD) of the best 5% training points (rounded, with minimum of 2), assuming that the coordinates are independently and identically distributed. The Gaussian mixing weight is set to 0.9. Other combinations involving a smaller mixing weight (0.5), a larger training set (10%) and a scale correction of the MAD estimate were also tested, but overall performed no better than CS+XA. In turn, GS+XA typically performed worse than S or C, as shown in Figures 12, 13 and 14.
CMA-ES Covariance Matrix Adaptation Evolution Strategy (Hansen, 2016) is a popular global optimization algorithm and is related to GS+XA, in that it samples from a Gaussian whose mean and covariance matrix are adaptive. In terms of implementation, In general, CMA-ES was outperformed by S, as shown in Figures 12 and 14. On Styblinski-Taun (Figure 13), however, it performed well, but no better than I+XA or I+XA+QM.
use pyema (Hansen et al., 2019) with a population size of 10.
I+XA+QM (Proposed) The uninformative constant mean I+XA is replaced by the quadratic mean from S+QM. In tests where S+QM performed comparatively well, e.g. MyylliskiiTaang (Figure 13), the combination I+XA+QM proved to be effective.
I+X++F (Proposed) I+XA uses empirical priors, optimizing the hyperparameters via marginal likelihood. However, if additional information about the objective is available, better performance can be achieved. For instance, in cases where the center is already a good solution, e.g. Levy (Figure 13) and Rosenbrock (Figure 14), it can be difficult to improve upon initial conditions. In such cases, we can commit to a more exploitative search around XA. Based on this intuition, the focused method I+XA+F uses relatively short fixed lengthscales An = 0.1vD to compute distances do, and a fixed ro = 0.1. The remaining hyperparameters are still Icarned from data. As future work, several priors can be used in tandem. A portfolio can then be managed using strategies similar to those proposed by Hoffman et al. (2011) and McLeod al. (2018).
++XA+F+TR (Proposed) In order to further encourage the search around XA, I+XA+F can also be com- bined with trust regions, achieving better performance in certain cases, c.g. Levy (Figure 13) and Rosenbrock (Figure 14).
24

---

D.3 Results
QBranin,50D
0Bramin,100p
SH
S+
S-QM -AA
1.2
GS+XA
CM ES
1.4
1-X0 (Prpossed)
1-XA (Proposed)
1-A TT shspssed)
2.2
1-X0- (Popseet)
1-X6-1 (hoposed))
1-XA-F-7A (roposed)
103
120
200
200
120
100
120
200
$QBranin,200
SQBranin,50D
SQBranin,1000
1.8
S-0M
S+
Sm-OM -AA
G5+XA
CMAVES
I-XA (Propased)
(Poppsed)
2.00
SSQBranin.200
$5QBBrain.33D
$SQBranin,1000
100
S+ON
5-TR
00
04
[Poposed]
Horma
(Papssdi
[hrsposed)
(Proposod)
(Prpossed]
002
Figure 12: Performance QBranin and variants, ranging from 20 to 100 dimensions. Solid curves and
shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different inl conditions. Solid vertical lines indicate the interquartile range.
Lovy,20D
Levy,50p
Levy,100p
S+ON
S+TR
Proposed
Normalied
Propeced
1-X0A-OM (Toossed)
(hosssed)
0:2
++0A-FTR Propoeed
3.0
230
220
26
26
Sty Stybilinkk-Tang.2000
Styblliskk-Tang.5000
Noma 20 Inprovemsss 50
0.
SHQM--X
DS+KA
H+00 (Poposed)
Proposed
0.4
11XA: TK (Hoposed)
1:XA: QM (Propose)
i-3hot (7 (Propseed
1-XA-P-TR (Poppeed)
Figure 13: Performance on Levy and Styblinski ranging from 20 to 100 dimensions. Solid curves and
shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lincs indicate the interquartile range.
25

---

Rosenbrock,200
Rosenbrock.500
Rosenbrock,
Normaizzd Improvvmmett
S-MM+XA
CMNES
1+X0 (Propesec)
(Proposec)
TR (Proposedi
1+XA+ON Proposed)
1434 (frpposee)
[+XA+F+ (Proposoc)
$35Rosenbrack,,00
535R0senbrock,1000
Normatzee Imprr Improvement
SITTR
S-4MM-XA
GS-XA
CMA-ES
1+X0 (Proposed)
[Proposes]
I-XA+TR (Proposed]
1-XA+OM (Proposed)
(Proposed)
0.0
160
155
360
530Rosenbrock.500
550Resenbreck.100
Normalied www.ddipppooooooooooo Improvereen
S+OM
S
S+
GS
CMA-ES
0.4
1+X0 (Proposee)
1+XA [Proposec]
1+XA+TR (Prpossd)
0.2
(Prpposdd)
+TX (Proposed)
0.0
40
60
150
123
:40
14c
180
300
300
130
160
300
180
700
60
60
130
120
140
165
160
230
1.0
565Rosenbrock.208
565Rosenbrock.500
$65Rosembrock.1000
Narnalzzed Improvement
0.5
S++M
+XA
0.8-
G5+XA
0.4
(Proposed)
1+XA [Propssec
1+XA+TH (Poposed)
(Preposed)
1-XA+f (Prsposed)
1-XA+F+TR (Proposed)
60
150
123
100
100
120
100
130
160
100
250
Figure 14: Performance on Rosenbrock objectives with 20, 50 and 100 dimensions. Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range.
26

---

E
Ablation Study
In this section, we investigate the contribution of spatially-varying prior variances (V) and lengthscales (L) toward more efficient high-dimensional optimization. For this purpose, we test these spatially-varying properties individually, referring to them as I(V) and I(L).
Overall, the combination of spatially-warying prior variances and lengthscales is most effective in higher- dimensional spaces (50 and 100), where I+XO and I+XA generally outperform their I I(V) and I counter- parts. In comparison to S, IVV and I()) prove to be effective their own, but I I seems to account, for most of the performance increase on Rosenbrock (Figure 15), where anchors are relatively close to the opti- mum. The shorter lengthscales around anchors promote local exploration, even after having acquired data in their neighborhood. Interestingly, when the anchor misspecified (Styblinski- Tang), the performance of (LL) +XO better than that of I+XO and similar to that of in 20 dimensions. However, despite anchor mis- specification, we observe the usefulness of spatially-varying prior variances the dimensionality increases to 100, likely because these mitigate the boundary issue. For adaptive anchors, it is less clear whether the performance of I+XA on Styblinski- Tang is mostly due to spatially-varying prior variances or lengthscales. Both seem to be important because I+XA consistently outperforms I(V) +XA and +XA)
Rotenbrock:20
Resenbrack,500
Rasenbrack, 100D
Normatized Improsarmen
(Popoieds (1 (roppeedd
Haakk (rrppaedd
Rosenbrrok.200
Resenbrack,500 55000
Rosenbrock.100D
S+0W
Bomm aed Indernumeee
(b0+-6 (Poporeds
0.2
S2yblieatl-lbnng.200
Ryyblommt-/hng.11000
5104
Womm and Impeuume
(00130 (hopesati
H1 (Popocer
0.0
Stybiinkh-hingg
StyblinskiTann..oo
StykiinckiTnng. LO0D
IS thapsedd 1LKISA (hopondt
Figure 15: Ablation study. Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range. Abbreviations: Stationary (S), Cylindrical (C) and Informative (I) covariances; Quadratic Mean (+QM); Origin (+x0) and Adaptive, greedily-chosen (+XA) anchors; spatially-varying prior Variances (V) and Lengthscales (L).
27

---

F
Sensitivity Analysis
The informative covariance functions used in our experiments introduce a ratio ro E (0, 1], which leads to an uninformative stationary GP prior as it approaches 1. In terms of optimization, this hyperparameter balances global and local exploration under the informative search model (Equation 11) by adjusting the magnitude of spatially-varying prior (co)variances via To (Equation 13) and lengthscales through ra (Equation 15). For example, compared the stationary case, a ratio of ro 0.1 indicates that the prior variance is up to 10x larger and the squared lengthscales are 10x shorter in the neighborhood of anchor Zn.
By default, the ratio ro is learned by empirical Bayes and its prior is an informative Kumaraswamy density function with a peak at 0.1, given by Kumar(3.100 164, 1000) and referred to as K3. 4 The Kumaraswanny distri- bution is related to the beta distribution and is similarly parametrized by two positive shape parameters, a and b, with probability density function
PKumar (r: a, b) = abr r9)a-1,
T (0, 1).
(42)
Notably, both [ 3) and Kumar[:; b) are special cases of the generalized (p: 7, 6)-bela distribution (Jones, 2009, Equation 2.1), given by (1, 3) and (a, b), respectively. However, the latter is more tractable because the beta function in the denominator simplifies as B(1, b) = b- For a list of similaritics and other advantages over the beta distribution, see (Jones, 2009, Section 7).
In the remainder of this section, we examine the performance of I+XA under different priors for the hyper- parameter ro- We tested uniform (U) and Dirac delta (D) priors, as well as several Kumaraswamy priors, with parameters chosen to yield densities that become increasingly narrower, while maintaining a peak at 0.1. In particular, Kumar(3. 164, 1000) denoted by K3. Among these, K1 is the broadest prior, shown in Figure 16. The use of priors that favor smaller values of ro is motivated by the belief that the informative search model is useful. In this case, strategies that promote local exploration under this search model are more likely to improve upon initial conditions than those using an uninformative, uniform search model (stationary GP prior).
tested Kumar(11467, 10) denoted by K1, Kumar(2.253, 100) denoted by K2 and
Figure 17 shows the performance of I+XA variants on QBranin, Rosenbrock and Styblinski-Tang. We begin by observing that, on QBranin, I+XA is relatively robust to the choice of prior over ro because the performance of all variants is similar and significantly better than that of S, S+QM and C. However, the diffcrences become more not iceable on Rosenbrock and Jtyblinski-Tang The uniform prior leads to the worst-performing variant, but I+XA+U can still outperform the other baselines on Rosenbrock. Conversely, narrower priors, such as K2 and K3, which promote exploration under the informative search model, lead to the best-performing methods. Interestingly, I+XA with a delta prior (I+XA+D) can outperform other variants on Rosenbrock. This result suggests that learning all parameters via marginal likelihood may not be the best approach. Alternative training objectives that also consider optimization performance may provide better results, but this study is left for future work.
2-
8.0
022
0.4
0.6
Figure 16: Probability density functions of the priors over To-
"The suffix +K3 omitted in other sections because it the default prior over ro- As discussed Appendix D, the only exception is I+XA+F that uses ro = 0.1, corresponding to the Dirac delta (D) prior that is tested in this section.
28

---

QBranin,20D
0Branin,50D
QBranin, 100D 1000
Normatzzd impprovamin hopoovemeen
I-XA+U (Propsed)
(Proposed)
X-XA+K2 (Proposedi
I-XA+K3 (Proposed)
I-XA+D (Propesec)
Rosenbrock,20D
Rosenbrock,500
Rosenbrock,1000
S+0M
Normatzee Imprr improvement
[+XA+U (Propesed)
(Proposed]
I+XA+K2 (Proposed
1+XA+K3 (Proposedi
(Propcsed)
180
100
135
360
20D
Styblinssk/Tang.500
Styblinski-Tnng.oor
S-QM
Normalieed Improoorrennn peproveeent
(Proposed)
1+XA-K1
(Proposed)
1+XA+K2 (Proposed)
1+XA+K3 (Proposed)
0.4
0.2
0.8
32
40
20
150
120
240
140
180
200
360
120
360
300
180
700
43
60 1.0
120
140
160
160
230
Figure 17: Performance of I+XA under different priors ro- Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid verlical lines indicate the interquartile range. Abbreviations: Stationary (S), Cylindrical (C) and Informative (I) covariances; Quadratic Mean (+QM); Adaptive, greedily-chosen (+XA) anchors; Uniform (+U), Kumaraswany (+K) and Dirac delta (+D) priors.
29

---

G Rover Trajectory Planning
In this section, we examine the rover trajectory problem proposed by Wang et al. (2018), which has also been adopted by more recent studies (Eriksson et al., 2019; Eriksson & Jankowiak, 2021). In brief, the goal is to optimize a 2D trajectory such that the rover starts as close as possible to Estart and stops near Bend: while avoiding collisions with objects. The trajectory is given by a B-spline that is fitted to 30 2-dimensional points, resulting in a 60D optimization problem. The original reward function is defined
((x) c(x:) + (I[(=)1; &start || + || Fend (11) b,
x E [0, 1
(43)
where y = 10, = and d(x) is a nonpositive function that penalizes collisions. For minimization, we turn the reward into loss function. In particular, we take the negative reward and apply a +b shift, yiclding a non-negative objective to be minimized. Models are again estimated log-transformed data. The domain is also adjusted to the centered hypercube X [-1, 1] by applying linear transformation to the inputs.
Figure 18 shows the original map layout with two overlaid trajectories found by S and I+XA (left), well as the corresponding performance curves (right). In this problem, we immediately observe that the trajectories do not pass through the 2D points. This outcome is due to oversmoothed splines, which explains the performance of S. Intuitively, given the map layout, we would expect points near the boundary to yield penalized trajectories, but this is only true if the trajectories pass through the specified points. Instead, points the boundary have a significant. effect the oversmoothed trajectories. The overexploration of boundarics that is characteristic of S in higher-dimension However, observe that I+XA is eventually able to find trajectories of similar cost.
problems confers an advantage in this situation.
In order to align the problem formulation with our expectations, we modify the procedure that fits B-splines to the collection of 2D points. In terms implementation, the function call that performs interpolation, scipy. interpolate splprep (Virtanen el al., 2020), is extended with a smoothing factor s that set 0, which forces the trajectories to pass through all points. An example is shown in Figure 19 (left), Additionally, we clamp all trajectories into the admissible range and increase the cost for deviating from the endpoints, 7 50. The latter promotes unfolded trajectories away from the initial incumbent solution (ma) 0, as the relative cost of collisions too steep otherwise. Performance curves for all rover tests are x best
shown in Figure 20. Notably, we observe that S is longer among the best-performing methods because the specification of points near the boundary correctly yields penalized trajectories.
16 898)
1+XXA Icost= (c05t=3.580)
OversmoothedRoverT..60
1.0
Normalizzd improvement
S+QM
0.8
S+TR
GS+XA
0.6
CMA-ES
1+X0 (Proposed)
I+XA (Proposed)
++XA+TR (Proposed)
++XA+QM (Proposed)
0.2
I+XA+F (Proposed)
++XA+F+TR (Proposed)
0.0
0
40
60
80
100
120
140
160
180
200
Figure 18: Original map layout with overlaid rover trajectories (left). Trajectories are given by B-spline that is fitted to the 30 2-dimensional points. Performance on the original rover trajectory problem (right).. Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range.
30

---

ShitteGGoalRovart
Fajectory
Figure 19: Original map layout with an overlaid rover trajectory (left). In RoverT, trajectories must pass through the specified collection of 2D points. The handcrafted trajectory demonstrates that not necessary to specify 30 different points. However, the reward function does not penalize less efficient trajectories, which leads an infinite number of equally-good trajectories. As an additional test, we consider the same map layout with shifted endpoint Lend (right).
OversmoathedRoverT.cc
RoverT,60D
ShilteddGollRovrr..600
5-
No malllazd Impro Improverrats
0.8
G5
0.5
CMA ES
(Proposed)
(Propesec)
0.4
1+XA+ ON (Propesed)
14XA-T8 (Proposedi
(Proposed)
(Proposed]
0.2
150
322 150
200
200
120 200
40
120
140
200
Figure 20: Performance on the rover trajectory problem. Solid curves and shaded regions represent the and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range.
As shown Figure 19 (left), an optimal trajectory does require the specification of 30 different points. However, the reward function does not take into account the distance traveled, which means that the is free to roam the environment, provided that it starts at 2strt and ends at Fend without colliding with other objects. For this reason, trajectories of similar cost can be significantly different and the choice of GP priors becomes even more important.
The importance of GP priors further demonstrated in Figure 21, which shows an example of the best RoverT solutions found by BO. In particular, the stationary covariance in (lop-left) does not promote exploration around any particular region. In this case, the best trajectories tend to cover the entire space, while trying to avoid collisions. In constrast, the quadratic mean prior (top-middle) and the cylindrical covariance function (top-right) promote exploration near the center. As quadratic weights can be nearly zero, some points along the trajectory can also be near the boundary. In both cases, the prior dominates to such an extent that the endpoints of the best trajectories are comparatively far from the target endpoints. Notably, I+XO (bottom-left) also encourages exploration around the center, but to lesser degree than S+QM and C, leading better solutions. Finally, I+XA (bottom-middle) and I+XA+F (bottom-right) find the most cost-effective and efficient trajectories because both promote exploration around the incumbent. solution, emphasizing the search for fine-tuned trajectories.
31

---

SccqqeD, costa50 600)
S+ (acquo, costu50 000)
cacque,, cost=90.000)
-accsm50.com costmtd.dcoooom 6621
5+-MM (scom500 costm90.com
(engm10, cost-Bt 9121
(acqm100, cosl-50.592
p010
S4-M (acqm100 100, costmcs.e90)
(acqm1on, 100, cosimta A85)
(aco= 150,
(acqe150,
(Acqr200,
(acgr200,
H+D0 lacouo, 90st=90 2001
1+X0 (acqu-) cost-500
1+XA++ (scquO, cost=90
#+00 (Sco=50, c35t=50.2999 298)
1+30 (acq-50
cost=
1eXA+F (aco-50. cost= 2271
H00 (aco-100, 5454-51.344)
1+30 (acq-2200 2ost- 5141
I+XA++ (acq-100. 2351-23.844
(aco-150, cost-31.431)
(SC3-150, cost-
[-XA-FF (Scq- 150--150.000 c051-27.7411
(acqu-200
Figure 21: Example of the best rover trajectorics found by BO using different GP priors. Top row: baselines s, S+QM and C. Bottom row: proposed methods I+XO, I+XA and I+XA+F.
32

---

H Sparse Axis-Aligned Subspaces (SAAS) Hyperprior
The Sparse Axis-Aligned Subspaces (SAAS) hyperprior, introduced by Eriksson & Jankowiak (2021), assumes that only a subset of dimensions affects the objective function. This sparsity-inducing hyperprior posits that the inverse squared lengthscales are distributed according to a half-Cauchy, 1/1 ~ HC(T), where 7 is a global shrinkage hyperparameter.
each step, multiple GP models are trained by empirical Bayes with 7 {1, 10-1 10-2, 3310-3. The GP model with the highest leave-one-out coss-validation likelihood then used for acquisition.
We evaluated modified version of S, which incorporates SAAS and follows the original training procedure described above. As shown in Figures 22 and 23, S+SAAS generally performs worse than 8, indicating that this hyperprior may adversely affect performance when the low effective dimensionality assumption does not hold. To further demonstrate the complementary nature of our proposed methodology. we tested I+XA+SAAS. This combination proved to be beneficial on QBranin, Rosenbrock and Levy. Recall that, in I+XA, the lengthscales also govern the nonstationary effects, and incorporating SAAS allows for longer lengthscales (c.f. Table 2). On the Styblinski- Tang function, however, this hyperprior continued to negatively impact performance, renewing the question of its suitability when the objective does not cxhibit low effective dimensionality.
Q0D
QBranIn,50D
QBranin,1000
2.0
S40M
Normaibed improvemmnt
0.8
S+54A5
H+X0 (Proposed)
KA (Proposed)
H K (Proposed)
0.4
(Proosed)
0.24
140
1.20
Rosenbrock, 20D
Rosembroct 500
Rosanbrock,
0.8.
HX00 (Proposed)
1AXA (Proposed)
SAAS (Proposed)
(Pooposel)
0.0g
30
60
80
100
320
150
350
200
50
330
100
350
180
200
20
40
87
80
133
340
180
3.00
300
Levy,20D
Levy,30D
Levy,100D
S+SAAS
Nomaized Imirovemerr
S
0.6
14X0 (Proposed)
(Proposed)
(Proposed)
43
300 140
230
1.0 248 185
100
340 200
Figure 22: Performance on test functions, ranging from 20 to 100 dimensions. Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over trials with different initial conditions. Solid vertical lines indicate the interquartile range, Abbreviations: Stationary (S), Cylindrical (C) and Informative (I) covariances; Quadratic Mean (+QM); Origin (+X0) and Adaptive (+XA) anchors; acquisitions within Trust Region (+TR); Sparse Axis-Aligned Subspace (+SAAS) hyperprior.
33

---

Styellinsk-Tangg2200
Styblliskk-Tann... 50D
Stybllinski-Tnng10000
S+QM
Normaizzo Improvement
0.8
5-5AA5
1+30 (Proposed)
I+XA (Proposed)
1+XA+ 5MAS (Proposed)
[4XA+TA (Proposed
0.2
300
130
502
Figure 23: Performance the Styblinkki-Tang function with 20, 50 and 100 dimensions. Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range. Abbreviations: Stationary (S), Cylindrical (C) and Informative (I) covariances; Quadratic Mean (+QM); Origin (+X0) and
Adaptive (+XA) anchors; acquisitions within Trust Region (+TR); Sparse Axis- ArisnAl Subspace (+SAAS) hyperprior.
34

---

Augmented Acquisition Functions
In order to incorporate beliefs about optimal solutions into BO, Hvarfner et al. (2022) proposed prior- weighted acquisition functions, which are obtained by multiplying an acquisition function a with a prior distribution over possible locations of the optimum. As oplimization progresses, the influence of the fixed prior T(x) decays according to ag,n & a( D): where the hyperparameter Ç controls the docay. As more acquisitions
made, the prior approaches uniform distribution, only tilting the acquisition step initially, when the exponent Ç/n is not small. By default, Hvarfner el al. (2022) used Maussian-Weighter ET (GWEI) and Ç ç N/10, where N is the evaluation budget. The default Ganssian prior is characterized by a diagonal covariance matrix, with standard deviation set to 25% of the domain.
We evaluated two additional baselines, S+GWEIXO and S+GWEIXA, which use Gaussian- Weighted EI with fixed and adaptive locations. However, the results on high-dimensional test functions indicate that these basclines perform poorly when compared to S, shown in Figures 24 and 25. We found that poor performance was also linked to machine precision because EI, whose values are already small, is multiplied by small factors, resulting in an acquisition function that is difficult to optimize in high dimensions. Tn order to improve performance, removed the normalizing constant in the Gaussian prior, obtaining Gaussian Kernel- weighted El (GKEI). The method I+XA+GKEI shows that combining informative covariance functions and belief-augmented acquisition functions can an effective strategy.
2r0
QBranin.SOD
OBranin,100D
1.0
S+0M
Improvaneets
S1GWEK0
SHGWEKA
0.0
(Propese5)
14 (roposee)
0.4
4+ GXTTA (Propssed)
14 X G (Proposed)
14 6- Eil (Poposel)
120
200
1.30
Rosenbrock.200
Rosembres 30D
Rosenbrock.100D
3.0
S+OM
S1GWBBK
S1GWEKA
(Proposed)
X+XA (Propssdd
Normaliee
0.4
XXTCXXIX (Proposed)
GXED (Peppsedd
(Propsedi)
0.2
0.0
L0V
Levy,50D
Levy,100D
S
Normalized Imprevemions
0.8
14X0 (Proposed)
1+XA (Proposed)
GXEIKO (Proposed)
GXEIA (Poposed)
(Poposed)
0.0
Figure 24: Performance test functions, ranging from 20 100 dimensions. Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range. Abbreviations: Stationary (S), Cylindrical (C) and Informative (I) covariances; Quadratic Mean (+QM); Origin (+x0) and Adaptive (+XA) anchors; acquisitions within Trust Region (+TR); Caassian-Weighhee Expected Improvement with fixed (+GWEIX0) and adaptive (+GWEIXA) location; Gaussian Kernel-weighted Expected Improvement with fixed (+GKEIXO) and adaptive (+GKEIXA) location.
35

---

Styellinskd-Tan.22000
Styblliskk-Tang..555 50D
Stybliinski-Tnng.1ooo
S+QM
Noraaizz Improvement
0.8
S-GNNIXO
1+X0 (Proposod)
[+XA (Proposed)
G (Proposed)
GXEIKA (Popozed)
02
60
300
130
202
Figure 25: Performance the Styblinski-Tang function with 20, 50 and 100 dimensions. Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range. Abbreviations: Stationary (S), Cylindrical (C) and Informative (I) covariances; Quadratic Mean (+QM); Origin (+X0) and Adaptive (+XA) anchors; acquisitions within Trust Region (+TR); Gaussian- Weigghtd Expected Improvement with fixed (+GWEIXO) and adaptive (+GWEIXA) location; Gaussian Kernel-weighted Expected Improvement
with fixed (+GKEIXO) and adaptive (+GKEIXA) location.
36

---

J High-dimensional Experiments with Larger Evaluation Budgets
Thus far, we have considered a budget of 200 acquisitions. As highlighted in Sections 2.1 and 3, the rationale is that BO is competitive when function evaluations are expensive, putting a practical constraint on the size of the budget, even for high-dimensional problems. Additionally, we have favored an analysis in terms of sample efficiency, demonstrating that the proposed methodology, which employs informative covariance functions, offers methods that discover superior solutions with fewer function evaluations.
Despite the rationale above, we increase the evaluation budget to 600 acquisitions. Figure 26 shows the results for two test functions that we consider representative (see Appendix B.3). The Rosenbrock function features banana-shaped valley, with the global optimum relatively close to the center of the search space. Conversely, the Styblinski- Tang is multimodal function, whose optimal solution is relatively distant from the center (local maximum). Once again, we observe that the proposed covariance functions can be combined with other methods, including quadratic mean functions and trust regions, to achieve higher sample efficiency and discover superior solutions.
As final note, while methodology provides scalable methods in terms of dimensionality (c.g., I+XO, I+XA), these still share some of the limitations of standard GPs. In particular, when dealing with large training sets, standard GPs become impractical because they require O(n") computation and O memory, where n is the number of training points. To address this limitation, one option is to use sparse GPs with inducing variables, as proposed by e.g. Titsias (2009) and Hensman et al. (2013). This allows for a significant reduction in computational and memory requirements due to low-rank matrix approximations.
boseebrock,100D
StyblinskiTang.1000
1.0
S+QM
Normal bed improvement
0.8
S+TR
[+X0 (Proposed)
(Proposed)
0.6
QM (Proposed)
1+XA+TR (Proposed)
0.2
0.0
0
50
100
150
200
250
300
350
400
450
500
550
600
50
100
150
200
250
300
350
400
450
500
550 600
Figure 26: Performance on the 100-dimensional Rosenbrock and Styblinski-Tann Tang functions. A value of indicates that the global optimum has been found (zero simple regret). Solid curves and shaded regions represent the mean and standard deviation of the normalized improvement, computed over 10 trials with different initial conditions. Solid vertical lines indicate the interquartile range. Abbreviations: Stationary (S), Cylindrical (C) and Informative (I) covariances; Quadratic Mean (+QM); Origin (+X0) and Adaptive (+XA) anchors; acquisitions within Trust Region (+TR).
3ased similar reasoning, Friksson & Jankowiak (2021) conducted 100-dimensional experiments with maximum budget. of 100 acquisitions.
37


Reviews:
[{"review_id": "0DR82fIgi6", "reviewer": "Reviewer_7J43", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, I think the authors' idea is great. This seems like a pretty nice way to get TR-like optimization behavior without necessarily having to use a TR (although the methods are complementary, as the authors point out), and in several experiments this seems to pan out in terms of getting fairly good optimization performance at lower evaluation budgets than we might expect using a TR. I especially like that a lot of the resulting parameters end up being learnable, which perhaps makes the authors' method require less hand crafting than some of the TR hyperparameters, which admittedly have reasonable out of the box settings.\n\nOn the weaknesses side, I didn't find the discussion around regret to be particularly useful. Some of the earliest approaches to provide regret bounds for Bayesian optimization in the setting where the hyperparameters aren't known in advance already view the objective function in that setting as drawn using a set of RKHS's indexed by a sequence of lengthscales, (e.g. Bull et al), and I think it's reasonably well known that the resulting regret bounds have to be pessimistic by analyzing under the shortest possible lengthscales. I understand that this section primarily serves as motivation for the authors' practical method, but the first bullet point at the end of Section 1 is perhaps a bit overstated as a result.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this paper, the authors introduce nonstationary GP priors useful for Bayesian optimization. In particular, the authors assume the existence of a set of anchors (most commonly L=1 anchor in the paper, which could be taken e.g. to be the current incumbent solution). The prior variance is then taken to be a function that smoothly inflates the correlation of two points when they are near an anchor, and the lengthscales are defined by an input warping (e.g. Snoek et al., 2014) that shrinks lengthscales in a neighborhood around the anchors (e.g., the incumbent). These factors lead to the influence of the anchors being greatly expanded, resulting in the optimization procedure being naturally encouraged to explore more closely around anchors. ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- The plots can become fairly illegible with as many as 5 of the authors own methods appearing in many of the result plots. The authors really should consider picking 1-2 instantiations of their method as exemplars, and then leave e.g. the comparison of XA to X0 to a separate ablation study. In half of these plots, I can't tell if the error bars e.g. between I+XA and S+QM overlap because there are 4 other instances of their own method on the plot.\n\n- The evaluation budget considered in every task is quite small. This is fine, because BayesOpt often targets the very low evaluation budget regime. However, I still think there's substantial value in at least *showing* us results as an ablation study with larger budgets on a few of these tasks. For example, S+TR does extremely poorly by N=200 observations in the authors' results on Rover, but obviously does fairly well on this task by N=7500 (Eriksson et al., 2018). Currently, the long term trade-offs of the authors' method is a fairly open question, and it's not obvious to me that the authors' method should have any scalability concerns, so the lack of these results is confusing to me.\n\n- [Minor] I'd like the authors to add a piece of clarification about the prior variance function. In particular, it looks to me as though the prior covariance function in (12) will inflate the covariance between two points that are close to anchors, *regardless* of whether they are close to the *same* anchor. Trivially this is a nonissue in the case of one anchor (equation (13)), which the authors primarily focus on, but I'm curious whether this implies that, if anchors aren't chosen to have very similar function values (e.g., if the second best observation in the domain far from the first best happens to be much worse), this might lead to somewhat strange model fits.\n", "broader_impact_concerns": "No real concerns here.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "2pReTkpuwu", "reviewer": "Reviewer_DmEg", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n\nThe writing is clear and articulate. The formulation of spatially varying prior variance and lengthscales given one or more anchor points is novel and intriguing. Although one might not always have prior knowledge of promising anchor points, a technique for adapting the anchor point as the current best incumbent is proposed and empirically validated. The regret analysis and discussion motivated the spatially varying kernel. Moreover the work highlights how the proposed kernel is complementary to prior-guided BO and trust region BO.\n\n\nWeaknesses\n\nThe main weakness is the empirical evaluation (see requested changes for comments). Compelling real world examples would greatly strengthen the paper. Moreover, since the empirical evaluation exclusively considers high-dimensional problems, additional high-dimensional BO baselines are warranted. It would strengthen the paper to include lower dimensional test functions---at a minimum to provide evidence for why non-stationary kernels are more beneficial/important for high-dimensional problems. If this paper is solely focused on high-dimensional BO, I think that decision should be better motivated (at least empirically) and that focus on high-dimensional BO should made clear early in the paper (and potentially in the title).\n\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper considers Bayesian optimization (B) with non-stationary kernels. Typically, stationary kernels are used in BO, but this work argues both empirically and theoretically that informative, non-stationary GPs can lead to improved sample efficiency. Specifically, this work introduces novel spatially-varying lengthscales and prior variance functions that can exploit prior knowledge about promising regions, provides a regret analysis, and demonstrates improved optimization performance on high-dimensional test problems. ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "* [Important] Why are there only high-dimensional test functions? An analysis of the performance of this kernel in lower dimensional problems would strengthen the paper.\n    * Since the evaluation is only on high-dim problems, how do other SOTA high-dimensional BO methods like SAASBO (Eriksson et al., 2021) compare?\n* [Important] More compelling real world examples would drive home why this is important\n    * Can the authors include a real world example a reasonable anchor point is known apriori?\n    * The NN problem is not very compelling real world use case\n* [Nice to have] How does performance vary with different nonstationary choices (e.g. no tying r,d,k)?\n* [Nice to have] How does performance vary based on the number of anchor points?\n* [Nice to have] What about vanilla input warping (Snoek et al., 2014)?\n\nAdditionally, the organization of the discussion of the results hinders readability.\n* In the fixed anchor section, it would be good to discuss poor performance on Styblinski Tang, before getting to the adaptive anchor point technique. I can see that the paper aims to craft a story here, but as a reader, I was wondering why the poor performance of the fixed anchor method was not discussed in the fixed anchor section. Perhaps starting with an introduction to all methods and then discussion the results in the figures would make the work easier to follow.\n\n\nNit: Page 4: \\citet at the end of Section 2.5", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Mx66p1VSHV", "reviewer": "Reviewer_vdWc", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n+ The paper investigates an important and interesting problem, that is, using spatial-varying nonstationary covariance functions in Bayesian Optimization, which might help to improve the performance of BO in high-dimensional or complex optimization problems.\n+ I like the idea of constructing spatial-varying nonstationary covariance functions by encoding the preferences for certain regions.\n+ The descriptions of the experiments are detailed and thorough. Highly dimensional test functions are used to compare the performance between the proposed method and other baseline methods.\n\nWeaknesses:\n\nEven though I like the general idea of the paper, I found there are various issues with the current version of the paper:\n+ One of the contributions listed by the paper is to provide a regret analysis that demonstrates the advantages of non-stationary covariance functions with spatially varying properties. However, Section 3 only provides some useful insights regarding the regrets of the spatial-varying nonstationary covariance functions, and there are no rigorous theorems proving that the spatial-varying nonstationary covariance function will tighten the regret bound. A lot of discussion/analysis in Section 3 should be formally formulated as theorems to guarantee their correctness. Lemma A.1 is not really helpful (1) it relies on an assumption that is not guaranteed to occur, (2) its conclusion does not help to understand the property of the regret bound as\n+ The paper aims to provide a general formula for spatial-varying nonstationary covariance functions with a set of anchors, but it eventually only focuses on covariance functions with only one single anchor – which is very limited. And there are various design choices regarding the set of anchors (or the single anchor) that are not mentioned. For example, how is the distance function d_l chosen?  What is the range of w_l? Does it need to be larger than 1? How can we ensure \\phi(x) to be positive? Why, in Eq. (11), x* is used, but not x? I feel really confused when reading the equation. What is the effect of choosing the wrong positions for the anchors? Currently, the experiment’s setup is that the global optimum is at the centre of the search space, and the anchor is chosen based on that. However, when we don’t know the location of the global optimum, how can we choose the anchors to ensure good performance?\n+ The adaptive anchor strategy is not clearly described in the paper. There is only one sentence mentioning it in the experiment section (adaptive anchor section). It should be described clearly in Section 4 if it is part of the proposed method.\n+ The experimental evaluation is quite limited, in my opinion. First, it is only conducted with one single anchor, which does not help to understand the generic method. Second, the number of iterations in each test case is too small (only 200). Normally for highly dimensional optimization problems, the number of iterations needs to be chosen much larger to understand the performance of the proposed method. Third, the normalized regret is used, which is a non-standard metric. Besides, in Section 3, the paper provides analysis regarding the standard regret, but in the evaluation, the normalized regret is used. Finally, the proposed methods used in the synthetic experiments are not the same as the proposed methods used in the real-world applications (in real-world experiments, there are I+XA+QM and I+XGood). The use of the proposed methods should be consistent in both synthetic and real-world experiments.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper investigates the strengths of using spatial-varying nonstationary covariance functions in Bayesian Optimization. The paper proposes a new class of covariance functions that leverage nonstationary to incorporate information about promising points in the search domain. The paper evaluates their proposed method on some high-dimensional test functions and real-world applications, under some prior information. It also shows that it can combine with other methods to improve its performance.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "+ More mathematical rigorous theorems to be included when analyzing the advantages of non-stationary covariance functions with spatially varying properties if this is listed as a contribution of the paper. Please see my comments in the previous section.\n+ The proposed method needs to be elaborated more. For example, more discussion on the use of a set of multiple anchors, and more discussion regarding the design choices of the proposed method. The adaptive anchor strategy needs to be described in detail. Please see my comments in the previous section.\n+ Experiments need to be improved. In particular, the number of iterations needs to be increased. The performance of all the methods with the standard regret also needs to be included in the experimental evaluation. The proposed methods to be used in synthetic and real-world experiments need to be consistent. Please see my comments in the previous section.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Based on the reasons above, it's not an easy decision. I'll just recap the main points: (+) useful idea, that works in practice, but (-) some (not particularly challenging) experiments are missing. All reviewers are hesitating, 1 recommends to reject [vdWc: "My opinion towards the paper is quite marginal (...) which is not enough to convince me towards the acceptance of the paper"] while 2 think that the pros outweigh the cons [7J43: "I'm leaning towards acceptance here (...) I think the actual method here is actually pretty sensible, and the kind of thing you might actually try adding to your bayesopt code in practice with pretty minimal overhead. To me, that makes it a lot harder to outright reject the paper. Honestly, if they'd just come back and run every result with budgets up to a couple thousand, I'd be championing this paper strongly at this point"; DmEg: "I lean slight toward acceptance, but I acknowledge the criticisms raised by Reviewer vdWc and I would not oppose a rejection", "Despite that, I do think this work would be of interest to the community and that its merits outweigh its weaknesses"]

Because I agree that this method should be in the toolbox of practitioners of Bayesian optimization, I will follow the opinion of Reviewers 7J43 and DmEg and recommend to accept the paper. We do not completely understand why the authors did not add more large budget experiments, this would have made the experimental section much stronger at little cost. I will recommend "accepted under minor revision" to allow the authors to include more experiments. Of course, at that stage, there won't be any more reviews, but I insist that the authors decision to include more experiments could make their own work more convincing. It's their choice ultimately -- on our side, we believe it is a useful contribution, and as such it deserves to be published.
Decision: Accept (poster)