The Behavior and Convergence of Local Bayesian Optimization
Authors: Kaiwen Wu
Kyurae Kim
Roman Garnett
Jacob R. Gardner
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Bayesian optimization
convergence rates
TL;DR: We provide the first convergence rate for a local Bayesian optimization algorithm in both the noisy and noiseless settings.

Abstract: A recent development in Bayesian optimization is the use of local optimization strategies, which can deliver strong empirical performance on high-dimensional problems compared to traditional global strategies. The "folk wisdom" in the literature is that the focus on local optimization sidesteps the curse of dimensionality; however, little is known concretely about the expected behavior or convergence of Bayesian local optimization routines. We first study the behavior of the local approach, and find that the statistics of individual local solutions of Gaussian process sample paths are surprisingly good compared to what we would expect to recover from global methods. We then present the first rigorous analysis of such a Bayesian local optimization algorithm recently proposed by Müller et al. (2021), and derive convergence rates in both the noisy and noiseless settings.

The Behavior and Convergence of Local Bayesian Optimization
Kaiwen Wu University of Pennsylvania kaiwenvudseas upenn. edu
Kyurae Kim
University of Pennsylvania kyrkimsseas upenn.
Roman Garnett
Jacob R. Gardner University of Pennsylvania jacobrgeseas upenn. edu
Washington University in St. Louis garnett@wustl.edu
Abstract
A recent development in Bayesian optimization is the use of local optimization strategies, which can deliver strong empirical performance on nigh-dimensional problems compared traditional global strategies. The "folk wisdom" in the liter- ature is that the focus on local optimization sidesteps the curse of dimensionality; however, little is known concretely about the expected behavior convergence Bayesian local optimization routines. We first study the behavior of the local approach, and find that the statistics of individual local solutions of Gaussian process sample paths are surprisingly good compared to what we would expect to recover from global methods. We then present the first rigorous analysis of such a Bayesian local optimization algorithm recently proposed by Müller et al. (2021), and derive convergence rates in both the noiseless and noisy settings.
Introduction
Bayesian optimization (BO) [11] is among the most well studied and successful algorithms for solving black-box optimization problems, with wide ranging applications from hyperparameter tuning [23], reinforcement learning [19, 27], to recent applications in small molecule [18, 26, 7, 13, 12] and
sequence design [25]. Despite success, the curse of dimensionality has been long-standing challenge for the framework. Cumulative regret for BO scales exponentially with the dimensionality of the search space, unless strong assumptions like additive structure [15] are met. Empirically, BO has historically struggled on optimization problems with more than handful of dimensions.
Recently, local approaches to Bayesian optimization have empirically shown great promise in addressing high-die sigh-dimensional optimization problems [10, 19, 20, 27, 17J. The approach is intuitively appealing: typical claim is that local optimum can be found relatively quickly, and exponential sample complexity is only necessary if one wishes to enumerate all local optima. However, these intuitions have not been fully formalized. While global Bayesian optimization is well studied under controlled assumptions [24, 3], little is known about the behavior or convergence properties of local Bayesian optimization these same settings the related literature focuses mostly on strong empirical results on real-world applications.
In this paper, begin close this gap by investigating the behavior of and proving convergence rates for local BO under the same well-studied assumptions commonly used to analyze global BO. We divide our study of the properties of local BO in two questions:
Q1 Using common assumptions (e.g., g., when the objective function a sample path from GP with known hyperparameters) how good are the local solutions that local BO converges to?
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

GIBO (0
GIBO (o 0.00)
GIBO =0
-5
Min
-10
-10
-10
-15
-15
10
20
30
50
1
5
10
20
30
50
5
10
20
30
50
Local Opt. Elfective Size
GP-UCB - 0))
Random Search
-5
-5
Min value
-10
-10
5
20
30
5
20
30
20 dimension df
dimension d
dimension d
The unreasonable effectiveness of locally optimizing GP sample paths. (Top row): Distributions of local solutions found when locally optimizing GP sample paths in various numbers of dimensions, with varying amounts of noise. (Bottom left): The minimum sample complexity of grid search required to achieve the median value found by GIBO (a 0) in expectation. (Bottom middle, right): The performance of global optimization algorithms GP-UCB and random search in this setting. See $3 for details.
Figure 1:
Q2 Under these same assumptions, there local Bayesian optimization procedure that guaran- teed to find a local solution, and how quickly can it do so?
For Q1, find empirically that the behavior of local BO on "typical" high dimensional Gaussian process sample paths is surprising: even single run of local BO, without random restarts, results in shockingly good objective values in high dimension. Here, we may characterize "surprise" by the minimum equivalent grid search size necessary to achieve this value in expectation. See Figure and more discussions $3.
For Q2, we present the first theoretical convergence analysis of a local BO algorithm, GIBO as presented in [19] with slight modifications. Our results bring focus number of natural questions about local Bayesian optimization, such as the impact of noise and the increasing challenges of optimization local BO approaches stationary point.
We summarize our contributions as follows:
1. We provide upper bounds on the uncertainty of the gradient when actively learned using the GIBO procedure described in [19] in both the noiseless and noisy settings.
2. By relating these uncertainty bounds to the bias in the gradient estimate used at each step, we translate these results into convergence rates for both noiseless and noisy functions.
3. We investigate the looseness in our bound by studying the empirical reduction in uncertainty achieved by using the GIBO policy and comparing it to our upper bound.
4. We empirically study the quality of individual local solutions found by a local Bayesian optimization, slightly modified version of GIBO, GP sample paths.
2 Background and Related Work
Bayesian optimization.
Bayesian optimization (BO) considers minimizing a black-box function f:
minimize fxx xEX
through potentially noisy function queries = -xx) + e. When relevant to our discussion, we
will assume an i.i.d. additive Gaussian noise model: (N0 o²). To decide which inputs to

---

be evaluated, Bayesian optimization estimates f via surrogate model, which in turn informs policy (usually realized by maximizing an acquisition function (xx over the input space) identify evaluation candidates. Selected candidates are then evaluated, and the surrogate model is updated with this new data, allowing the process to repeat with the updated surrogate model.
Gaussian processes. Gaussian processes (GP) are the most commonly used surrogate model class in BO. A Gaussian process G (j, k) is fully specified by a mean function μa and a covariance function k. For a finite collection of inputs X. the GP induces a joint Gaussian belief about the function values, f(X) ~ ((((X, A:X,, ((XXX)) X)). Standard rules about conditioning Gaussians can then be used to condition the GP on dataset D, resulting in updated posterior process reflecting the information contained in these observations:
fo ~ GP(DD, KD).
Existing bounds for global BO. Srinivas et al. [24] proved the first sublinear cumulative regret bounds of a global BO algorithm, GP-UCB, in the noisy setting. Their bounds have an exponential dependency the dimension, result that is generally unimproved without assumptions function structure like an additive decomposition [15]. This exponential dependence can be regarded consequence of the curse of dimensionality. In one dimension, Scarlett [21] has characterized the optimal regret bound (up to a logarithmic factor), and CP-UCB is mear-optimal in the case of the RBF kernel. In the noiseless setting, improved convergence rates have been developed under additional assumptions. For example, De Freitas et al. [6] proved exponential convergence by assuming, f is locally quadratic in neighborhood around the global optimum, and Kawaguchi et al. [16] proved an exponential simple regret bound under an additional assumption that the optimality gap f(xx - S is bounded by a semi-metric.
the
Gaussian process derivatives. A key fact about Gaussian process models that we will use is that they naturally give rise to gradient estimation. f is Gaussian process distributed, I ~ ((, k), and k is differentiable, then the gradient process f also is also (jointly) distributed as a GP. Noisy observations y at arbitrary locations X and the gradient measured at an arbitrary location x are jointly distributed as:
Vo(x) -~ (((42))) ( (
X) + o²1 VA(X, X)
x(X,x)VT) V
This property allows probabilistic modeling of the gradient given noisy function observations. The gradient process conditioned on observations 'D is distributed as Gaussian process:
V D(VpD, V
3 How Good Are Local Solutions?
Several local Bayesian optimization algorithms have been proposed recently alternative to global Bayesian optimization due to their favourable empirical sample efficiency [10, 19, 20]. The common idea is to utilize the surrogate model to iteratively search for local improvements around some current input X, at iteration t. For example, Eriksson et al. [10] centers a hyper-rectangular trust region on x, (typically taken to be the best input evaluated far), and searches locally within the trust region. Other methods like those of Müller et al. [19] or Nguyen et al. [20] improve the solution x locally by attempting to identify descent directions d, that x, can be updated as X,+1 X, m.d,-
Before choosing a particular local BO method to study the convergence of formally, we begin by motivating the high level philosophy of local BO by investigating the quality of individual local solutions in a controlled setting. Specifically, we study the local solutions of functions f drawn from Gaussian processes with known hyperparameters. In this setting, Gaussian process sample paths can be drawn as differentiable functions adapting the techniques described in Wilson et al. [29]. Note that this is substantially different from the experiment originally conducted in Müller al. [19], where they condition a Gaussian process on 1000 examples and use the posterior mean as an objective.
To get a sense of scale for the optimum, we first analyze roughly how well one might expect simple grid search baseline to do with varying sample budgets. A well known fact about the order statistics of Gaussian random variables is the following:
3

---

Algorithm 1: A Local Bayesian Optimization Algorithm
Input: A black-box function, over convex compact set X.
for 1 = 1,2, , _o
X = argminz atrace (x,, Z) where Z E 1
// bt is the batch size evaluate the black-box function f on X, obtaining (possibly noisy) measurements // add yX.y) to the training data
D. = D yX)
X/+1 X2 Tie (UD,
end
Remark (e.g.. [8, 14]). Let X2... Xn be (possibly correlated) Gaussian random variables with marginal variance s², and let = max ( (i)). Then the expected maximum is bounded by
E S 2
This result directly implies upper bound on the expected maximum (or minimum by symmetry) of grid search with n samples. The logarithmic bound results in an optimistic estimate of the expected maximum (and the minimum): as n - oc the bound goes infinity as well, whereas the maximum (and the minimum) of GP is almost surely finite due to the correlation in the covariance.
With this, paths starting from x = from a centered Gaussian process with a RBF kernel (unit outputscale and unit lengthscale) using variety of observation noise standard deviations We then iterations of local Bayesian optimization (as described later in Algorithm 1) to convergence or until an evaluation budget of 5000 is reached. In the noiseless setting (o O), we modify Algorithm I to pass our gradient estimates to BFGS rather than applying the standard gradient update rule for efficiency.
evaluating the quality of local solutions. To do this, we optimize 50 sample
We repeat this procedure for d {1, 10, 20, 30, 50) dimensions and E {0, 0.05, 0.2) and display results in Figure 1. For each dimension, we plot the distribution of single local solutions found from each of the 50 trials local optimization (i.e,, without random restarts) is able to find a median objective value of 12299 corresponding to grid size of at least
box plot. In the noiseless setting, find that by d = 50, single run of
1036 achieve this value in expectation!
For completeness, we provide results of running GP-UCB and random search also with budgets of 5000 and confirm the effectiveness of local optimization in high dimensions.
4 A Local Bayesian Optimization Algorithm
The study above suggests that, under assumptions commonly used to theoretically study the perfor- mance of Bayesian optimization, finding natural next question is whether and how fast we can find them. To understand the convergence of local BO, describe prototypical local Bayesian optimization algorithm that will analyze in Algorithm 1, which is nearly identical to CIBO as described [19] with one exception: Algorithm allows varying batch sizes b, in different iterations.
single local solution is surprisingly effective. The
Define iruzz (X,, xt) (x,,,, x,) kD (x,, Z) [Z Z) + - xt). Namely, iruZZ : is the posterior covariance function conditioned the data D as well new inputs Z. Crucially, observe that the posterior covariance GP does not depends on the labels y. hence the compact notation kpuz. The acquisition function Otrrace of GIBO is defined
Oftrace (x, Z) ir V kp (x,, (x,, V
(2)
which is the trace of the posterior gradient covariance conditioned on the union the training data D and candidate designs Z.
The posterior gradient covariance Vkp((x, x,) V quantifies the uncertainty in the gradient fx(x), By minimizing the acquisition function Osrace (X,, Z) over Z. we actively search for designs Z that minimizes the one-step look-ahead uncertainty, where the uncertainty measured by the trace the posterior covariance matrix. After selecting set of locations Z, the algorithm queries the (noisy) function values Z, and updates the dataset. The algorithm then follows the (negative) expected gradient V UDD, (xx) to improve the current iterate X,, and the process repeats.

---

Convergence Results
In light of the strong motivation we now have for local Bayesian optimization, we will now establish convergence guarantees for the simple local Bayesian optimization routine outlined in the previous section, in both the noiseless and noisy setting. We first pause to establish some notations and few assumptions that will be used throughout our results.
Notation. Let X C Rdd be convex compact domain (e.g., X' [0, 119). Let H. be reproducing kernel Hilbert space (RKHS) on X' equipped with a reproducing kernel k(-,.) -). We use ||-|| to denote
the Euclidean norm and I| denote the RKHS norm. The minimum of continuous function defined on compact domain is denoted as f = minxex f(x).
In our convergence results, make the following assumption about the kernel k. Many commonly used kernels satisfy this assumption, e.g., the RBF kernel and Matérn kernel with v 2.
Assumption 1. The kernel k is stationary, four times continuously differentiable and positive definite. The following a technical assumption that simplifies proofs and convergence rates.
Assumption 2. The iterates {xt] stay in the interior of the domain X.
Intuitively, with the domain X' large enough and appropriate initialization, the iterates most likely always stay in the interior. In the case when the iterate X,+1 slides out of the domain X, we can add projection operator project it back to X' which will be discussed in $5.2.
We also define a particular (standard) smoothness measure that will be useful for proving convergence.
Definition (Smoothness). A function is T-smooth if and only for all X1 E X, we have:
- S L||x1 X2||-
Next, we define the notion of an ervor function, which quantifies the maximum achievable reduction in uncertainty about the gradient at x = 0 using set of b points Z and no other data, i.e., the data 'D an empty set in the acquisition function (2).
Definition 2 (Eroo function). Given input dimensionality d, kernel i and noise standard deviation we define the following error function:
Ed (b)
inf tr (V =
(3)
ZeRbxd
5.1 Convergence in the Noiseless Setting
We first prove the convergence of Algorithm 1 with noiseless observations, ie.., a 0. Where necessary, the inverse kernel matrix kp (Z, Z)-1 should be interpreted as the pseudoinverse. In this section, we will assume the ground truth function is a function in the RKHS with bounded norm. This assumption is standard in the literature [e.g.. 3. 24], and the results presented here ultimately extend trivially to the other assumption commonly made (that f a GP sample path).
Assumption 3. The ground truth function f is in H with bounded norm |f < B.
Because & is four times continuously differentiable, f is twice continuously differentiable. Further- more, on compact domain X, f is guaranteed to be L-smooth for some constant T. (see Lemma 9 for complete explanation). To prove the convergence of Algorithm 1, will show the posterior mean gradient VPD approximates the ground truth gradient V.. By using techniques meshless data approximation [e.g., 28, 5], we prove the following error bound, which precisely relates the gradient approximation error and the posterior covariance trace.
Lemma 1. For any E H, any E X and any D, we have the following inequality
[|V - VHD(x)]1 < tr (x) 11
(4)
Since / has bounded norm ||/|| < B, the right hand side of (4) is a multiple of the posterior covariance trace, which resembles the acquisition function (2). Indeed, the acquisition function (2) can be interpreted as minimizing the one-step look-ahead worst-case gradient estimation error the RKHS H. which justifies GIBD from a different perspective.
5

---

Next, we provide characterization of the error function Ed under the noiseless assumption o 0. This characterization will useful later to express the convergence rate.
Lemma 2. For = O, the error function is bounded by bb S C max 1 + d b}, where
8² C' = maxxisise Oz, k(0,0) 0) the maximum of the Hessian's 's diagonal entries at the origin.
Now we are ready to present the convergence rate of Algorithm under noiseless observations.
Theorem 1. Lei f & H whose smoothness constant is L. Running Algorithm with constant batch size b and step size Th. = for iterations outputs sequence satisfying
min
x (2L(f(x))) + B² Ed,
(5)
As T + co, the first term in the right hand side of (5) decays to zero, but the second term may not. Thus, with constant batch size b. Algorithm converges region where the squared gradient norm upper bounded by B² F.,,,o(b) with convergence speed O ((7). ). An important special case occurs when the batch size b, d in every iteration. In this case, F,, (b) 0 by Lemma 2 and the algorithm converges stationary point with rate (++).
Corollary 1. Under the same assumptions of Theorem 1, using batch size b, = d + 1, we have
min
[ ( f
SKT
Therefore, the total number of samples = (dT)) and the squared gradient norm [V converges to zero at the rate (dd/n).
We highlight that the rate is linear with respect to the input dimension d. Thus, as expected, local Bayesian optimization finds stationary point in the noiseless setting with significantly better non- asymptotic rates in d than global Bayesian optimization finds the global optimum. Of course, the fast convergence rate of Algorithm is the notable cost of not finding the global minimum. However, as Figure 1 demonstrates and as discussed in $3, single stationary point may already do very well under the assumptions used in both analyses.
We note that the results in this section can be extended to the GP sample path assumption as well. Due to space constraint, we defer the result to Theorem 4 in the appendix.
5.2 Convergence the Noisy Setting
We now turn our attention to the more challenging and interesting setting where the observation noise o >> O. We analyze convergence under the following GP sample path assumption, which is commonly used in the literature [e.g.. 24, 6]:
Assumption 4. The objective function f is a sample from GP, f GP k) with known mean function, covariance function. and hyperparameters. Observations of the objective function are made with added Gaussian noise of known variance, y = f(x) + E, where e (0,0a²)
Since the kernel k is four times continuously differentiable, f is almost surely smooth.
Lemma 3. For 0 < 8 1, there exists constant I 0 such that f L-smooth w.p. at least 1 6.
Challenges. The observation noise introduces two important challenges. First, by the GP sample path assumption, we have:
V (xxx) (VV (x:), VkD (x,, V.)
Thus, the posterior mean gradient V VHD, (x;) is approximation of the ground truth gradient V fxx)). where the approximation error is quantified by the posterior covariance Vkp, (x,, For any f. we emphasize that the posterior mean gradient VIDD (x.),(x)) is therefore biased whenever the posterior covariance is nonzero. Unfortunately, because of noise a > 0, this is always true for any finite data DD Thus, the standard analysis of stochastic gradient descent does not apply, as it typically assumes the stochastic gradient is unbiased.

---

The second challenge is that the noise directly makes numerical gradient estimation difficult. To build intuition, consider finite differencing rule that approximates the partial derivative
dr. 2h Øf 1 [f(x he;) f(x - he,)] + O(h))
where e; is the i-th standard unit vector. order to reduce the Ohh²) gradient estimation error, we need h + 0. However, the same estimator under the noisy setting
as - (ex) + EI) (f(x he,) + E2) + O(h²), 2h
where E1, E2 are ii.d. Gaussians, diverges to infinity when h + 0. Note that the above estimator is biased even with repeated function evaluations because cannot go zero.
We first present a general convergence rate for Algorithm that bounds the gradient norm.
Theorem 2. For 0 ô <, suppose f is a GP sample whose smoothness constant is L w.p. least 1- 6. Algorithm 1 with batch size by and step size Ttt = 7, produces sequence satisfying
Vf(xx)11 & ((2[(f(x)) f*) + + C,Eax, ISST
min
(6)
with probability at least 1 26, where C2 2log( (1o (t) /6))
1 The second term 7. Ei= C.E (b.) in the right hand side of (6) is the average cumulative bias of the gradient. To finish the convergence analysis, must further bound the error function Ed, (bi). For the RBF kernel, we obtain the following bound:
Lemma 4 (RBF Kernel), Let k(x), = eep((=2lxx - x2ll²) be the RBF kernel. We have
(2md)) S d 1 W
(01++(W- (-e(mm o2)
m
)) O(adm +)),
where E N and W denotes the principal branch of the Lambert W function.
The error function (3) is infimum all possible designs, which intractable for analysis. Instead, we analyze the infimum over a subset of designs of particular patterns (based on finite differencing), which can be solved analytically, resulting the first inequality. The second equality is proved by further bounding the Lambert function by its Taylor expansion at -1/e.
In addition, obtain similar bound for the Matérn kernel with 2 by slightly different proof.
Lemma (Matern Kernel), Let k(.,)) be the
Matérn kernel. Then, we have
dm o S adm 1 + o dm 1om-i 1 O(odm= +).
Interestingly, Lemma 4 and Lemma 5 end up with the same asymptotic rate. Writing the bound in terms of the batch size b, we can see that E (6) b à) for both the RBF kernel and the 2.5 Matérn kernel. Coupled with Theorem 2, the above lemmas translate into the following convergence rates, depending on the batch size b, in each iteration:
Corollary 2. Let k(.,-) -) be either the RBF kernel the 2.5 Matérn kernel. Under the same conditions as Theorem if
d logg t;
O(I/T) + (Odd); O log T) O n log n); O O log2 T) O log" n))
b. =
dt; di
then
min ISSS
[Vx)x
with probability at least 1 - 26. Here, T is the total number of iterations and is the total number of samples queried.
We suspect similar bound holds for the entire Matérn family.
7

---

With nearly constant batch size b, dlog" t. Algorithm 1 converges to region where the squared gradient norm is ((ad), With linearly increasing batch size, the algorithm converges stationary point with rate 0(ad1125n-0022 log n), significantly slower than the (ddnn) rate in the noiseless setting. The quadratic batch size is nearly optimal up to a logarithm factor increasing the batch size further slows down the rate (see Appendix D for details).
To achieve convergence to a stationary point using Corollary the batch size must increase as optimization progresses. We note this may not be an artifact of our theory or our specific realization of local BO, but rather general fact that is likely to be implicitly true for any local BO routine. For example, a practical implementation of Algorithm 1 might use a constant batch size and a line search subroutine where the iterate x; is updated only when the (noisy) objective value decreases otherwise, the iterate X, does not change and the algorithm queries more candidates to reduce the bias in the gradient estimate. With this implementation, the batch size is increased repeatedly on any iterate while the line search condition fails. As the algorithm converges towards a stationary point, the norm of the ground-trut gradient | fx+)) decreases and thus requires more accurate gradient estimates. Therefore, the line search condition may fail more frequently with the constant batch size, and the effective batch size increases implicitly.
We provide two additional remarks on convergence in the noisy setting. First, the convergence rate significantly slower than in the noiseless setting, highlighting the difficulty presented by noise. Second, when the noise is small, the convergence rate is faster. If o + O, the rate is dominated by lower-order term in the big O notation, recovering the ((+) rate in the noiseless setting (see Appendix D for details). This is in sharp contrast with the existing analysis of BO algorithms. Existing convergence proofs in the noisy setting rely on analyzing the information gain, which is vacuous when o + 0, requiring new tools. It is interesting that no new tools are required here.
Finally, we revisit Assumption 2. In the case when it does not hold, we use modified update:
X++1 = projx(xx - 178 VyD, (x:)).
where the projection operator projæ (.) is defined as proj (x) = argmingexI1z -x|,, i.e., the closest feasible point to x. When the iterates stay in the interior of the domain X, the projection operator is an identity map and the update rule simply reduces to X++1 X, - mh VHD,
Define the gradient mapping
G(x,) 1 (x: projx Ne fxx))),
which is generalization of the gradient in the constrained setting: when x 1), fxx:) lies in the interior of X' the gradient mapping reduces to the (usual) gradient (xx,) The following gives convergence rates of [|G(x))]l-
Theorem 3. Under the same conditions as Corollary 2. without Assumption 2, using the projected update rule (7), Algorithm 1 obtains the following rates:
dt; dt?
(0 O(odán log n + of atd log n); (odd (0 n log² n + a d - logn),
if b. with probability at least 1 25. Here, n is the total number of samples queried.
then
min 151ST
[xn_11((x))1==
These rates are slower than Corollary 2. We defer more details to Appendix C.
Additional Experiments
In this section, investigate numerically the bounds in our proofs and study situations where the assumptions are violated. We focus on analytical experiments, because the excellent empirical performance of local BO methods high dimensional real world problems has been well established in prior work [e.g., 10, 19, 20, 27]. Detailed settings and additional experiments are available in Appendix E. The code is available https:/ //gi /gi thub. woon/kayween/non lorieeeen/oollllleeee bo- convergence.
6.1 How loose are our convergence rates?
This section investigates the tightness of the bounds on the error function F.g,, - key quantity in our convergence rates. We plot in Figure 2 the error function Ed for the = 2.5 Matérn
8

---

10
10a
bound
o=
bound
6=0.1
error function Eako (b)
err. func. d=0.1
error function (b)
err. func. 0=0.1
10
10"
10
10-1
10-1
bound 0.05
bound 0.05
10-2
err. func, o= 0.05
err. func. a- 0.05
10a
10a
10
10P
10"
10
batch size b
dimension
(a) lixed dimension, varying batch sizes
(b) fixed batch size, varying dimensions
Figure 2: Compare the error function (3) of the v 2.5 Matérn kernel and our upper bound in Lemma 5. The error function (d) is approximated by minimizing (3) with L-BFGS. Both plots are in log-log scale. Left: The slope indicates the exponent on b. Since the slope magnitude of the error function is slightly larger, the error function might decreases slightly faster than ((b=z) asymptotically. Right: The slope indicates the exponent on d. Since all lines have roughly the same slope, the dependency on the dimension in Lemma seems to be tight.
probability density
best value found
-19
-20
-20
19
10
60
80
100
local solution value
repeats
Figure 3: The performance of random restart on GP sample path in 100 dimensions. Left: a density plot for the minimum value found single restart (compare with Figure 1). Right: the median and a 90% confidence interval for the best value found after given number of random restarts.
kernel and our bound O(odàb-)) implied by Lemma The error function is empirically estimated by minimizing the trace of posterior gradient covariance (3) using L-BFGS. Figure 2a uses fixed dimension d 10 and varying batch sizes, while Figure 2b uses fixed batch size 500 and varying dimensions. Both plots are in log-log scale. Therefore, the slope corresponds to the exponents on b and d in the bound O b- ), and the vertical intercept corresponds to the constant factor in the big notation.
Figure 2a shows that the actual decay rate of the error function (t),, may be slightly faster than (( 3). the slopes of the (empirical) error function have magnitude slightly larger than / Interestingly, Figure 2b demonstrates that the dependency d on the dimension is quite tight all lines in this plot share similar slope magnitude.
6.2 What is the effect of multiple restarts?
In $3, we analyze local solutions found by a single run on different GP sample paths. Here, we investigate the impact of performing mulriple restarts on the same sample path. In Figure 3, we plot (left) kernel density estimate of the local solution found for series of 10000 random restarts, and (right) the best value found after several restarts with a 90% confidence interval. We make two observations. First, the improvement of running 10-20 restarts over single restart is still significant: the Gaussian tails involved here render difference of +1 in objective value relatively large. Second, the improvement of multiple restarts saturates relatively quickly. This matches empirical observations

---

2.8
2.8
18
24
24
1.0
15
1.0
2.0
1.0
2.0
0.
12
1.6
0.3
1.6
0.9
1.2
1.2
0.6
subgradient data
0.0
subgradient data
0.8
data
005
0.3
-0.5
0.4
-0.5
-1.0
00.5
0.0
0.5
1.0
-1.0
0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
(a) (x) (0.039, 1.018).
(b) Vpp(x)= (0.005, 0.318).
(c) VDp(() = (-0.031, 1.053). Subgradient [-1, 1] {1}-
The ground truth (0, 1).
Subgradient [-1,
Figure 4: Estimating the "derivative" at x (0, 1) with Matérn Gausssian process (v 2.5) in samples, the posterior mean gradient is close to the ground truth. Middle: f(x) = ||x||1 The & norm is not differentiable at (0, 1). With exactly the same samples as the left panel, the posterior mean gradient has higher error. Right: f(x) = |||||_ Increasing the sample size to = 10 decreases the estimation error.
three different settings. Left: [(x 2 ||| With
made when using local BO methods on real world problems, where optimization is often terminated after (at most) a handful of restarts.
6.3 What if the objective function non-differentiable?
In this section, we investigate what happens when f is not differentiable setting ruled out in our theory by Assumption 1. In this case, what does the posterior mean gradient ViDD learn from the oracle queries? To answer this question, we consider the & norm function x +> ||x|| in two dimensions, which is oin-differentiable when either 2 or 22 zero. The l norm is convex and thus its subdifferential is well-defined:
0||||| {g IR; |l 1, ||x|11]
In Figure 4, we use the posterior mean gradient Vμr to estimate the (sub)gradient at (0, 1). As comparison, we also show the result for the differentiable quadratic function + ||x||22 Comparing Figure 4a and 4b, we observe that the é norm has higher estimation error than the Euclidean norm despite using exactly the same set of queries. This might suggest that non-differentiabiiiyy increases the sample complexity of learning the "slope" of the function. Increasing the sample size to n = 10 eventually results in Vud (x) much closer to the subgradient, as shown in Figure 4c.
7 Discussion and Open Questions
For prototypical functions such as GP sample paths with RBF kernel, the local solutions discovered by local Bayesian optimization routines in high dimension are of surprisingly high quality. This motivates the theoretical study of these routines, which to date has been somewhat neglected. Here have established the convergence of recently proposed local BO algorithm, GIBO, in both the noiseless and noisy settings. The convergence rates are polynomial in both the number of samples and the input dimension, supporting the obvious intuitive observation that finding local optima is easier than finding the global optimum.
Our developments in this work leave a number of open questions. It is not clear whether solutions from local optimization perform so well because the landscape is inherently casy (e.g., most stationary points good approximation of the global minimum) or because the local BO has unknown algorithmic bias that helps it avoid bad stationary points. This question calls for analysis of the landscape of GP sample paths (or the RKHS). Additionally, we conjecture that our convergence rates are not yet tight; $6.1 suggests that there is likely room for improvement. Further, it would be interesting to establish convergence for trust-region based local BO [e.g. 10, 9].
In practice, while existing lower bounds imply that any algorithm seeking to use local BO as a subroutine to discover the global optimum will ultimately face the same exponential-in-d sample complexities as other methods, our results strongly suggest that, indeed as explored empirically in the literature, local solutions can not only be found, but can be surprisingly competitive.
10

---

Acknowledgments and Disclosure of Funding
KW would like to thank Thomas T.C.K. Zhang and Xinran Zhu for discussions on RKHS. RG is supported by NSF awards OAC-21 18201, OAC-1940224 and IIS-1845434, JRG is supported by NSF award IIS-2145644.
References
[1] Ahmad Ajallocian and Sebastian U Stich. On the convergence of SGD with biased gradients. arXiv preprint arXiv:2008.0000 2020.
[2] Sébastien Bubeck et al. Convex optimization: Algorithms and complexity. Foundations and TrendsO in Machine Learning, 8(3-4):231-357,20) 2015.
[3] Adam D. Bull. Convergence rates of efficient global optimization algorithms. Journal of Machine Learning Research, ((88)):8779-900) 2011.
[4] loannis Chatzigeorgiou. Bounds on the Lambert function and their application to the outage analysis of user cooperation. IEEE Communications Letters, 17(8):1505-15000 2013.
[5] Oleg Davydov and Robert Schaback. Error bounds for kernel-based numerical differentiation. Numerische Mathematik, 132:243-269, 2016.
[6] Nando De Freitas, Alex Smola, and Masrour Zoghi. Exponential regret bounds for Gaussian process bandits with deterministic observations. In Proceedings of the 29th International Conference on Machine Learning, 2012.
[7] Aryan Deshwal and Jana Doppa. Combining latent space and structured kernels for Bayesian optimization over combinatorial spaces. In Advances in Neural Information Processing Systems, volume 34, pages 8185-8200, 2021.
[8] Luc Devroye and Gábor Lugosi. Combinatorial Methods in Density Estimation. Springer Science Business Media, 2001.
[9] David Eriksson and Matthias Poloczek. Scalable constrained Bayesian optimization. In Proceedings of The 24th International Conference Artificial Intelligence and Statistics, volume 130, pages 730-738. PMLR, 2021.
[10] David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local Bayesian optimization. In Advances in Neural Information Processing Systems, volume 32, 2019,
[11] Roman Garnett. Bayesian Optimization. Cambridge University Press, 2023.
[12] Antoine Grosnit, Rasul Tutunov, Alexandre Max Maraval, Ryan-Rhye Griffiths, Alexander I Cowen-Rivers, Lin Yang, Lin Zhu, Wenlong Lyu, Zhitang Chen, Jun Wang. et al. High- dimensional Bayesian optimisation with variational autoencoders and deep metric learning. arXiv preprint arXiv:2106. 106.03609, 2021.
[13] Rafael Commez-Bombarelf Jennifer Wei, David Duvenaud, José Miguel Hernández-Lobatt Benjamín Sännchez-Lengeling Dennis Sheberla, Jorge Aguilera-Iparragee Timothy D. Hirzel, Ryan Adams, and Alán Aspuru-Guzik. Automatic chemical design using a data-driven continuous representation of molecules. ACS Central Science, 4(2):268-22) 226,8 2018.
[14] Gautam Kamath. Bounds the expectation of the maximum of samples from Gaussian. 2015. URL http: //www/ wtt.://aattttaaaaaaaaa.. com weam/writiin www.tiigst/i/eeeeeeeeei /gaps//gsssiiis an_maxx____aaaaaaa pdf.
[15] Kirthevasan Kandasamy, Jeff Schneider, and Barnabas Poczos. High dimensional Bayesian opti- misation and bandits via additive models. In Proceedings of the 32nd International Conference on Machine Learning, volume 37, pages 295-304, 2015.
[16] Kenji Kawaguchi, Leslie Pack Kaelbling, and Tomás Lozano- -érez. Bayesian optimization with exponential convergence. In Advances Neural Information Processing Systems, volume 28, 2015.
[17] Natalie Maus, Haydn Jones, Juston Moore, Matt J Kusner, John Bradshaw, and Jacob Gard- ner. Local latent space Bayesian optimization over structured inputs. Advances in Neural Information Processing Systems, volume 35, pages 34505-34518, 2022.

---

[18] Natalie Maus, Kaiwen Wu, David Eriksson, and Jacob Gardner. Discovering many diverse solutions with Bayesian optimization. In Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, volume 206, pages 1779-1798, 2023.
[19] Sarah Müller, Alexander von Rohr, and Sebastian Trimpe. Local policy search with Bayesian optimization. In Advances in Neural Information Processing Systems, volume 34, pages 20708-207200 2021.
[20] Quan Nguyen, Kaiwen Wu, Jacob Gardner, and Roman Garnett. Local Bayesian optimization via maximizing probability of descent. In Advances in Neural Information Processing Systems, volume 35, pages 13190-13202, 2022.
[21] Jonathan Searlett. Tight regret bounds Bayesian optimization one dimension. In Proceed- ings of the 35th International Conference on Machine Learning, volume 80, pages 4500-4508, 2018.
[22] Yao Shu, Zhongxiang Dai, Weicong Sng, Arun Verma, Patrick Jaillet, and Bryan Kian Hsiang Low. Zeroth-order optimization with trajectory derivative estimation. In The Eleventh International Conference on Learning Representations, 2023.
[23] Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical Bayesian optimization of machine learning algorithms. In Advances in Newral Information Processing Systems, volume 25, 2012.
[24] 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 on Machine Learning, pages 1015-1022, 2010.
[25] Samuel Stanton, Wesley Maddox, Nate Gruver, Phillip Maffettone, Emily Delaney, Peyton Greenside, and Andrew Gordon Wilson. Accelerating Bayesian optimization for biological se- quence design with denoising autoencoders. Proceedings of the 39th International Conference
on Machine Learning, volume 162, pages 20459-20478, 2022.
[26] Austin Tripp, Erik Daxberger, and José Miguel Hernández-Lobato Sample-efficient optimiza- tion in the latent space of deep generative models via weighted retraining. In Advances in Neural Information Processing Systems, volume 33, pages 11259-11272 2020.
[27] Linnan Wang. Rodrigo Fonseca, and Yuandong Tian. Learning search space partition for black-box optimization using Monte Carlo tree search. In Advances in Neural Information Processing Systems, volume 33, pages 19511-19522, 2020.
[28] Holger Wendland. Scattered Data Approximation, volume 17, Cambridge University Press, 2004.
[29] James T Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Peter Deisenroth. Pathwise conditioning of Gaussian processes. The Journal of Machine Learning Research, 22(1):4741-478 2021.
12

---

A Technical Lemmas
Lemma 6. Lef (2) Eizo 2or be the degree-n Taylor polynomial of e around x: O. Then for any h 2 and any I E RR we have
2 (2)-
Proof. The proof is an induction on k. The base case k: = is trivial: e > 1 x.
Consider the general case 1- Define g(x) ex p (x). It easy to that g (a) e" (2) and g"(x) e" (2 By the induction hypothesis, g" 2 0 and therefore is convex. Thus, the minimum of g is given by its stationary points. easy to observe that
indeed stationary point. Thus, min g(x) = g(0) = 0, which finishes the proof.
Lemma 7. Let o > 0 and W be the principal brunch of the Lambert W function. For any m E N, have
20
1 W(
)
e(m + r2)
+ 02
Proof. The problem reduced to proving
1 W(x) < 2 er)
for - S x S 0. In fact, the right hand side is exactly the first-order Taylor expansion of the left hand side at x = [e.g.. 4].
Square both sides. It is equivalent to prove (1 +1+ ((x))² < 2(1 + +ex for æ E [- Define g(x) = (1 + ((x)) (1)
+z).
The derivative of g in (-1 is
( = ((x)
By the definition of the Lambert W function, we have
W(x)e < W(x)
since - < ((x) for x E (- 0). Thus, gxxx for x (- 0) and g is decreasing
function in |- -- 0]. Observe that g = Therefore g(ø) S 0 for all x [- 0j, which
completes the proof.
Next, we present lemma which states that the error function bounds posterior covariance trace in each iteration of Algorithm 1.
Lemma 8. In the t-th iteration of Algorithm we have
Lr LrkkD, (X1, xx) VT) S (ba),o(b)))
Proof. Without loss of generality, we assume X₃ = O. Otherwise, shift the data D; and x; by -Xt- which does not change the value of the left hand side becauce stationarity of the kernel Let Z E R bb be arbitrary candidates. Then, we have
tr(Vko, Luz (0, O) < (r(Vk(0, T
Vk(O0, (), 2oo O)V
Because the LHS conditions on both D.-1 and Z but the RHS only conditions on Z. Now, minimize on both sides.
On the one hand, the LHS becomes tr (Vkp, (0, 0)V This is because D, is the union of D,-1 and the minimizer of the acquisition function Ottrace (0, Z) tr(VkD,. 0)V T).
On the other hand, the RHS becomes E (b.) by definition of the error function, which completes the proof.
Note that the edge case, where the minimizer the acquisition function argminz Otrace (0, Z) does not exist (e.g., when 0), can be handled by careful limiting argument using the same idea. D
13

---

A.1 Lemmas for the RKHS Assumption
By Assumption the kernel function is four times continuously differentiable. The following lemma asserts the smoothness of E H.
Lemma 9, Suppose f E maps from a compact domain X' R. Then f is I--moooh for some
Since the kernel k is four times continuously differentiable, f is twice continuously differen- tiable. On a compact domain X, the spectral norm of the Hessian [[V]())]] has a maximizer. Define L = maxxex [x)]l- Then f is L-smooth.
Proof.
Lemma 1. For any f E H. any x E X and any D, we have the following inequality
VV - [up(x)]| & tr (Vkp(x, V))
(4)
Proof. This is simple corollary of standard result in meshless scattered data approximation [e.g-- 28, Theorem 11.4]. The main idea is to express the estimation error as a linear functional, and then compute the operator norm of that linear functional.
Let = Dx : H + be the composition the evaluation operator and differential operator, i.e. \ʃ (x).). Wendland [28, Theorem 11.4] provides bound on (λʃ \pp)
(Af - >pd) 2 s (2) (1l0)))111
where (1) applies À to the first argument of kD and ( applies À to the second argument of kp. Pick the linear functional À f 2 f(x) where 1 S S d. Then, the left hand side becomes @
(Af - a (x)
b2 oz (x))) The right hand side ( (-)(s) is exactly the i-th diagonal entry of VkD (x, x)V V. For each i, use the above inequality, and then summing all coordinates finishes the proof.
A.2 Lemmas for Convergence on GP Sample Paths
In this section, we provide few lemmas for the GP sample path f ~ GP(O, k). By Assumption 1, we have the follow lemma which asserts that f smooth with high probability.
Lemma 3. For O < 8 < 1, there exists constant I, >> O such that I is L-smooth w.p. at least 1 8.
Proof. The proof uses the Borell-TIS inequality, Let f~ CP(O, k) be a Gaussian process. Provided that supxex (x) is almost surely finite, the Borell-TIS inequality states that
u
Pr(sup fx(x) E sup f(x))) I exp(- xEX
2s²
xEX
where u 0 is arbitrary positive constant and s² supxex [x) Namely, if the supremum of f is almost surely finite then the supremum of f is bounded with high probability.
Since k is four time continuously diffcrentiable, the second-order derivative
dz,dz, f exists and is
8²
almost surely continuous. On compact domain X,, the supremum suPxxxx 0z,02, fxx) is almost surely finite. By the Borell-TIS inequality, the expectation supxxx lā, [x) finite and the
2²
a²
supremum supxex dr;do, f(x) is bounded with high probability. Thus, the Frobenius norm of the Hessian V fx)) is bounded with high probability by a union bound. Since the spectral norm of (x(x) can be bounded by its Frobenius norm, the spectral norm of the Hessian [x)]] is also bounded with high probability, which gives the smoothness constant.
d²
The next two lemmas bound the gradient estimation entor with high probability.
Lemma 10. Let u ~ (0, E) be a Gaussian vector. Then for any >
t²
Pr[[|u[| > t) S 2 2exp
2
14

---

Proof. This is a standard concentration inequality result, but we give self-contained proof here for completeness. Denote the spectral decomposition QAQ". By Markov's inequality, we have
Pr(I|u] > t) = Pr(ee > e*t)
S < Belllull
where is arbitrary positive constant. Let c (0,0 I) be standard Gaussian variable. Then is easy to see that u|] = ||A Then, replacing |ul with | e gives
Pr([|u|| > 1) S e Be A ttll
A Aeel
S e st Be
d
= =e
st
Ee
i=1
2e -st
IE EVa
1=1
2e at 202 Ei A
where the first line plugs in the second line uses || 112 S || I1 the third line is due to independence of e;; the forth line removes the absolute value resulting an extra factor of the last lines uses the moment generating function of E,. Optimizing the bound over gives the desired result
Pr([|u|| > t) 5 2e
2 b
2exp(
2 Σ
Lemma 11. For any 0 <6<,, ler C. = 2 2log( 66
=
Then, the inequalities
[v,();) VuD, (xx) S C. tr Vhr, (x,, x,) V
hold for any 2 with probability least 8.
Proof. Since f((x)) ~ NNVV (x:), VkD, (x,, (x. V applying Lemma 10 gives
1
Pr(IVf((,) VHD, (x:) C. tr(Vko (,, x, [T)) s 2 cxp( Ct).
The particular choice of Ct makes the probability the right hand side become 68 union bound over all 2 1 and using the infinite sum Σ
Using the
x² finishes the proof.
n
We provide an important remark. The probability in Lemma 11 is taken over the randomness of and the observation noise. On the other hand, the posterior mean gradient Vup, is deterministic, since it is conditioned on the data Dr-
B Bounds the Error Function Eak..
This section is devoted to bounding the error function (a) in terms of the batch size b. The results in this section immediately give a bound on the posterior covariance trace by Lemma 8.
Before diving into the proofs, we present some immediate corollaries of Assumption I on the kernel. Because is stationary, the kernel can be written as k(x, x') o(x x') for some positive- definite function d. Observe that Vk(x, Denote the first-order partial derivative 0x Bag (xx and the second-order partial derivative a² o(o(x) az? (x). It is easy that is function and V is odd function. In addition, 0 is maximum of Φ. Therefore, V (0) 0 and the Hessian 2 (00) is negative semi-definite.
Vo(x x') and Vk(x, x')
V p x').
15

---

B.1 Noiseless Setting
The following a bound for the enror function Ed, for arbitrary kernels satisfying Assumption in the noiseless setting d = O.
Lemma 2. For O, the error function bounded by Ed S C max (0, + d b), where
8² C' = sigd Ta, dw,ō (0, 0) is the maximum of the Hessian's 's diagonal entries at the origin.
Proof. The bound holds trivially for 0, and thus proof is only needed for b 2 2, which we split into two cases < f d + 1 and > + 1.
We first focus on the case 2 < b < d + 1. Let Zo = 0 and = he, where i = 1,2, ... 1, where e, is the i-th standard unit vector and the definition of the error function, we have
is constant. Define (zo 21 26-1) T. By
(b) & (V((0,0)) -
A
b-1
C(1 + d - -b)
Σ
An
where define A Vk(0, 0) (Z) -0) S C for < i S d.
and use the inequality An s
Let us focus on the i-th diagonal entry A; where 1 sis b- - 1. Then we have
A# - (0) - (0 0,0(-hes) (
(0)
o(he,))
0
p(he,) 600)
(2,((hee))
-7)(00) (Φ))) - (o(he,))
1
(0
(-o()))((h))))
6(0) (hhh,)
-000))
(ao((noa))
x(0))
= =0,0(00) (0)0,,(())))
(o(he))(((h)))))
where the first line is because conditioning on the subset Zo and z, does not make the posterior smaller. Now let + and compute the limit by L' L'Hôpital's rule. We have
(()h(e))))
lim (h) = lim h
0g ((0)
h-+
(á(0))2 - (o(he,))2
((0)
(ho(hee)0))(() -,,(he)
lim h=0
- x0)(0) (00) + (he;;)
Thus letting -> gives the inequality E. (b) S C(1 + - b) for 2 S 6 d
When d >> + 1, note that (b) is an decreasing function in b and thus Fa (b) S F. (d
1) 0. Both cases can be bounded by the expression C maax (0, + d b).
B.2
Noisy Setting
This section proves bounds on the error function EaR..r for the RBF kernel and the v = 3 Matérn kernel in the noisy setting. The lemmas in this section will implicitly use the assumption that k(0, 0) = 1. This assumption is indeed satisfied by RBF kernel and the Matérn kernel, which are of primary concern in this paper.
Before proving the bound on Ed we need one more technical lemma:
Lemma 12 (Central Differencing Designs). Consider the 2md points Z B 2mdxd defined
(t)
-hei, 11,,2, ... 17
hei j j= m - 1, 1,m- + 2 .. 2m.
16

---

where 1
S 2m and e; is the i-th standard unit vector. Define
A = Vk(0,0)V (, (k; + +
Then, we have
2
Au S -80((0) (1 ac) +
for all 1 sis and thus
233
tr (r) S Σ (0f6(0) (1- a;) +7
where a; @(2he;), 8r = 0,0(-hee) and d²
Proof. Note that A is the posterior covariance the origin 0 conditioned Z. Denote Z the subset of 2m points that lie on the i-th axis. Then,
(1)
(i)
2
22
m+1
2m2m
the i-th diagonal entry Au can be bounded by
Au -ff 0, Z (k(( Z + - T
since conditioning on only a subset of points z() would not make the posterior variances smaller (e.g., the posterior covariance is the Schur complement of positive definite matrix). The remaining proof is dedicated to bounding the right hand side.
First, we need to compute the inverse of the 2m x 2m kernel matrix
K
(8)
Z
( + O²1
11 O; 1,111
+
a,11
11
where =<(hee; is nonnegative constant and 1 is m dimensional vector (we drop the index i in the matrix K for notation simplicity). We compute the inverse analytically by forming its eigendecomposition
K = QAQ
where A= = diag(A, A1, "., .. A2m) and Q = (q)
(2m). Observe that:
K
1 a
11 + +
where & denotes the Kronecker product. Because the eigenvalues (vectors) of Kronecker product equal the Kronecker product of the individual eigenvalues (vectors), and because adding diagonal shift simply shifts the eigenvalues, the top two eigenvalues of K are a₁ m(1 a:) + o² and d2 = m(1 - a;) + o The remaining 2m 2 eigenvalues are a². The top two eigenvectors are
(1)
1 (1))
q1
22m
2mm
Next, we cope with the term 0,k(0, ()), where the partial derivative taken w.r.t. the first argument's s i-th coordinate. Denote v = (0, Z ) Then it easy to see that
v (
1
where 3; = =0,ō(-he,) Note that v happens to be an eigenvector of K as well, because v I q2- As
a result, Q V has a simple expression Q (0 22mm, 0 0) Thus, we have
A -0)((0) v QA -
2m.3?
-8,6(0) m(1 a;) + +
287
-0)((0) (1 - ax) +
Summing over the coordinates < i < d finishes the proof.
17

---

With Lemma 12, we are finally ready to present the bounds for the RBF kernel and the Matérn kernel. Lemma (RBF Kernel), Let x(X1,XX) = exp(=
-- [|| - x2ll²) be the RBF kernel. We have
B (2md) &
- (
m
))
(oddm b);
o²)
where mn E N and W denotes the principal brunch of the Lambert W function.
Proof. For any Z E R2ndxd the following inequality holds by the definition of error function
B (b) S (r(Vkk0,00) - Vk(00, 7
Consider the following points
( ) defined as
( hei, j= 2.
he,, j == m - 1,m + .. 2m,
where isis The total number of points is exactly 2md. By Lemma 12, we have
2/9 Ed,k,o (b) S (0)((0) + (1 - at) +7
d
For the RBF kernel, the values of a; and Bi the for each coordinate isotropic:
sis S dl since is
0; a= ((2hcc = cxp(-2h²). == 0;0( (- =2)
-xp
Plugging the value of a, B, y and (0) into the bound on Ed,kios we have
2m (ineexpphh?)
Aa S
mm1 expl
-2)²2) + o²
2m exp(
222
< m expl
22h2) o²
where the second inequality replaces exp(- with (-p(-212) in the numerator. Because the
bound holds for arbitrary h, we apply the transformation h Ja h, which gives the inequality
exp(-h3)) m(1 exp(-h?)) 0²
A S
Our goal is to bound Ag terms of and o². Therefore, we minimize the right hand side over h. me Define g(z) = m(1- =)- where x 2 0. The derivative is given by:
m(m (m + 0²) 0²)e" (x-1)) 1))
g(xx) = (m(e= - 02e=)2
The unique stationary point is x = 1 W( (- where W is the principal branch of the Lambert W function. It is easy
the stationary point x" is the global minimizer of g(x) over Rzo. Plug æt into the expression of g. Coincidentally, we have g(x*) = 1 W(- c(m++o?) (=e(mm+ as well x² is fixed point of g.
"((++o=]) for each coordinate Summing As over In summary, we have shown A. 1 W (V(--ntfor)
all coordinates proves the first inequality (2mx) d + W( (((=((faa))) The second inequality is a direction implication of Lemma 7, which completes the proof.
Lemma (Matern Kernel), Let k(..-) be the v = 2.5 Matérn kernel. Then, we have (akkoo2md) S adm σ dm & Olodm=
18

---

Proof. The proof is similar to Lemma 4. The difference is that we need to upper bound 00) - 28
(1-0:)+y by a rational function. Otherwise the expression is intractable to minimize.
The Matérn kernel with half integer v can be written as a product of an exponential and a polynomial. In partieular, for = we have
(xx) = (1+ + v5| | +
exp(
When x is in the nonnegative orthant, the gradient is
v5 x
10 x)
v5
exp(- exp(-v5| +1 + 5 + =1
V (x(x) = exp( ( x
||x||
+ /
3
Since the Matérn kernel is isotropic, the as and B3 in Lemma 12 the same across different coordinate i, and their values are
20
a = a; = 6(2he,) = exp((22/5h) 2v5h +
3
5 exp(- (-v v5h)h,
B 3 0,0(- -0,ō(he,,)
3
In addition, -876(0)
3. By Lemma 12, we have
2,32
Ag S -85(10) (1 (1-a) +
10
+ h)
3(3 exp( 23 + 6h + 4h2)) + 97
Next, we approximate the exponential function exp(-2h) by its Taylor polynomials. By Lemma 6, use the inequality cxp(-2h) 2 1 2h for the numerator and the inequality cxp( -22) 2 2h
2h h for the denominator. Applying these two inequalities gives
5 10(1 21)(1 + h)2
Aa S
23 8h²) + 9y
Let - and thus 7 = h4, Then we have
5 10(1 - 2h) + +
Aa S
3
2h²(3 + 813) 9h
5h²(27 + 28h)
3(6 + 9/2 16h²)
112 28h)
18
15
70 a
15
70 2 m 1
om
where the third line drops 1 and 1 the denominator; the forth line plugs in the value h 7 back and drops the constants. Summing over the coordinates 1 S gives the first inequality:
Ea,k. (2md) S ddm 1 + o dm
For large enough m, the bound is dominated by the first term odm- 1 which completes the proof.
We end this section with a short summary. Lemma and Lemma 5 happen to end up with the same rate Ea (2md) = O(adm )). Replacing 2md with the batch size b, we have shown that (b) a b-1) for the RBF kernel and y 2.5 Matérn kernel.
19

---

B.3 Discussion: Forward Differencing Designs
This section explores an alternative proof for the error function based on forward differencing designs, as opposed to the central differencing designs in Lemma 12. Similar to the previous section, we assume (0, 0) 1, which is indeed satisfied by the RBF kernel and Matérn kernel.
Lemma 13. Consider following (d + 1)m points Z E [Rd+1)mxd defined as
(0)
2
0,
j == 1,2,
(3) hei i2 1, j j= 1.2 ,.,
2
where e; is the i-th standard unit vector. Define
A = Vk(0,0)J - Vk(0, (k + (A ( 0)
Then, we have
(1 (1+
tr < - Σ 030(0) +
(1 + 7)2 a7
where @(he,), 3, = (he) and o²
Proof. The proof is similar to Lemma 12.
Denote (6) be the subset of points consisting of z and 2 where j j= 1,2, m. Notice that the i-th diagonal entry Ay can be bounded by
Au < - o. 0, 2 (k:(())) (1) + (O.K(,
We need invert the 2m 2m matrix
K = k(Z (6) 2 + +
11 T
a.l
+ -
051 11
Again, we resort to the eigendecomposition of K QAQT. The top two eigenvalues of K are
A = m1 + a;) - d² and 22 = m(1 - az)) + a². The remaining 2m - 2 eigenvalues are σ². The
top two eigenvectors are
1
V2m (4)
q1
q2
22m
Denote T = -0, ) Note that v can be written as linear combination of qi and q2:
0 = マ 3.v 3.V2mlqq -q2)-
v v= 3
-1 0 O)T. Then, straightforward calculation gives Q fBv 33.22m 2m
Then, we have
Aa < -8t0(0) v QA QTv
-70(())
1 2 m
m(1 + as) 02 +
1
1
m(1 as) 0
-P0(0) I 32
1
1+aity + 1- a; +7
8 + 7) (1 + 7)2 a
-P0(0)
Summing over all coordinates finishes the proof.
20

---

Lemma 14. Let x2) = exp [1 -x2 be the RBF kernel. The forward differencing designs give a decay rate of o. O :).
Proof. Deñie = <(he;) and B (:0(-hee; Plugging the values of a and 3 into the bounds in Lemma 13 yields
122 exp((11)) + 7)
Au S 1
(1 -1xp))-exp(-b3)
The bound holds for arbitrary h. Applying the transformation h + vh gives
A S1 hc +
(1+ exp(-h)
h(1 hh
S 1
(1+ 7)² (1- h)
h(1 - h) +
1
72 + 27 h
h(1 - h)
72 + 2y + h
where the second line uses the inequality exp(-h) 2 h in the numerator and the denominator; the last line is because y is nonnegative. Let h y so that y h2. Then we have
h(1 h)
Au S
112 + 21 + h
13 + 3h
h3 + 2h +
S h2 + 3h
7 37
$ +
amm"
where the first line plugs in 7 = n2; the third line drops the h 2h in the denominator; the forth lines plugs 72; the fifth line is because 72 dominates the bound when m large. Thus, we have shown = (dm + m) = (addm
The above lemma shows that forward differencing designs achieve the same asymptotic decay rate as the central differencing designs for the RBF kernel. Though, the leading constant the big O notation is slightly larger.
C Convergence Proofs
The following is a useful lemma for biased gradient updates, which bounds gradient norm via the cumulative bias. The proof is adapted from lemma by Ajallocian and Stich [1. Lemma 2].
Lemma 15. Let J be L-smooth and bounded from below. Suppose the gradient orucle gs has bias bounded by & in the 1-th iteration. Namely, we have l/g. - gil! 11 S & for all 0, where g: fxx) is the ground truth gradient. Then the update Xx+1
7gg with n 7 produces
a sequence ( satisfying
2(f(x)) 12x
(-)
IT 1= E Σr= 1)t
min ISST
+
Σt= 11
21

---

Proof. By L-smoothness, we have
f(x(+1) s + (x:+1 x,) + 2 x x111?.
Plugging in the update formula X(+1 = X; 7).gr, we have
f(x(+1)= - f(x)) V f(xx)) T 82 + LiT |
$ f(x)) fx())) T &1 + T/ 1l
2
$ f(x) TN vf(x), T 8: + 2 vf((x,) + V
fx1) + (0 v [Vf(x;))
f(x)) [ + 7cE..
where the first inequality uses L-smoothness; the second inequality uses n) S 1, the fourth inequality expands the squared Euclidean norm: the last inequality uses the definition of bias. Summing the inequalities for t 2, I'' and rearranging the terms, we have
Th [V ((((x))) ((XT+1)) + Σ ne
T
$ 2(([(x)) - -1)) + Σ nche.
Dividing both sides by Σ ?) gives
Ee me f((x))1 2(f(x)) -r)
Σt= m+Et
S
El
Et 1/t
E
The left hand side is a weighted average, which is greater than the minimum over 1 S t S T, which completes the proof.
Next, we prove a variant of Lemma 15 for the projected update x (x. - - fxx,)) For the ground truth gradient g: fx.). define the gradient mapping
G(x:) (x:
projx 1//g:))
T/t
For the approximate gradient g: = VpD, (x1), define the gradient mapping
(xx,) (xi - projox T)g:
T)
In the following, we introduce two lemmas characterizing the projection operator pr) Lemma 16 (e.g., Lemma 3.1 of Bubeck et al. [2]). Let be convex and compact. Let x E X and z E RR Then we have
(projx(z) T (x proix(z)) 2 0.
As result,
|x- 2|22 [projx(z) 2112 ||x (z)1x
Lemma 17. The following holds:
1. [G(x;)] l|ga
2. [G(xx)) - gll S |gg]]
3. [[(xx)) G(x,)) [I
l|g. -g.l).
22

---

4. g' gg G'xx)) 2 |
Proof. (1-2): The first two inequalities are direct corollaries of Lemma 16.
(3): This is proved as follows:
[lū(x.) G(x))l = Irrojx (x: Thg:)) - propoxxxxxx 17:8:)
7/t
g Bell
where the second line is because the projection operator is non-expansive.
(4): By Lemma 16, we have
(xx) (n:)((x)
(T) - -Tgg:)) 2 O.
Rearranging the terms finishes the proof.
Now we give lemma proving biased gradient update with the projection operator. The proof is adapted from a lemma by Shu et al. [22].
Lemma 18. Let be L-smooth over convex compact set X. Moreover, assume the gradient norm fx)]l is bounded by L' on X. Suppose the gradient oracle &t has bias bounded by & the t-th iteration: ||g: gel? for all 2 O where gi fx(x))- Then the update Xi+1 = (x meSe) with m7 S LLP produces a sequence {x, satisfying
[Gxx) ((f(x))-ff +
ET ThE L'ET nEVE
min 151ST
Σtm A1
Etm The
Eim Me
where G(x,) 1 dt
pxojx (xooixx m(Bi)) is the gradient mapping.
Proof. By L-smoothness, have
f(x(+1) f(xi) S [x - x1) + L |||x++1 - x111
1
ng ((x,) + LnZ |(((x))11
2
1
- G(xx) + 2 Thl ||G((x)))1?
where the second line plugs in the update X (xx) and the third line is due to m/ S +. Now we analyze the two terms separately. For the first term, we have
-eB. T G(x,)
n. T (ā(x,) G(x,)) ngg. G(x,) S <-n.g, T (ā(x.) G(xx)) - n/ (xx)|
where the second inequality uses Lemma 17. For the second term, we have
1 2 nh G(x))]
n ((()) G(x²) G(xx)]l
nh
n (G(x)) - G(Xi)) G(x,) + 1 (G[())))))
no [|(((x)) G(xx))l + n
< 2 Th m:
1l.. - B:112 + ne (G(x)) G(xi)) G(xx) + [nl[G((x))]
Summing them together, we have
1
gi112 2 n (G(x)) G(x:)) T (G(x.) - gr) 1 31
f(x+11) - f(xe)= 5
Th 18:
5
1
gil22 n ((x:) - G(xx)) ||g:11 7 G(x)))1
2
1
1
2 nc5 + L' L'vEe=
2 G(
23

---

where the second line uses Cauchy-Schwarc inequality.
telescoping gives
T
T
T
Σ
The GG(xx)I &
Σ Σ
nhE + 2L
TKe
-
t=1
which results in
2(S(x1) -n)
Σ
ii= Mest
T
L'E
noEE
min ISTT
|G(x)
Ei=1 11
E 1/1
Σi
The rest of this section proves all theorems and their corollaries in the main paper.
Theorem 1. Lei f E H whose smoothness constant is L. Running Algorithm with constant batch size be b and step size 1J. I for T iterations outputs sequence satisfying
min
B² ff(x)) > (2L(((x)) ff)) - B². Ed,,,o(b). Eb)
(5)
Proof. By Lemma 1 and Assumption 3, we can bound the bias in the iteration & as
[Vx(x V.D, (xx) I B LrkV, (Xe, V
By Lemma 8, the trace in the RHS can be bounded by the error function Ed (b). Thus, the gradient
bias is B Ek (b). Applying Lemma 15 with 7) T. and & = B² E (b) finishes the proof. O
Corollary 1. Under the same assumptions of Theorem 1, using batch size b, = d + 1, we have
min [V =
Therefore, the total number of samples n = OdT') and the squared gradient norm [Vx(())] converges to zero at the rate (d/n)..
Proof. By Lemma 2, we have Ea (d + 1) 0. Plugging it into Theorem 1 gives the rate in the
iteration number T. To get the rate samples n,, note that n Σin (d + 1) = (d + 1)T. Plugging T Z+I into the rate finishes the proof.
Theorem For
suppose a GP sample whose smoothness constant L w.p. least 1 - 8. Algorithm 1 with batch size by and step size T)t + produces sequence satisfying
[xxx)))11 < (2L((((x)) (*)) + E C.Ea,, ISIT
min
(6)
with probability al least 1 26, where C. = ( (=22
Proof. By Lemma 11, we have
[V)( (xD, < C. C,tr (x,,
with probability at least 1 5. The trace on the RHS can be further bounded by the error function E (bi) by Lemma 8. Applying the union bound, with probability at least - 28, the inequality
[VV Vpp, (x:)) < C.Fa,k,o
holds for all 2 and f is T-smooth. Applying Lemma 15 with n = and g C2 C.Fa,,, (b:) finishes the proof.
24

---

Corollary 2. Let k(.--) be either the RBF kernel or the V/ = 2.5 Matérn kernel. Under the same conditions as Theorem if
d doog t;
O(/TT) + +
b.
dt;
then
min 151ST
(odT log T) O log n)): O log" T) O n log n)
de2
with probability at least 1 1 28. Here, the total number of iterations and n the total number of samples queried.
Proof. By Theorem 2, with probability at least - 26, we have
min /V(xx S (21(f(x)) f*)) 7 ET= C.E.kk,o
The proof boils down bounding the average cumulative bias. The full details in Appendix D.
Theorem 3. Under the same conditions as Corollary 2, without Assumption 2, using the projected
update rule (7), Algorithm 1 obtains the following rales:
{ dt; di2,
O(odin logn + +oid logn); Olodin log" + o log n),
if b.
then
min 1S1ST
[[(x)))
with probability at least 26. Here, n is the total number of samples queried.
Proof. Since is twice differentiable on compact set X. its gradient norm [x)]] attains maximum. Thus, there exists a constant L' such that 2 |V((x)] for all x E X. By Lemma 18 and similar argument in Theorem 2, with probability at least 1 26, we have
min [[(xx))1 < (2L(f(x))) - ft) + + Cf2 (b)) + Et= 0.2)) (b.), ISST
where cl and C constants growing in (logt). t). By Lemmas 4 and 5, plug in the error function k (b) (ad2bb The rest of the proof follows similar argument in Corollary as shown in Appendix D.
Finally, we present convergence result for GP sample path under noiseless assumption.
Theorem 4. For 0 < 1, suppose f GP sample whose smoothness constant is with probability at least 1 - 8. Assuming O, Algorithm I with batch size bt = d + 1 and step size T] = 2 produces a sequence satisfying
2L(f(x1) - [V
min SSST
T
with probability at least 1 28.
Proof. By Theorem 2, we have
min fxx)) - SEST
T
+
Σ
C.Eak (d
T
T
with probability least 25. By Lemma 2, F.a (d + 1) O. the second term essentially zero, which finishes the proof.
25

---

D Optimizing the Batch Size
In this section, we optimize the batch size in Theorem and give explicit convergence rates. The discussion in this section will give a proof for Corollary 2.
For the RBF kernel and 2.5 Matérn kemel, by Theorem 2, Lemma 4 and Lemma 5, we have shown the following bound
min ISrT
[V,,)
- T
Σ
Ed or (bx) logt
T
t=1
od
T
-
bt
log t.
T
We discuss polynomially growing batch size by = dla, where > O. Then we have
T
min ISTST
[Vx(x
2 7 + od
da
log t.
T
t=1
We discuss three cases: 0 < a < 2, a = 2 and a > 2.
Case When S 2, the infinite sum Σ ta logt diverges. Its growth speed is on
the order of O(T 1a log T'.. The total number of samples = Ei= b. Oddrat1), and thus Od मों n Thus, the rate is
min ISST
Vf(x) < T od T logT
3ul2 d क ad
n 2αाh log n
Jo+2 S od 20+1)
logn,
where the last inequality the faet that the second term dominates the when 0 a
Case 2. When = 2, the infinite sum Σ 1= t- logt diverges. Its growth speed is the order of
Ologg2 T). The total number of samples is n = d' t = -(dTT), and thus = Old 3n3). Then the rate is
1
ad log T
min ISST
1
T
T
x adT log log21 T
S 0d - log
n.
Case 3. When a > 2, the infinite sum Ei= t ja logt converges to a constant when T - cc. The total number of samples n (dd'*+] Thus, the rate is
1=
min Vx(x) T odd
ad n +|.
Note that a = 2 achieves the fastest rate O(n in terms of samples n.
Now we discuss batch size with logarithmic growth b, dloog t. We have
T
min
S T
+ od
Σ
(1)
T 1=1
x + ad.
26

---

E Additional Experiments
This section presents additional experimental details and additional numerical simulations.
Details of Figure 2. In Figure 2a, we plot the error function starting from b = 20 to make sure b 2 2d so that Lemma 5 indeed applies. The decay rate ((dd b b- has (leading) hidden constant of 15y2 inside the big O notation (see the proof of Lemma 5), and thus the bounds plotted in Figure are multiplied by this constant. Otherwise, the expression ad b alone is not valid upper bound.
2
The ReLU function max 3 is non-differentiable at x convexity the subdifferential at = 0 is defined as [0, 1]. We estimate the "derivative" at = by minimizing the acquisition function. The estimate (0) and queries are plotted in Figure 5. The estimated derivative (0) is always in [0, 1]. Thus, the posterior mean gradient y (0) produces subgradient this case.
ReLU.
Nevertheless, thanks to
1.0
ReLU
1.0
ReLU
data
0.8
data
x
0.6
0.4
0.4
0.2
0.2
0.0
0.0
-2.0 -1.5 -1.0 -0.5
0.0
0.5
1.0
-2.0 -1.5
-1.0
-0.5
0.0
0.5
1.0
(a) n - 2, (0( = 2.06 x
10-5
(b) n 4, vb(0) = 0.49
Figure 5: Estimating the "derivative" of ReLU at 0 with noisy observations (o 0.01).
27


Reviews:
[{"review_id": "zrfG6ZqEzA", "reviewer": "Reviewer_aL9q", "summary": "This paper gives theoretical considerations on local Bayesian optimization methods for black box optimization. In particular, for gradient-based local BO, when the objective function satisfies the smoothness assumption, the rate of convergence to the stationary point is derived for the noiseless and noisy cases, respectively. \nIn the noiseless case, local solutions can be found at a much faster rate with respect to the input dimension d than with conventional global BO. \nAlthough the convergence rate is worse in the presence of observation noise than in the noiseless case, it returns to the noiseless rate in the limit of zero noise variance, suggesting that it is a reasonable bound.\nThe tightness of the derived convergence rates is validated by numerical experiments.\n", "questions": "- I think one advantage of local search is that it is compatible with parallel distributed processing (such as TurBO). Is it possible to extend the gradient-based local BO treated in this paper to parallel distributed processing as well? And if so, would some theory of  convergence for such methods be obtained by a direct extension of the results derived in this paper? \n\n- Theorem 3 derives the upper bound of the error function for the squared exponential kernel. Can we easily obtain the upper bound of the error function for other kernels (e.g. Matern kernel)?\n\n- In Lemma 2 and line 226, is $E_{d, k, \\sigma}(b_t, d)$ mistake for $E_{d, k, \\sigma}(b_t)$?", "limitations": "The authors appear to have adequately addressed the limitations of the study and, where applicable, the potential adverse social consequences.", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The rate of convergence of the gradient-based local BO algorithm to the local optimal solution is theoretically derived, and the rate is expected to be somewhat tight (empirically validated).\n\n- It is clear how the observational noise affects convergence in the local BO. Furthermore, by taking the limit of the noise variance parameter, the convergence rate connects the noiseless and noisy cases.\n\n- The convergence rate results for the noiseless case are given for both cases where the objective function is an element of the RKHS (in the main text) and a sample path of a Gaussian process (in the appendix).\n", "weaknesses": "- The local BO discussed here is a gradient method type algorithm called GIBO in particular, and approaches that construct the GP itself locally, such as TurBO [Eriksson+ 2019], have not been considered. \n\n- It would be nice to have a theory that connects convergence to a stationary point and regret analysis. In Section 3, the results suggest that local solutions are good enough compared to the global solution, and if this is a property that holds to some extent universally, then it seems reasonable to expect that the regrets of local search also have good properties (e.g., sublinear properties).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments. We provide our response as follows.\n\n**Comment:** It would be nice to have a theory that connects convergence to a stationary point and regret analysis.\n\nWe agree this is an interesting question to convert our convergence rates to global regret bounds. However, this is actually a fairly challenging problem that would probably constitute a large technical lift in its own right. We’d need to analyze the regret of finding a sequence of stationary points using random restarts. Essentially, this is a question about the relationship between stationary points of GP sample paths and the global optimum. To the best of our knowledge, there is little research on this question, although some of the theories by Adler (1981) may be relevant. We think this is interesting future work.\n\n**Comment:** Is it possible to extend the gradient-based local BO treated in this paper to parallel distributed processing as well?\n\nFirst we note that, while it’s definitely not obvious, search direction based BO methods actually do effectively use batch evaluation already! Most of the evaluation budget of methods like GIBO and MPD is spent acquiring batches of data Z to learn about the gradient. Optimizing the acquisition function for Z can and often is done for relatively large (tens to hundreds) batches of points at a time. Thus, even the sequential single restart version of GIBO makes reasonably effective use of parallel evaluation capabilities. If additional parallelization is required, considering multiple restarts simultaneously is possible. \n\n**Comment:** Upper bound of the error function for other kernels.\n\nWe do now have an additional result on the Matern kernel ($\\nu = 2.5$). Interestingly, we are able to prove the Matern kernel has the same rate as the RBF kernel, i.e., $E_{d, k, \\sigma}(b) = \\mathcal{O}(\\sigma d^{\\frac32}b^{-\\frac12})$. This new proof slightly adapts the techniques in the paper. The main new technique is to use a rational function upper bound on the posterior covariance trace for the Matern kernel, which allows us to analyze the error function by minimizing this rational function upper bound.\n\nWe conjecture that our proof can be generalized to the **entire** Matern family (widely used in practice) and the rational quadratic kernel. Additionally, we remark that the error function we proposed in this paper is a quantity analogous to the information gain in global BO analysis — bounds on the error function for **any** kernels immediately result in convergence rates by our proof.\n\n**Comment:** In Lemma 2 and Line 226…\n\nCorrect. Indeed, the dependency on the dimension is through the subscript, not the argument. Thanks for spotting this typo.\n"}, {"review_id": "9nyMppLVEI", "reviewer": "Reviewer_kgRX", "summary": "The paper investigates the behavior and convergence property of the local Bayesian optimization approach, in particular, the method GIBO [16]. The paper provides the convergence rates of GIBO in both the noiseless and noisy settings. The paper also performs various empirical experiments to understand the tightness of the derived bounds. It also conducts some experiments to understand some other properties of local BO like multiple restarts and when the objective function is non-differentiable.", "questions": "Please answer my comments in the Weaknesses section, and the following additional questions:\n+ In the experiments (Section 6.1), how are the hyperparameters of the bounds or the algorithms set (like \\beta) ?\n", "limitations": "The paper seems to not have a dedicated section describing the limitations of the work presented in the paper.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "+ The paper tackles an important and difficult problem, which is to perform theoretical analysis for the local BO approach\n+ The paper writing is very clear and easy to understand. The paper is really technical; it has a lot of theory, but thanks to the clear writing, all the theoretical analysis is much easier to understand.\n+ The theoretical analysis on the convergence rates of GIBO seems to be sound to me (based on my overall understanding regarding the theoretical analysis of BO – I couldn’t dig into the proof in a very detailed manner).\n+ Some empirical experiments are also conducted to understand the tightness of the derived bounds.\n", "weaknesses": "+ The paper actually focuses only on one of the local BO approaches (the GIBO method) and this approach is actually quite different compared to some other popular local BO approaches like TurBO. The GIBO approach is based on gradient descent, but the other local BO approach like TurBO is still based on BO. The paper should make very clear regarding the focus of the theoretical analysis.\n+ Another weakness of the paper is that it still relies on the assumption that are used to analyze the global BO including the assumption that the objective function is a sample from a GP with known hyperparameters. Actually one of the main motivation of the local BO approach is based on the observation that the objective function can’t be model using one single global GP in the whole search domain. Therefore, the assumption of the theoretical analysis in this work seems to be quite strict, and therefore, the theoretical analysis (the bounds) may not shed too much light on the performance of the local BO approach.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments. We provide our response as follows.\n\n**Comment:** The paper actually focuses only on one of the local BO approaches (the GIBO method).\n\nSee the global response above.\n\n**Comment:** Sample path assumption with known hyperparameter.\n\nWhile this assumption is standard and widely used in the theoretical BO literature, there is also potential hope for relaxed assumptions for local BO. For local BO, we don’t need the GP to be a “potentially reasonable” surrogate model of the function everywhere in the input space. Intuitively, we only need sufficient smoothness to be able to estimate the gradient at the current iterate $x_t$ – all other areas of the input space are irrelevant at iteration t. It’s interesting future work to consider what the weakest assumptions are that enable provable GP-based gradient estimation. We plan to explore this possibility further in future work.\n\n\nNevertheless, we believe that our results already reveal some advantages of local BO algorithms. Our convergence rates (for local stationary points) are **polynomial** both in the dimension $d$ and the number of samples $n$. This is in sharp contrast to any global BO algorithm, e.g. GP-UCB. In this sense, our results provide partial justifications for local BO algorithms and their fast convergence in practice (at least for a particular algorithm, GIBO).\n\n**Comment:** In the experiments (Section 6.1), how are the hyperparameters of the bounds or the algorithms set (like \\beta) ?\n\nFigure 2 studies the tightness of the bound on the error function (Definition 2), which depends only on the kernel and not the objective function, and therefore makes no use of $\\beta$. The kernel used in Figure 2 is a RBF kernel with a unit length scale, as in Theorem 3. The purpose of Figure 2 is to compare how quickly the posterior derivative variance can be decreased empirically to our theoretical upper bound on that rate. Since the posterior variance doesn’t depend on the actual labels collected (with fixed hyperparameters), and since the error function is measured at the origin with only the collected batch of data, the results in Figure 2 are independent of collected labels or any particular objective function.\n"}, {"review_id": "QndLhLIER8", "reviewer": "Reviewer_GVja", "summary": "Local Bayesian optimization methods are widely used to cope with the curse of dimensionality when optimizing high-dimensional black-box functions. Although these methods have very good performance empirically, virtually all theoretical results focus on the global optimization setting by proving bounds on the optimality gap. This paper shows that under typical assumptions a particular variant of local Bayesian optimization (GIBO) provably converges to a stationary point in polynomial time. ", "questions": "- based on the description in line 113 it seems the x axes in Figure 1 are mislabeled, is that correct?\n\n- I found the following statement on line 220 a bit vague: \"However, because observations of f are made with noise, the noise inherit too [sic] both evaluations of f is amplified by 1/2h as h -> 0.\" Can you please elaborate?\n\n- I found the second paragraph of the discussion intriguing. Do you think the noise in the BayesOpt gradient estimate may play a similar role as the minibatch gradient noise in SGD? My understanding of the argument is minibatch gradient noise tends to cause SGD to find wider basins in the loss surface, which some have argued is linked to better generalization. In this context we aren't concerned with generalization per se, but perhaps one could argue that wider basins of f are preferable since they will be more forgiving if the estimate of the stationary point is slightly off. \n\n- in practice it is quite common to estimate the GP surrogate hyperparameters through the log marginal likelihood, can you comment on how this might affect your analysis and empirical behavior?\n\n- how might one start thinking about \"convergence\" in the multi-objective case? Could we prove that local BayesOpt converges to a locally non-dominated point?\n\n- although this paper did not propose GIBO, can you comment on its numerical stability in practice? In my experience look-ahead acquisition functions tend to exhibit poor numerical behavior, which I imagine is only amplified by introducing an additional derivative operator.", "limitations": "I see two primary limitations\n\n- the analysis is restricted to GIBO. It would be interesting to know if any similar arguments can be made for TurBO [1], which is quite widely used. \n\n- the analysis relies on fairly strong (albeit common) assumptions on the nature of f and the specification of the surrogate. It seems plausible to me that (for example) BayesOpt could find a stationary point even if the surrogate is misspecified, assuming f is not too pathological.\n\n[1] Eriksson, David, et al. \"Scalable global optimization via local Bayesian optimization.\" Advances in neural information processing systems 32 (2019).\n", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "I greatly enjoyed this paper. It is well-written and easy to follow, and answers important questions with proof techniques that are novel to me. I am very happy to see that the authors included some empirical studies to verify the integrity of their results. ", "weaknesses": "This work is ready for publication.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments. We provide our response as follows.\n\n**Comment:** The x-axis of Figure 1.\n\nThanks for spotting this typo! Indeed, the x-axis in Figure 1 is the indices of the dimension list {1, 5, 10, 20, 30, 50}. We will replot the x-axis of the figure.\n\n**Comment:** Line 220.\n\nWe will rewrite this sentence in the revision. The central differencing scheme has $O(h^2)$ truncation error, and thus $\\frac{\\partial f}{\\partial x} = \\frac{1}{2h} \\big( f(x + h) - f(x - h)\\big) + O(h^2)$. To reduce the estimation error, we need a small $h$. However, the function evaluation is corrupted by a Gaussian in the presence of noise. The right hand side becomes $\\frac{1}{2h} \\big( f(x + h) - f(x - h) + 2 \\epsilon\\big) + O(h^2)$, where $\\epsilon$ is Gaussian ($2\\epsilon$ is due to two noisy function evaluations). When $h \\to 0$, the denominator $\\frac{1}{2h}$ blows up the noise. This illustrates that noise creates a sort of general hardness in estimating gradients, whether you’re using a GP or not.\n\n**Comment:** Second paragraph of discussion. Do you think the noise in the BayesOpt gradient estimate may play a similar role as the mini batch gradient noise in SGD?\n\nThat’s an interesting question. We think the answer might be yes potentially. Besides the case you mentioned, there is another intuitive case. Due to the estimation error, local BO may be less likely to be trapped by saddle points (a particular type of bad stationary points). The intuition is that convergence to saddle points is unstable — a slight perturbation may let the algorithm find a decent direction and therefore escape the saddle point.\n\nHere we use the term *algorithmic bias* to broadly refer to this phenomenon: local BO algorithms (e.g. GIBO), are able to find good solutions (which is demonstrated by their strong empirical performance in previous research), despite the existence of potential bad local solutions/stationary points. It is not completely well-defined unfortunately but may be worth studying in the future.\n\n**Comment:** Estimating the hyperparameters through maximizing log likelihood.\n\nOur theoretical analysis assumes a kernel with fixed hyperparameters, which is relatively standard in the theoretical BO literature. Extending to unknown hyperparameters is an interesting future work. In practical settings, we should of course generally expect that estimating the hyperparameters improves the model fit, which typically in turn improves the convergence. The empirical performance of local BO methods in settings where hyperparameters are learned is of course well established in prior work.\n\n**Comment:** Convergence in the multi-objective case.\n\nIt might be possible to extend our results in multi-objective settings. One of the intermediate steps in our proof is to show that the objective decreases in every iteration if the gradient estimation error is small (see the equations below Line 395). If we can show a similar property holds in the multi-objective setting (e.g., each iteration decreases all objectives), then it might be possible to prove convergence. Note that decreasing objectives locally is a very natural requirement for local BO algorithms. We think it is interesting to check if certain multi-objective BO algorithms satisfy this property (or propose new algorithms satisfying this property).\n\n**Comment:** Numerical stability of GIBO.\n\nOptimizing the acquisition function of GIBO is generally stable with good initialization. The resulting designs typically center around the current iterate $\\mathbf{x}$ (e.g. Figure 4). In practice, however, updating the iterates $\\mathbf{x}_{t + 1} = \\mathbf{x}_t - \\eta \\mathbf{g}_t$ may require careful attention. The estimated gradient $\\mathbf{g}_t$ is only an **approximation** of the true gradient. While in theory the approximation error can be bounded by our theorems, in practice the approximation error is affected by two additional complications: a) model misspecification and b) small batch size in the noisy setting. Thus, it is recommended to use small step sizes to avoid overshooting, or even line search strategies. Oftentimes normalizing the estimated gradient $\\mathbf{g}_t$ is helpful (Muller et al., 2021).\n"}, {"review_id": "lRsb6E4Z0P", "reviewer": "Reviewer_pxRe", "summary": "Local Bayesian optimization has been popular high-dimensional optimization methods, which achieves the state-of-the-art empirical performance over various benchmarks. This paper investigates the behavior of local Bayesian optimization (BO) methods. Under a prototype algorithm, the authors empirically  demonstrate superior performance of the stationary points found by local BO compared to vanilla BO and grid search. Under the assumption that the function generating parameter are known, they theoretically and experimentally analyze the convergence rate of local BO methods. ", "questions": "1. Many experiment in this paper use GP sample path as the objective function, with a dimension up to 100 (Figure 3). How to sample a function from GP under such high-dimensional input space? Any discretization is used during the function sampling?\n", "limitations": "The authors have addressed the limilations. I did not find potential negative societal impact of this paper.", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. I think the study problem in this paper is important and the theoretical experiment results partially explains the good performance of local BO.\n\n2. The theoretical analysis is easy to follow, and gives insight of how to design local BO hyperparameters according to the problem dimension.", "weaknesses": "1. The prototype local BO algorithm samples points to reduce the gradient estimation unceratinty, and move search center according to the estimated mean gradient, which may not be representative of all local BO methods such as trust region Bayesian optimization (TuRBO) [1].\n\n2. The assumptions in theoretical analysis may not hold in many real-world problems.\n\n[1] Eriksson, David, et al. \"Scalable global optimization via local Bayesian optimization.\" Advances in neural information processing systems 32 (2019).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments. We provide our response as follows.\n\n**Comment:** The assumptions in theoretical analysis may not hold in many real-world problems.\n\nWe definitely agree that these assumptions may not hold in real-world problems – most real world functions are not GP sample paths. The assumptions we use (RKHS/GP sample) are the most commonly used assumptions in the BO theory literature because proving convergence in the black-box setting with arbitrary model misspecification is generally a pretty daunting task.\n\nWith that said, interestingly, there actually is hope for relaxed assumptions for local BO. For local BO, we don’t need the GP to be a “potentially reasonable” surrogate model of the function everywhere in the input space. Intuitively, we only need sufficient smoothness to be able to estimate the gradient at the current iterate $x_t$ – all other areas of the input space are irrelevant at iteration t. It’s interesting future work to consider what the weakest assumptions are that enable provable GP-based gradient estimation.\n\n**Comment:** How to sample a function from GP under such a high-dimensional input space? Any discretization is used during the function sampling?\n\nWe did not use discretization in our experiment -– see Line 99. We adapted the pathwise sampling technique of Wilson et al. (2021) as implemented in BoTorch, but using exact conditional calculations rather than variational / inducing point conditional calculations (since our goal was not improved asymptotic complexity).\n"}]

Author rebuttal:
We thank all reviewers for their encouraging comments and helpful feedback on our submission.

All reviewers share a common question of whether and how our theory might be extended beyond descent direction approaches like GIBO to trust region algorithms like TuRBO. In general we agree it would be very interesting to analyze TuRBO theoretically, and we’re interested in investigating this as well. While at this point we think it’s probably possible, the analysis is slightly more involved than GIBO. In particular, we believe TuRBO must be more significantly modified than GIBO in order to analyze it. In particular, the rules for shrinking and expanding the trust region will likely be different than simply “grow after k successes, shrink after k failures.”

We will revise relevant sentences in the paper to make it clear that we have only proved the convergence rate for search direction style local BO methods and list the convergence of TuRBO as important future work.


Meta-review: This paper investigates a theoretical behavior of the local Bayesian optimization, presenting its convergence rate for both noiseless and noisy cases. All reviewers agree that the paper has a solid contribution, touching an important and timely subject in BO. The authors did a good job in the author response. I believe that the paper is ready for publication.
Decision: Accept (spotlight)