This Too Shall Pass: Removing Stale Observations in Dynamic Bayesian Optimization
Authors: Anthony Bardou
Patrick Thiran
Giovanni Ranieri
Venue: NeurIPS 2024 poster (2024)
Keywords: Dynamic Bayesian Optimization
Black-Box Optimization
Zeroth-Order Optimization
Primary area: probabilistic_methods

Abstract: Bayesian Optimization (BO) has proven to be very successful at optimizing a static, noisy, costly-to-evaluate black-box function $f : \mathcal{S} \to \mathbb{R}$. However, optimizing a black-box which is also a function of time (*i.e.*, a *dynamic* function) $f : \mathcal{S} \times \mathcal{T} \to \mathbb{R}$ remains a challenge, since a dynamic Bayesian Optimization (DBO) algorithm has to keep track of the optimum over time. This changes the nature of the optimization problem in at least three aspects: (i) querying an arbitrary point in $\mathcal{S} \times \mathcal{T}$ is impossible, (ii) past observations become less and less relevant for keeping track of the optimum as time goes by and (iii) the DBO algorithm must have a high sampling frequency so it can collect enough relevant observations to keep track of the optimum through time. In this paper, we design a Wasserstein distance-based criterion able to quantify the relevancy of an observation with respect to future predictions. Then, we leverage this criterion to build W-DBO, a DBO algorithm able to remove irrelevant observations from its dataset on the fly, thus maintaining simultaneously a good predictive performance and a high sampling frequency, even in continuous-time optimization tasks with unknown horizon. Numerical experiments establish the superiority of W-DBO, which outperforms state-of-the-art methods by a comfortable margin.

This Too Shall Pass: Removing Stale Observations in Dynamic Bayesian Optimization
Anthony Bardou IC, EPFL Lausanne, Switzerland anthony. bardou@epfl.ch
Patrick Thiran IC, EPFL
Giovanni Ranieri IC, EPFL
Lausanne, Switzerland patrick. patrick.thirangepie. thiran@epfl.ch
Lausanne, Switzerland giovanni.. iewranni.raiiee er ereerriieee @epfl. . ch
Abstract
Bayesian Optimization (BO) has proven to be very successful at optimizing static, noisy, costly-to-evaluate black-box function S E. However, optimizing a black-box which is also a function of time (i.e., a dynamic function) f : x - R remains challenge, since dynamic Bayesian Optimization (DBO) algorithm has to keep track of the optimum over time. This changes the nature of the optimization problem in at least three aspects: (i) querying an arbitrary point in S T is impossible, (ii) past observations become less and less relevant for keeping track of the optimum as time goes by and (iii) the DBO algorithm must have a high sampling frequency it can collect enough relevant observations to keep track of the optimum through time. In this paper, we design a Wasserstein distance- based criterion able to quantify the relevancy of an observation with respect to future predictions. Then, we leverage this criterion to build W-DBO, a DBO algorithm able to remove irrelevant observations from its dataset the fly, thus maintaining simultaneously a good predictive performance and high sampling frequency, even in continuous-tin time optimization tasks with unknown horizon. Numerical experiments establish the superiority of W-DBO, which outperforms state-of-the-an methods by comfortable margin.
Introduction
Many real-world problems require the optimization of a costly-to-evaluute objective function f Rd + with an unknown closed form (i.e., either the closed form expression of f exists but remains unknown to the user, does not exist). Such setting occurs frequently, and examples be found in hyperparameters tuning [ networking 2 3 robotics 4 or computational biology [5]. In such applications, f can be Bayesian Optimization (BO) is an effective framework for black-box optimization. Its core idea is to leverage surrogate model, usually a Gaussian Process (GP), query specific inputs. By doing so, BO algorithm is able to simultaneously discover and optimize the objective function.
black-box and cannot be optimized by usual first-order approaches.
Since its inception, BO has proven be very effective at optimizing black-boxes in variety of contexts, such as high-dimensional input spaces 6 7 g]. batch mode [9] or multi-objective optimization [1O]. However, few works study BO in dynamic contexts (i.e., with time- varying objective function), despite its critical importance. Indeed, dynamic black-box optimization problems arise whenever optimization task is conducted within environment that comprises exogenous factors that may vary with time and significantly impact the objective function. Dynamic black-boxes found in network management (T1). unmanned aerial vehicles tasks [2]. hyperparameter tuning in online deep learning [13] online clustering [14] or crossing waypoints location in air routes [15]
In dynamic context, ; S T - IR is tame-varying, black-box, costly-to-evaluate objective function with spatial domain C Rd and temporal domain T C IR. Unlike common black-box
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

objective functions that take as imput a point a E S in a spatial domain only, the function f takes as input point (x. t) É S x T in space and time. Since the problem of optimizing iis still addressed under the BO framework, the framework is called dynamic Bayesian optimization (DBO).
Taking time into account does not boil down to merely adding extra dimension to S. It changes the nature of the optimization problem on at least three aspects: (i) at present time to. a DBO algorithm can query any arbitrary point æ S but cannot query past (i.e., points in time, (ii) as time (only) moves forward, a previously collected observation (12 t), f(x, becomes less and less informative about the future values of f as time goes by and (iii) the response time (i.e., the time required for hyperparameters inference and acquisition function maximization) becomes a key feature of the DBO algorithm since has a direct impact on the sampling frequency of the algorithm and, consequently, on its ability to track the position of the optimum as f changes.
tq) nor future (i.e. to)
Interestingly, (ii) and (iii) imply that each observation eventually becomes stale and, such, computational burden if kept in the dataset. Since a DBO task might require the optimization of an objective function over arbitrarily long periods of time, the ability to remove observations from the dataset on the fly, as soon as they become irrelevant for future predictions, is essential to prevent a prohibitive growth of the response time of DBO algorithms. Overall, (i), (ii) and (iii) require to address dynamic (i.e., space-time) problems differently from usual (i.e. space-only) problems.
The main contributions of this article twofold. First, we propose fast and efficient method able to quantify the relevancy of an observation. Second, we leverage this method to build a DBO algorithm able to identify and delete irrelevant observations in online fashion
2 Background
Let us start by describing the BO framework with a GP prior, as introduced by [[6]. Given an objective function f S C Ed R. BO assumes that f is a GP (0, k(x, x')). For any æ E S, and given a dataset of observations D = {(z,, 8:) Be))ien..n where x; E S is a previously queried input with (noisy) function value yı- [ax)) + e N (0, the posterior distribution of f(x) is N(((()), (z)), where
m(x)) k. (xx, X) A y
(1)
(²()x) k(2,, x) k (ar, x)A--k((s X) with X = (x), xn). y = (V1,-.. .Vn), A = k(X, X) + 0 and k(X,)) (k (2) x3)) 2,, (X,"E)))
(2)
=
To find In+1, the input to query at the (n + 1)th iteration, a BO algorithm exploits an acquisition function S + R. The acquisition function o quantifies the benefits of querying the input æ in terms of (i) exploration (i.e., how much it improves the GP regression of f) and (ii) exploitation
(i.e., how close it is to the optimum of f according to the GP). A large variety of BO acquisition functions have been proposed, such as GP-UCB [17, Expected Improvement [8] or Probability of Improvement
19 Formally, the BO algorithm determines its next queried input by finding
2n+11 = arg max (xx)
BO extends naturally to dynamic problems. by adapting the covariance function A to properly capture temporal correlations (discussed later in this section). The resulting inference formulas are very similar to (I) and (2). Given a dataset of observations D = {((x,, fz) (.),)))))) be))ee where (ax; t,) S T is previously queried input with the (noisy) function value y; f(x,, t,) + c, the posterior distribution of [(x, 1) is N (u(æ, l), ax, ()) for any (æ, () E S x T, with
μ( t) = k. T
((x, t). -)-- -yy,
(3)
0² t) = k((xx) t), (a,t)) t)) k ((x,t) ()) *k(xxt) (x, X) ,yn), with X = ((11, (In, tn),,, y (y1, ,.))), A k(X, X) + I and k(X
(4)
=
(Ie ((x,,tt) (z),,t))) (.)) (
At least, not immediately.
Iss Python documented implementation can be found at httpu//////////gttttooo/o/o///// Httm///DDDOALLORTTH///AU
algo A PyPI package can be quickly installed with the command pip install adbo-algo.

---

Exploiting the usual acquisition functions in dynamic context is also straightforward. Since DBO algorithm can only query an input at the current running time to. the next queried input is simply (In+1. fo). with xn+1 = arg m sp(æ, lo). Some DBO algorithms (e.g., (20]) extend the querying horizon to the near future, that is, from in to a time interval (Lo, to + dz]. In that case, the next queried input is (Im+11 f1) = arg (x. ((xx, t).
BO is an active field of research, but relatively few works address DBO, despite the natural extension of BO dynamic problems described above. We conclude this section by reviewing them. In (2], the objective function is allowed to evolve in time according to a simple Markov model, controlled by hyperparameter e E [0, 1]. On the one hand, the authors propose R-GP-UCB, which handles data staleness by resetting the dataset every ((e) iterations. On the other hand, the authors also propose TV-GP-UCB that incorporates data staleness by weighing the covariance of two queries (1-)) c) In [22], the authors use the same model with an event-triggered reset of the dataset. Although less relevant to this work, let us mention [23 and [24] for the sake of completeness. Under frequentist assumptions, they also propose DBO algorithms that forget irrelevant observations by either resetting their datasets or by using decreasing covariance weights. However, they assume that the variational budget of the objective function is fixed, which
qi (2; 12) and q3 = (zz,. t2) by (1 -
has the drawback of requiring the objective function to become asymptotically static. This is a very different setting than the one of interest in this paper, which does not make this requirement.
The aforementioned algorithms all work with discrete, evenly-spaced time steps. This setting simplifies the regret analysis of DBO algorithms through the use of proof techniques similar to the ones used for static BO. However, also overlooks a critical effect of the response times of their algorithms. In fact, the response time of BO algorithm heavily depends BO inference is in O(n)). Although it is reasonable to ignore this for classical BO because the objective function is static, DBO algorithms cannot make this simplification as it directly impacts their ability to track the optimal argument of the objective function through time. Many algorithms (e.g-, see (1 23 24] recommend keep all the collected observations in their datasets, whereas in practice, their response times would asymptotically become prohibitive. Other al algorithms (e.g.,
its dataset size n, since
see [21 22 23 propose to reset their datasets, either periodically or once an event is triggered. This probably deletes some relevant observations in the process. More importantly, these algorithms necessarily estimate their covariance function hyperparameters beforehand and keep them fixed during the optimization. This lack of adaptivity of the estimation might lead to severely under-optimal characterization of the function by the hyperparameters, especially when optimizing an ever-changing objective function on an infinite time horizon.
To the best of our knowledge, only one work acknowledges these problems. It proposes ABO 20 an algorithm that uses a decomposable spatio-temporal covariance function k((z,, f), (x',²") E'l) to accurately model complex spatio-temporal correlations and samples
ks<[ æ' (()))) the objective function only when deemed necessary. Although this reduces the size of ABO's s dataset, ABO does not propose a way to remove stale observations, it only adds new observations less frequently. Therefore, using ABO will still become prohibitive in the long run.
The most relevant methods to quantify the relevancy of observation can be found in the sparse GPs literature (e.g., see [25] (6)). However, they require non non-trivial adjustments to account for the particular nature of the time dimension. As far as we know, there is method in the DBO literature able to quantify the relevancy of an observation in an online setting. As mentioned before, such a method is much needed as it would allow a DBO algorithm remove stale data on the fly while preserving the predictive performance of the algorithm. We bridge this gap by providing a sound criterion to measure the relevancy of an observation and an algorithm exploiting this criterion to remove stale data from its dataset.
3 A Wasserstein Distance-Based Criterion
3.1 Core Assumptions
To address the DBO problem under suitable smoothness conditions, let us make the usual assumption of BO, using a Gaussian Process (GP) as surrogate model for f (see (16).
Assumption 3.1. f is GP (0, ((,,, 1) (æ', I')), whose mean 0 (without loss of generality) and R..
whose covariance function is denoted by k S x x S
3

---

In order to accurately model complex the decomposition and isotropy in time and space of the covariance function is as 22] Assumption 3.2.
dynamics, we make the same assumption on
k((x,, (2', {)) = Aks x' 12 (t)
(5)
with O, kg R+ + [0, 1] and by .R+ + [0. 1] two positive correlation functions, parameterized by lengthscales Is > 0 and 1T > O, respectively. The factor À > scales the product of the two correlation functions and hence, controls the magnitude of the covariance function k. The lengthscales Is and Ir control the correlation lengths of the GP (see 27 for more details) in space and in time, respectively.
Although the covariance function k is able to model temporal correlations with kr. it does not accurately measure the relevancy of an observation. The next section addresses this question.
3.2 Measuring the Relevancy of an Observation
By definition, when an irrelevant observation gets removed from the dataset, the GP posterior experiences hardly any change. Therefore, we propose to measure the relevancy of an observation o, = ((z,; t4), 9i) by measuring the impact that the removal of has on the GP posterior.
Let GPD be the GP conditioned the dataset D {((2,; (4), U/ [e[L... with (at, 8 T and = f(æ;, (i) + E, ~ N (0, Without loss of generality, let us measure the impact of removing ((11,t)) 31) from the dataset on the GP posterior. Clearly, the measure should be defined on the domain of future predictions at time lo. denoted by Ftg. which must include the whole space and only the future time interval (Lo, +o0):
Fro S +00).
(6)
We compare a GP conditioned on the whole dataset, denoted by GPp, with a GP conditioned on D, the dataset without 101,41,11 denoted by gP For an arbitrary point (x, E Ftg. GPp provides a posterior distribution (D( = N (1D( (c), and so does GP with (z, t) = NV (#5( t), (2x, We compare these two distributions by using the 2- Wasserstein distance [28 given by
oB(a)))))) W2 (N (), ((ND( ((B((, t))² -
(7)
A natural extension of the 2-Wasserstein distance from point (æ, t) Ftg to the domain Fi, is
W₂ (CP (PP,,Pp)) - (fs
w2 (No(a, dadd
dardi
(8)
Observe that (8) is a criterion that effectively captures the impact of removing the observation
01 = ((x1, (1.,11 from the dataset on the GP posterior. However, discussed in Appendix (F] the
covariance function hyperparameters 0 = (A, ls, !T) control the magnitude of g This is illustrated by Figure which depicts couples of GP posteriors that achieve the same value d8). Depending on the lengthscale magnitude, the posterions may be quite different or, conversely, very similar. As a result, ( cannot be directly used gauge of observation relevancy.
To remove this ambiguity, we normalize (8) by W2 GP&) (i.e., the 2- 2asssrsstni distance between the GP conditioned on D and the prior GP), and obtain the ratio
We(PPo.Ppp)
R(GPD.9P3) = Wa(CPp.GPo)
(9)
Intuitively, W2(GPD, GPs) measures the impact of resetting the whole dataset D on the GP posterior and serves as a baseline that puts into perspective the distance measured by (8) Technically, taking the ratio (9) successfully cancels out the influence of the covariance function hyperparameters on the magnitude of the Wasserstein distances, as further discussed in Appendix F As direct consequence,

---

Wa(GPo, GPp) 0.46
WsGPPn GP 0.46
1.3
Pn
GP(
1.5
1.0
Pi
1.0
0.5
0.5
0.0
0055
0.5
-1.0
1.0
Up
WxGPP, GP
-1.
1.5
Pp
0.0
2.5
3.0
7.5
10.0
0.0
2.5
5.0
7,5
10.0
Figure 1: Similar values of Wasserstein distance, different effect on posteriors. For visualization purposes, only the posterior means of two posterior GPs (blue for μp and orange for PD) are depicted, along single dimension (e.g., time). The Wasserstein distance between the two posteriors shown by the green shaded area. The GPs have small lengthscale (left) or, conversely, a large lengthscale (right) for the chosen dimension.
Small lengthscale
Large lengthscale
C.P,
0.23
WaGPP 0.23
W.GPPo.GF3)
(C (1A)
1.5
MD
W/[GPo,
1.5
Pb
Wy(GPp, GPs)
Observation is mildlly relevant
1.0
1.0
0.5
0.5
0.0
0.0
-0.5
-1.
1.0
WsGPPo,GD)) WsGGPo,
-.5
2.5
5.0
7.5
10.0
0.0
2.5
5.0
7.5
10.0
WaIGR, GPa)
0.79
0.79
CP.
1.5
pD
(IGPo,GPD) W/G
1.5
PD
Wg.GPb.GPp
Observation is very relevant
1.0
Hb
1.0
Pb
WyGPPo.PP)
0.5
0.5
1.0
0.0
-05
-1.0
1.0
-1.5
-1.5
2.5
5.0
7.5
0.01
2.5
5.0
7.6
10.0
Figure 2: Normalized Wasserstein distances. Similarly to Figure ! a few couples GP posterior means (HD. P.B) depicted. The top (resp., bottom) row depicts couples of posteriors that yield small (resp., large) ratio (9). The left (resp., right) column depicts couples of posteriors controlled by small (resp., large) lengthscale. The prior GP mean A4 0 is shown black dashed line, and the Wasserstein distance between the posterior and the prior as a gray shaded area.
5

---

p is an unambiguous indication of how relevant an observation is. This is illustrated by Figure 2 which depicts couples of GP posteriors under different contexts. When (9) is small (resp., large), the posteriors are similar (resp., dissimilar) regardless of the magnitude of the lengthscale.
Exploiting this criterion is straightforward. When p is small, can infer that the observation Oj can be safely removed from the dataset since it will have virtually no impact on the posterior Conversely, when ( is large, one can infer that removing or would alter the posterior too much, and conclude that it is a relevant observation that must remain in the dataset. The exploitation of ( is discussed in more details in Section
The criterion ( is useful for DBO algorithm if and only if it can be computed on the fly, in an online setting. In the next section, we show that (9) can be approximated efficiently, and we describe a DBO algorithm able to exploit the criterion.
Using the Criterion in Practice
4.1 Computational Tractability
In [29], the authors provide an algorithm to approximate the 2-Wasserstein distance between two GPs up to arbitrary precision level. However, the computational cost of this algorithm is too expensive in an online setting, where it is crucial to keep the per- per-iteration cost as small as possible to ensure high sampling frequency. In this section, put this issue to rest by deriving explicit approximation (g. These formulas are computationally cheap enough to be exploited on the fly.
In Appendix A we show that (7) can be computed efficiently. Next, in Appendix B we apply these results to obtain an upper bound of (8) The key observation for deriving this result approximate the integrals in (8) by a convolution of the covariance functions with themselves in space and time. The same trick can be used for approximating W2(GPp, pP.)).
Theorem 4.1. Let to be the present time and D {((x,, t4), M))aee be dataset of observations made before to. Let D = {((x;, ti), y.) [.).e be the dataset without the first observation and Fr, =Sx (00) be the domain future predictions. Then, an upper for (GPpp. on Fi, is
W3 (CPD. = x² 2 +
+ >2 + tc)(((x).))),))
(10)
+ + tr (bb T + 3) C(D,D))
where C(X,)) = ((ks * k - xx,) (lT kr 1t 10,,,,0 (t, t.)) (mx,)))xx where (f * g) denotes (æz d))
the convolution between f and g, and (f 9)% denotes the convolution between and g restricted to the interval [a, b). The terms a, b, E and M are explicited in Appendices A and B
Moreover, upper bound for (GP
on Fio is (D,
W2 (GPD. = X2
(11)
This theorem provides the analytic form of an upper bound for the Wasserstein distance between GPD and GP z and the Wasserstein distance between GPD and GP, the domain of future predictions Fro Using it, we can compute an approximation Át of the relative criterion (
that
W2(5PD,9Pb) MISPD.9Pb) We(SPo.GPu
(12)
In Appendix C we study the error between the criterion 2 and the approximation 12) In essence, we bound the approximation eттor caused by estimating the integrals in (8) by self-convolution of the covariance functions ks and ky (i.e., kg kg and ky k.). Furthermore, evidence that the approximation errors in the numerator and the denominator of (12) compensate
provide numerical
each other at least in part, making (1) decent approximation for (9).
In practice, the upper bounds 10) and TT) can only be computed efficiently if the convolutions the covariance functions can themselves be computed efficiently. The analytic forms for the convolution

---

Table 1: Usual covariance functions. 1' is the Gamma function and K2 is a modified Bessel function of the second kind of order v.
Covariance Function
ax)))
1| i
Squared-Exponential (I)
22
Matérn (V,1)
2 r(w)
[20l ls
K2 [220||a]12
of two usual covariance functions listed in Table T namely Squared-Exponente (SE) and Matérn 30 are provided in Appendix [ together with Tables 3 and 4 that list the self-convolutions for the spatial (resp.. temporal) covariance function. Their detailed computations are also provided in this appendix. In nutshell, the formulas are obtained first in the Fourier domain by computing the square of the spectral densities of the covariance functions, and next by computing their inverse Fourier transform.
Together, Tables 3 in Appendix D and Theorem 4.1 show that the approximation of 2 given by (12) be computed efficiently in online setting. In effort to generalize our results to class of covariance functions that extends beyond Assumption 3.2 we also discuss how to compute the self- convolution of an anisotropic spatial SE covariance function in Appendix B We now leverage (12) to propose a DBO algorithm able to pinpoint and remove irrelevant observations in its dataset.
4.2 W-DBO
The metric (9) and its approximation a percentage, that separates GPD and GP z.. Indeed, the Wasserstein distance W (GPp, GP&)) is scaled by the Wasserstein distance W (GPp. gPp). that is, the distance between GPv and the prior. In other words, (9) and its approximation (12) measure the relative drift from GPp to GP caused by the removal of one observation. When removing multiple observations, the relative drifts naturally accumulate in a multiplicative way (similarly to the way relative errors accumulate). As consequence, removing multiple observations could, in the worst case, make GPp b drift exponentially fast from GPp. To keep this exponential drift under control, one can use removal budget b(t) (1 + a)' that allows a maximal relative drift from GPD of a per unit of time (e.g., if a = 0.1, the allowed maximal drift is 10 % per unit of time). The cost of removing an observation is given by (
12] can be seen relative ertor (or drift), expressed as
Algorithm describes W-DBO, a DBO algorithm exploiting (12) remove irrelevant observations on the fly. As described above, the removal budget is controlled by single hyperparameter and grows exponentially as time goes by (see line 24). At each iteration, [12 used to compute the relevancy of each observation in the dataset (see lines 10-13). The relevancy of the least relevant observation then compared to the removal budget, and the observation removed if the budget allows it (see lines 14-17). This process is repeated until all the budget is consumed. Such greedy observation removal policy causes W-DBO to overestimate the impact of removing multiple observations We discuss and motivate the expression of the removal budget in Appendix G The sensitivity analysis conducted Section 5.1 supports this removal budget, by showing that the same value of the hyperparameter is valid for large set of different objective functions.
Finally, note that using 12 to remove irrelevant observations on the fly can be performed in conjunction with any BO algorithm, because it can be appended at the end of each optimization step as a simple post-processing stage. This agnostic property of W-DBO is supported by the ability of Algorithm|| to take as inputs any GP model GP and any acquisition function şp. Similarly, observe that lines 5-8 in Algorithm describe the usual BO optimization loop, without any modification.
Numerical Results
In this section, study the empirical performance of W-DBO. To measure the quality of the queries made by the DBO solutions, we compute the average regret (lower is better). For the sake of realistic evaluation, two iterations of solution are seperated by its response time (i.e., the time taken to infer its hyperparameters and optimize its acquisition function). Furthermore, all covariance function
'To prevent this, the criterion 12 could be computed on every element of P! at each iteration. Unfortu- nately, this policy does not scale well with the dataset size (D).
7

---

Algorithm W-DBO
1: Input: GP, acquisition function so, hyperparameter a, clock C, initial observations D.
2: Get current time from clock
3: b =
4: while true do
5: Find æ; arg max, (xx,
6: Observe 1/, [(xt () +e, e ~ (2),
7: D= DU ((22, yb)]
8:
Condition GP on and get the MLE parameters (A, Is, 62)
9:
while b do
10:
for all (x,, ti, 81) E D [in parallel] do
11:
D = D {(2,, 14. v/.]
12:
Compute R2 using (12
13:
end for
14:
i' arg min,e]L, R
15:
i 1+ Rş- then
== {( De-)
6/ + R₂-)
18:
else
19:
break
20:
end if
21:
end while
22:
t'
23:
Get current time from clock C
24:
b= b(1 + a)(ttt))/zz
25:
end while
Eggholder2d
Sensitivity Analysis
1.0
StyblinskiTangtdd
Rosenbrock3d
g.0.331
a=1.0
Epgholde
Hortmann3dd
0-00622
0-2.0
Resirigin5d
Powellid
0-0.25
a-4.0
0-8.0
Normalized Average Ragret
0.8
Ack ey4d
Schwerel4d
Hartmannöd
Tamparaturedd
Average Regret
350
Q-0.5
Sheiel
Aug. Perf.
300
250
0.2
200
0.0
2
2
20
2
2-4
2-2
22
Hyperparametera
Hyperparameter
Figure 3: (Left) Sensitivity analysis on the Eggholder function. (Right) Aggregation of sensitivity analyses of W-DBO made on 10 synthetic functions and real-world experiment. For aggregation purposes, the average regrets in each experiment have been normalized between 0 (lowest average regret) and 1 (largest average regret). The average performance of W-DBO over all the experiments is shown in black. Standard errors are depicted with colored bars (left) and shaded areas (right).
parameters and the noise level are estimated on the fly. Please refer to Appendix H.I for further experimental details, and to Appendix Hor for detailed description of the dynamic benchmarks.
5.1 Sensitivity Analysis
We start by studying the impact of the W-DBO hyperparameter on the average regret. Recall that we take into account the response time of algorithm. This evaluation protocol reveals that trade-off must be achieved between having accurate model of the objective function (which requires large dataset) and being able to track the optimal argument the function as it evolves (which requires a high sampling frequency, thus moderately midzed dataset).

---

Table 2: Comparison of W-DBO with competing methods. The average regret over 10 independent replications is reported (lower is better). The performance of the best algorithm is written bold text. The performance of algorithms whose confidence intervals overlap the best performing algorithm's s confidence interval is underlined.
Experiment (d + 1)
GP-UCB
R-GP-UCB
TV-GP-UCB
ET-GP-UCB
ABO
W-DBO
Rastrigin (5)
17.81
53.67
26.50
19.54
36.16
18.54
Schwefel (4)
469.10 10
954.03
520.40
428.97
662.32
290.34
StyblinskiTang (4)
18.83
45.82
15.74
22.16
58.40
13.04
Eggholder (2)
542.53
273.60
287.01
559.61
256.92
225.68
Ackley (4)
4.10
4.45
3.27
3.96
3.63
2.24
Rosenbrock (3)
31.37
25.99
17.55
28.79
171.04
3.81
Shekel (4)
2.56
2.21
2.03
2.70
2.06
1.72
Hartmann3
1.17
0.26
0.82
1.06
0.55
0.35
Hartmann6 (6)
1.33
1.25
0.44
1.46
0.61
0.32
Powell (4)
1992.1
1167.6
1223.4
534.2
9888.6
428.1
Temperature (3)
1.02
0.69
1.36
1.25
1.21
0.68
WLAN (5)
1.46
4.84
1.33
4.98
12.94
1.19
Avg. Perf.
0.48
0.47
0.29
0.54
0.62
0.01
Ackley4d
Ackley4d
5-
ft
Battaee Size
10a
101
GF-UCB
R-GP-UCB W-DBO
GP UCB
G.
ABO
GVB
ET-GP.UCB
ET GP UCB
W DBO
100
200
300
400
500
600
100
200
300
400
500
600
Time
Duration (S)
Figure 4: (Left) Average regrets of the DBO solutions during the optimization of the Ackley synthetic
function. (Right) Dataset sizes of the DBO solutions during the optimization of the Ackley function.
To study this trade-off, we compare the performance W-DBO with several values of its hyper- parameter a, illustrated by the left side of Figure 3 We apply this protocol on 11 different benchmarks (described in Appendix H.2). The aggregated results (see the right side of Figure 3) show that achieving trade-off between the size of the dataset and the sampling frequency can significantly improve the performance of W-DBO. Clearly, the sweet spot reached for a 7. This hyperparameter value is used to evaluate W-DBO in the next section.
5.2 Comparison with Baselines
The competing baselines are the relevant algorithms reported in Section namely R-GP-UCB and TV-GP-UCB (2I). ET-GP-UCB (22 and ABO [20]. We also consider vanilla BO with the GP-UCB acquisition function. which only considers spatial correlations. For comparison purposes, all the results are gathered in Table 2 The benchmarks comprise ten synthetic functions and two real-world experiments. All figures (including standard errors) are provided in Appendix HH.2 and the performance of each DBO solution is discussed at length in Appendix H.3 Furthermore, the
provided supplementary animated visualizations discussed in Appendix H.4
this section. we only depict the performance in Figure 4 because it illustrates best the singularity of W-DBO. The Ackley function is known for its almost flat outer region (with lots of local minima) and deep hole the center its domain. Observe that most DBO solutions miss that hole, as their average regrets skyrocket between 200 and
DBO solutions on Ackley synthetic function

---

Aggregated Results
1.0-
StybllinkTTangd
Y0 Sheleel
Hertmannde
Egghoider22
Rasenbrock3d Hertmennndd
lemperature.dd
RastriginSd
WLAN5d
Schweft 4d
Powell4d
Aug. Port.
AckleyAc
Normllzedd Aperage
0.6
0.4
0.0
GP-UCB
ABO
ET-GP-UCB R-GP-UCB TV-GP-UCB
W-DBO
Figure 5:
Visual summary of the results reported in Table For aggregation purposes, the average regrets in each experiment have been normalized between O (lowest average regret) and 1 (largest average regret). The average performance of the DBO solutions is shown in black.
400 seconds. In contrast, W-DBO manages to rapidly exploit the hole the center of the function domain, thereby maintaining a low average regret.
This performance gap be explained by studying the dataset size of W-DBO (see the right side of Figure At first, the dataset size increases since most collected observations are relevant to predict the outer region of the Ackley function. After 200 seconds, the dataset size plateaus as W-DBO begins to realize that some previously collected observations are irrelevant to predict the shape of the hole that lies ahead. Between 300 and 400 seconds, the dataset size is halved because most previously collected observations are deemed irrelevant. Eventually, after 400 seconds, W-DBO explores the flatter outer region of the Ackley function again. Consequently, its dataset size increases again.
For summary of the performance of the DBO solutions across all our benchmarks, please refer to the last row of Table and to the visual summary in Figure [ In our experimental setting. ABO and ET-GP-UCB obtain roughly the same performance as vanilla BO. R-GP-UCB shows slightly better average performance than GP-UCB, while TV-GP-UCD appears significantly better than the aforementioned algorithms. Remarkably, W-DBO shows significantly better performance than TV-GP- PCB and outperforms the other DBO solutions by comfortable margin. In fact, it obtains the lowest average regret on almost every benchmark.
6 Conclusion
The ability to remove irrelevant observations from the dataset of DBO algorithm essential to ensure a high sampling frequency while preserving its predictive performance. To address this difficult problem, we have proposed (i) a criterion based on the Wasserstein distance to measure the relevancy of an observation, (ii) a computationally tractable approximation of this criterion to allow its use an online setting and (iii) DBO algorithm, W-DBO, that exploits this approximation. We have evaluated W-DBO against the state-of-the aae of DBO on variety of benchmarks comprising synthetic functions and real-world experiments. The evaluation was conducted the most challenging settings, where time is continuous, the time horizon is unknown, as well as the covariance functions hyperparameters. We observe that W-DBO outperforms the state-of-the-arr margin. We explain this significant performance gap by the ability of W-DBO to quantify the relevancy of each of its observations, which is not shared with any other DBO algorithm, to the best of our knowledge. As a result, W-DBO can remove irrelevant observations in a smoother and more appropriate way than by simply triggering the crasure of the whole dataset. By doing so, W-DBO simultaneously ensures a high sampling frequency and very good predictive performance.
DBO by a comfortable
In addition its impact on DBO itself, we believe that W-DBO can have significant impact on the fields that make heavy use of DBO (e.g., computer networks, robotics). As future work, we plan on W-DBO. Furthermore, we plan on better understanding the excellent
exploring these applications performance of W-DBO by addressing the difficult problem of deriving regret bound that holds in continuous time setting and incorporates the effect of the sampling frequency of the DBO algorithm well as the deletion of irrelevant observations.
10

---

References
[1] James Bergstra, Daniel Yamins, and David Cox. Making science of model search: Hyper- parameter optimization in hundreds of dimensions for vision architectures. In International conference on machine learning, pages 115-123. PMLR, 2013.
[2] Gregory Hornby, Al Globus, Derek Linden, and Jason Lohn. Automated antenna design with evolutionary algorithms. In American Institute of Aeronautics and Astronautics, 2006.
[3] Anthony Bardou and Thomas Begin. INSPIRE: Distributed Bayesian Optimization for Improv- iNg SPatlal REuse in Dense WLANs. In Proceedings of the 25th International ACM Conference on Modeling Analysis and Simulation of Wireless and Mobile Systems, pages 133 133-142, 2022.
[4] Daniel Lizotte, Tao Wang, Michael Bowling, and Dale Schuurmans. Automatic gait optimization with gaussian process regression. In Proceedings of the 20th International Joint Conference on Artifical Intelligence, IJCAI'O,, 07, page 944 949, San Francisco, CA, USA, 2007. Morgan Kaufmann Publishers Inc.
[5] Javier González, Joseph Longworth, David C. James, and Neil D. Lawrence. Bayesian opti- mization for synthetic gene design. In NIPS Workshop on Bayesian Optimization Academia and Industry, 2014.
[6] 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. Curran Associates, Inc., 2019.
[7] Johannes Kirschner, Mojmir Mutny, Nicole Hiller, Rasmus Ischebeck, and Andreas Krause. Adaptive and safe bayesian optimization in high dimensions via one-dimensional subspaces. In International Conference on Machine Learning, pages 3429-3438. PMLR, 2019.
[8] Anthony Bardou, Patrick Thiran, and Thomas Begin. Relaxing the additivity constraints in decentralized no-regret high-dimensional bayesian optimization. In The Twelfth International Conference on Learning Representations, 2024.
[9] Javier González, Zhenwen Dai, Philipp Hennig, and Neil Lawrence. Batch bayesian optimization via local penalization. In Artificial intelligence and statistics, pages 648-657. PMLR, 2016.
[10] Samuel Daulton, David Eriksson, Maximilian Balandat, and Eytan Bakshy. Multi-objective bayesian optimization over high-dimensional search spaces. In Uncertainty in Artificial Intelli- gence, pages 507-517. PMLR, 2022.
[11] Seokhyun Kim, Kimin Lee, Yeonkeun Kim, Jinwoo Shin, Seungwon Shin, and Song Chong. Dynamic control for on-demand interference-managed wlan infrastructures. TEEF/ACM Trans- actions on Networking, 28(1 97, 2019.
[12] Aurelio G Melo, Milena F Pinto, Andre LM Marcato, Leonardo M Honório, and Fabrício O Coelho. Dynamic optimization and heuristics based online coverage path planning in 3d environment for uavs. Sensors, 21(4): 108, 2021.
[13] Doyen Sahoo, Quang Pham, Jing Lu, and Steven CH Hoi. Online deep learning: learning deep neural networks the fly. In Proceedings of the 27th International Joint Conference on Artificial Intelligence, pages 2660-2666, 2018.
[14) Charu C Aggarwal, Jiawei Han, Jianyong Wang, and Philip S Yu. A framework for projected clustering of high dimensional data streams. In Proceedings of the Thirtieth international conference on Very large data bases- Volume 30, pages 852-863, 2004.
[15] Xiao Mingming, Zhang Jun, Cai Kaiquan, Cao Xianbin, and Tang Ke. Cooperative co-evolution with weighted random grouping for large-scale crossing waypoints locating in air route network. In 2011 IEFE 23rd International Conference Tools with Artificial Intelligence, pages 215-222. IEEE, 2011.
[16) Christopher Williams and Carl Rasmussen. Gaussian processes for regression. Advances in neural information processing systems, 8, 1995.
[17] Niranjan Srinivas, Andreas Krause, Sham M. Kakade, and Matthias W. Seeger. Information- theoretic regret bounds for gaussian process optimization in the bandit setting. IEEE Transac- tions on Information Theory, 58(5):3250-3266 2012.
[18] Jonas Mockus. Application of bayesian approach to numerical methods of global and stochastic optimization. Journal of Global Optimization, 4:347-365, 1994.
11

---

[19] Donald R. Jones, Matthias Schonlau, and William J. Welch. Efficient global optimization of expensive black-box functions. Journal of Global optimization, 13(4):455-492, 1998.
[20] Favour M Nyikosa, Michael A Osborne, and Stephen J Roberts. Bayesian optimization for dynamic problems. arXiv preprint 18X3::8000 2018.
[21] Iija Bogunovic, Jonathan Scarlett, and Volkan Cevher. vinying gaussian process bandit optimization. In Artificial Intelligence and Statistics, pages 3 314-323. PMLR, 2016.
[22] Paul Brunzema, Alexander von Rohr, Friedrich Solowjow, and Sebastian Trimpe. Event- triggered time-varying bayesian optimization. arXiv preprint arXiv:2208. 10790, 2022.
[23] Xingyu Zhou and Ness Shroff. No-regret algorithms for time-varying bayesian optimization. In 2021 55th Annual Conference on Information Sciences and Systems (CISS), pages 1-6. IEEE, 2021.
[24] Yuntian Deng, Xingyu Zhou, Backjin Kim, Ambuj Tewari, Abhishek Gupta, and Ness Shroff. Weighted gaussian process bandits for non- non-stationary environments. In International Conference on Artificial Intelligence and Statistics, pages 6909-6932. PMLR, 2022.
[25] Joaquin Siinonero-Candele and Carl Edward Rasmussen. A unifying view of sparse approxi- mate gaussian process regression. The Journal of Machine Learning Research, 6:1939-199 1959, 2005.
[26] Henry B Moss, Sebastian W Ober, and Victor Picheny. Inducing point allocation for sparse gaussian processes in high-throughput bayesian optimisation. International Conference on Artificial Intelligence and Statistics, pages 5213-5230. PMLR, 2023.
[27] Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
[28) Leonid V Kantorovich. Mathematical methods of organizing and planning production. Man- agement science, 6(4):366-422, 1960.
[29] Anton Mallasto and Aasa Feragen. Learning from uncertain curves: The 2- wasserssteii metric for gaussian processes. Advances in Neural Information Processing Systems, 30, 2017.
[30] Emilio Porcu, Moreno Bevilacqua, Robert Schaback, and Chris J Oates. The mat) 'ern model: A journey through statistics, numerical analysis and machine learning. arXiv preprint arXiv:2303.0275:1 2023.
[31] Izrail Solomonovich Gradshteyn and Iosif Moiseevich Ryzhik. Table of integrals, series, and products. Academic press, 2014.
[32] Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton. Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy. BoTorch: Framework for Efficient Monte-Carlo Bayesian Optimization. In Advances in Neural Information Processing Systems 33, 2020.
[33] Wenzel Jakob, Jason Rhinelander, and Dean Moldovan. pybind1 seamless operability between c++1l and python, 2017. https//ggitubbccom////ii//ybbbdd
[34] JHB Kemperman. On the shannon capacity of an arbitrary channel. In Indagationes Mathemati- cae (Proceedings), volume 77, pages 101-115. Forth-Hollane 1974.
12

---

A Wasserstein Distance at a Point in Fro
At the core of 8 lies 7 In this appendix, we provide explicit expressions for 7 proving the following lemma. Lemma A.1.
Let start by
E H
G F
(13)
with AD = D(D,D) + I and E, F, G, H defined as
= (A+ +2 0² - T (101,11),) A (111,
-1
(14)
F A A+ (101, ((0,)(101))))))))) kJ (mi, ())
(15)
G EkT (101, A
(16)
H = x Fk (1),
(17)
where &s (X,J) = (k: (i), (x), and Ab = k(D,D) 12
(x,t)ex (nyY)eyy
Proof.
The proof is trivial using the inverse of block matrix:
A B
AA
BD CD
(D-CA-1D -1 (-CA-1 I
I -BD -BD-1)
C
D
C)
- BD (A-BD--C)- BD
(18)
CD-CAA CA-1
D CA
Note that A and D must be invertible.
We
18 to write AD' function of A b since
Ap' D = (1) + o² k
4), (T(11,11))
A((x11t1),
A o -1 4 (1x1,,11)
X
(19)
k1)
A
Note that, 19 X + x² and Ad are invertible. Therefore, yields the destred result.
18
can be applied
19 and this
L
To compute 7 we provide the following results. Proposition A.2.
pD(x, t) 1,5(xxt) = ali((ar, 1), (m1, (1)) + bk ((x, f), D)
(20)
with a Eyi + Gý. bT = H y1 + "j (F-A and ÿ = (y22 (. ,yn).
Proposition A.3.
(ao(x, 1) (o:0))opp(ppp) (2l 1))2 Ek((o, (xi, (21,,)))
k((x,, 1), (x), (1),11))ck(((a(()) T ((x, (), ()
(21)
with c G+H and M = F A
We now prove Proposition A.2
13

---

Proof. According to 13)
(B( t) = k ((x, ()x,t), D),
(22)
Applying the same definition for (2 f) we have
t (21, (z1, ((z, t), D)) (A((z, (21, ((z, (), D))
Ap'y
E G H F y
(23)
(k((z, (21. t), D))
Eyy + Gý
Нyı + Fÿj
((( f), (21, t1)) Ey: + Gÿ) + k ((x, (), (Hy: + Fø),
(24)
where (23) follows from Lemma A.T
Finally, we have
( 1) (0 () (1,, (Ey) + Gtj) + k ((x1 (), D) (H + Fū)
k ((x, 1), D) 5'
which can be reduced to
pD(a, t) - (b(xx, t) = ak((2, t), (x1, + bk( ((x t), t).
(25)
with a Ey1 + Gý and bT = = y1 + "J F (F- 5') This concludes the proof.
Finally, we prove Proposition A.3
Proof. As (op(z, t) - (z, t))² is hard to integrate, to get 8 we upper bound it by
(oo(a, og(x, t)) oz t) + oz 20 (2, - c 1) + (z,.1)) - (0 1)
(26)
= - (x, 1) (z)
(27)
where 26 follows from oD(2,)) t) S (.r,t))
Now, yields that for x= D,
(x) = = - k- T=TTTT(ae ((x, t), D)) ()((e))
k((x,t), D).
(28)
and for X = =D = Du (z),
o (xx = x - (M((z,)) t), (x1 (x1, t1)), kT ((z,((z t(z, A
(((((, (11)
(k((a,,),
(),( D)
> (M((,,)), (1), (mi, T ((z, t), D))
(
E G
(29)
F)
k((x,(),
x (M((a,, (), 4)), TT ((0, 1), D)
(E,((,,)) (2),t))) + Gk(( t).
Hk((x, t), (2),))) Fk(( t), D)
(30)
where 29 follows from Lemma A.I Developing the dot product in (30). we have
(x 1) = À 2 ((x, 1), (1) - 1), ( D) <k(xx, t), (#1, H1)) k((x,t), D) - k (T t), F) D).
(31)
Combining (88) and 31). we get
( (x, t) Ek((x, (t), + k((x, t), (x), (c)c (xx,ttt) t).
(32)
+ k ((x, t), Dk) t), D)
where C-G+ H and M F F-
Combining (27 and 32 concludes the proof.
14

---

B Wasserstein Distance on Fro
In this appendix, we extend the computation of the Wasserstein distance between posterior Gaussian distributions at an arbitrary point in Ft, (see Propositions A.2 and A.3 to the Wasserstein distance between GPD and GPB on Fro. We also provide an upper bound for the Wasserstein distance between a posterior GP conditioned on D (i.e. GPb) and the prior GP (i.e, GP
Let us start by proving the following lemma.
Lemma B.1. Let to be the present time and = {((x;; (:)) be a dataset of observations before to- Let D {((2,, 14), Vi) set2.n and Fro Sx to: +00) be the domain future predictions. Then, for any pair of observations (Ii,, li, (2:j. 1,, from D we have
ff I 1: ((22,)) (a)e,
Is.
1:)) (2n,, ((x,, t), (e,,t))) dxdt S X2 (ks ka) xx) ky tor to- (ti
t.p)
(33)
with (ʃ * g) the convolution between and and (ʃ * g)a the convolution between I and g restricted to the interval (a, b).
Proof. For the sake brevity, let us denote by the LHS of (33 According to Assumption 3.2
Fubini's theorem and since Fr.- S x |to, +oo), we have
I T= x² f ks( - x. 22ll2)ks (11a: - - x3l12)d
/
kr - 41)r((1 tyl)dt.
(34)
Let us focus the integral on the spatial domain S. Depending on the expression of the covariance function, this integral can be quite difficult to compute. Let us turn this expression into a more pleasant upper bound
p ks([l x. x.(12)bs(11x xyll2)dx < of ks(11 x.ll2)ks(1 x,112)dx
(35)
= f ks
xx; - x--ull2)dde ull2)dd
(36)
35 holds since ks is positive and (36) comes from the change of variable u x; Obviously, this upper bound remains interesting because ky is usually an exponentially decreasing function (see Table We discuss this point in more details in Appendix[[
where
Observe that (36) is actually the convolution (ks kx) xx). A similar change of variable can be made regarding the time integral, with v = 44:
f +oc
kT (It 44 (e(((((( &
(T(l)))T((t - vl)dv.
(37)
Note that (37 is also convolution, but restricted to the interval [to 44, t00). We denote this convolution on a restricted interval (KT kT)
+ - (t, - 44).
Combining 34
36 and 37) yields Lemma B.1
We now prove the following lemma. Lemma B.2. Let to be the present time and D = {((x,, tx), be) be dataset of observations
before lo. Lei De [((0;, 4:), M:) vi))ie and F., = S x [to, +00) the domain of future predictions. Then, we have
$ 1 (uD(a,,) t) - (x((x,t)) dadt & XaaC(((1, (m1,t)) + 2X 222
(38)
+ x tr (bb ((D, D)).
with a = Ey1 + Gy and bt H y1 + T (F-A5 C(X =
((ks ks)(z, - (kr # kr trr (t, - t,)) (m))))) (E)
1 the conformable vector of ones,
15

---

(f 9) the convolution between f and 9 and (f * 9)a the convolution between f and g restricted to the interval [a, b).
Proof. By Proposition
A.2
f [ (uD(ar, () (3( dadd
dadd f { () + bk: ((e, 1), D)) dadl
(39)
with Eyi Gý and bJ H y1 H
+ 3" (F-
Expanding the square, we get three different integrals, denoted A, B and C with
f L d2x ddd,
A = Is
B $s I/ 2abk((æ 1), (m ()) t), D)dadt,
{
{ { bk (zx,t). k_ ((z,)), b² dedt.
C fs
Using Lemma B.1
computing upper bound for A, B and C is immediate. In fact,
A tt),
(40)
B S (1),))
(41)
C < X²1T
6b o c(D,D)
= X2 bb CDD,
(42)
where O is the Hadamard product and 1 the conformable vector of ones.
Adding 40
and
42 together concludes our proof.
To get the first part of Theorem 4.1 we prove the following lemma.
Lemma B.3. Let to be the present time and D {((2,, tx), yelvee be dataset of observations before lo. Lei D [((x,, 4:), and Fi. S x [lo, +00) the domain of future predictions. Then, we have
$ 51 (op(z,t) ()) dxdt & X² EC((
ti), (zn + X (CC((E))
(43)
xr (MC(D,D) D)).
with (Oks ks)( x. (Ko # kr tta
G + H
M
F A.. C(X
(2.))))) (x)
( the convolution between
and g and (f * g)2 the convolution between and g restricted the interval [a, 6).
Proof. This proof is conceptually identical to the previous one. By Proposition A.3
dxdd < o 1 (z) 21))² + k((x,,)), t), (21, f), D) 8s 1 (op(a, t) (2x ddxdt
$
+ + ((,,t) t), dredt
(44)
where G + H and M F -
16

---

By linearity of the integral, the RHS of (44) can be split into three different integrals A, B and C where
A=
Ek((z, (1) dxdt dxdt
Is t
B $ { k((2, t), (2), (21,41)ck((as; t), D)
S
C / I k ((a,
D)
Once again, Lemma B.1 can be used to compute an upper bound for A, and fact,
- X (1), B $ (())((
(45)
(46)
c A ((D,
tr (MC(D,
(47)
where 3 is the Hadamard product and 1 is the conformable vector of ones.
Adding 45
46 and
47 together concludes the proof.
B. and B.3 yield the first part of Theorem For the second part of Theorem we prove the tollowing lemma.
Together, Lemmas B.1
Lemma B.4. Let la be the present time and = (((a,, (e), Wi))ie( be a dataset of observations before to Let Fio =S x foo, +o0) the domain of future predictions. Then, we have
(GPP. < X2 (a (D) tr C
(48)
with A- and C(X,)) ((kg kss) ( a.:((r 1e) (1) (4)) (m.t))x (f g) (a))))y
the convolution between and g and (f g) the convolution between f and restricted to the interval [a, b].
Proof. Recall that, according to GPo. f(x,t) N (0, A) for any point (z, t) E Fio- Consequently,
(f:00-(ff.. (x l)dxdt + 111 (Vā- dwdt
Wi(PP,.PPo)
(49)
These two integrals in (49) can be computed with the same techniques as above. For the mean integral, we have
$ 1 (2, t)dxdt of 1 y A D). - "ydxd
[1 k
TT ()k((x) D)) A A-ly
to
< X C((,D))-
(50)
Regarding the variance integral in 49
we have
1110000 (V-ooo(ztt) dxdt
$ pc
f (1 2Voo ob, dædt
J 1
< $
---bbla 1)) dadd
(51)
Js
f Js
kT (((, t), D)A *k((x, t), D)
(52)
< X tr A -IC(D,D))
(53)
17

---

where (51) holds because (a, t) 2 (2, and 52 holds because of
Together, 50 and (53) conclude the proof.
By combining Lemmas B.2 and B.3
the proof of Theorem 4.1 is immediate.
C Approximation Error
In Appendix B we provided computationally tractable upper bound of W2 9P and W2(GPD, GPo) In this appendix, we provide the expression of the corresponding approximation error and we study its magnitude with the Squared-Exponentiaa (SE) covariance function.
First, and without loss of generality, assume S [0, _J In Appendix B have to approximate the Wasserstein distance because the integration over S in 34 is difficult to compute. The upper bound proposed in (35)
fo.1 ks(1l== xill2)ss(11a =zyl2)dx s f ks(l 2.1la)ks((1 xyy 2zy12)dæ.
(54)
Recall that the upper bounded quantity is product of functions which are decreasing exponentially (see Table []. As consequence, their product decreases exponentially as well, so that extending the integration from (0. 1] to Rd has bounded impact on the result. Clearly, a first absolute
approximation error for (54) is
ks(ll= - 1=ll2)ks(11== 2.y 12)122)x.
Because the upper bounded quantity product two correlation functions kis on a hypercube of volume 1, the upper bound can be capped to 1 as well. This leads to the more refined absolute approximation error:
A(x, ss)ls) = min f ks(l a=lle))s((12 12)12
(55)
o fo.1ja
- [1.1===;]la))[[(( x;il2)dx
Obtaining closed-form for the approximation error (55) is difficult. However, because the spatial lengthscale controls the correlation lengths in the spatial domain, it is clear that the left term in (55) is an increasing function with respect ls. Conversely, the right term in (55) is decreasing function with respect to (s. This observation allows us to derive the spatial lengthscale for which (55) is maximal.
Proposition C.I. Let (æ;, E S², with S = [0, 11 Let ki be SE kernel with lengthscale ls. Then,
arg max A(x;, x;; (s) Wo leERt
(56)
with Wo the principal branch of the Lambert function.
Proof. Because the two terms in 155 respectively increasing and decreasing with respect to the spatial lengthscale Is. (55) is maximal when both terms are equal. Therefore, from (55). we have the following relation
L { ks(ller 21112)bs(11z 2zyll)ddx
{ Pu. ks(1 az ksll2)ss(1 in 22il2 da.,
ks([læ a: 2ss12) -1,l12)da
-11m;-0,118
rd n
e
1,
(57)
18

---

Kernel
SE Kamai
[d-1. Kk
0.53
(dm3. 1P0 Xila
1011
109
10
10"1
10
10
10'
Aboo Emer
Absokrt Ener
10:1
Eor Temm
Em torn
1st tro Temn
2nd Eror lomn
2nd Erro Tarm
2nd Enor Turn
Ueal
Critical 5
Crucclly
0.4
1.
1.
0.2
0.4
0.66
0.6
1.0
1.2
1.4
02
0.4
0.6
12
1.4
Sporial Longthssalee it
Scatial Lengthhcale
Spatial Lamgthwoalk
SF Kernel
SE Kernel
SF Kernel
[|4 4/lle
X1l -0.51
100
10-
10"
10
Ranatvee
12%
Relit Errrr
Relaive
RE WEP 1st Temn
1011
R5 with term
10-3
RE W 2nd Tmm
Tam
BE with 2nd Torn
ocalls is
lcalls
0.2
0.6
D.
1.0
14
02. 0.0
1.0
1A
02
0.0
1.2
Spall al Lengthscale
Spetial Lengithceee
Spaltel Langtnacele
Figure 6: (Top row) Absolute approximation error (55) with respect to the spatial lengthscale (g for 1, 3 and 5-dimensional spatial domain. Both error terms in (55) are shown in orange and green dashed lines, respectively. Finally, the critical lengthscale (56 is shown as a red vertical line. this example, ks is SE correlation function. (Bottom row) Relative approximation error with respect to the spatial lengthscale Ig. The color codes are the same.
Kernel
Kernel
Kernal
Relerive Ertr
Retative Error
Ralatvo Brar BE with 1at Tarm with 2nd Term
hedebre Ere Critern
RC with 2et Turn with lkm
with 1s: Rim
Dalatbe Erooe on Citarton
0.4
0.
022
64
0.0
0.0
-0.2
-0.
Spatial Lengthscale
Spelial Lingifscale.
Spatial Lengtrsca
Figure 7: Relative error between the criterion (9) and its approximation [12 with respect to the spatial lengthscale (s for 1, 3 and 5-dimensional spatial domain. The relative error computed with both terms in (55) are shown as orange and green dashed lines, respectively. In this example, kg is a SE correlation function.
where 57 uses a result derived in Appendix and reported in Table 3 Solving for is concludes the proof.
Let us illustrate (55) and Proposition C.I by plotting the relative approximation error (55) with respect the spatial lengthscale 1s. The integral over S is computed with Monte-Carlo numerical integration technique. The results are shown in Figure 6 Looking at the top row, we see that the absolute approximation error peaks at spatial lengthscale 13 given by (56), as anticipated above. For Is < 15. the error is given by the first error term in (55). and conversely, the error is given by the second error term in (55) for Ig I3. The bottom row of Figure shows that the same observations apply to the relative errors.
Figure 6 shows that even though it is bounded, the approximation error of the upper bound (54) is non-negligible. This is particularly noticeable when looking at the relative errors in the bottom row of Figure 6 which clearly increases in magnitude with the dimensionality of the spatial domain.
Nevertheless, recall that we seek to approximate the ratio 9 with a ratio of upper bounded Wasserstein distances (12) that involve (54) both in the numerator and in the denominator: Because the spatial lengthscale T. does not vary when computing the numerator and the denominator, the errors are of similar magnitude and point in the same direction (both the numerator and the denominator are upper
19

---

Table 3: Analytic forms for the convolution of usual spatial covariance functions. I is the Gamma function, Ja is Bessel function of the first kind of order a, Ka is modified Bessel function of the second kind of order a.
Covariance Function ky
(ks (g)(
-lleli 41g
(Is)
Matérn (s)
2 -ut 2
V2v
2
|
2
K2++4
K
[|l||
2/
((r)T(2vdd)
4
t
Table 4: Analytic forms for the convolution of usual temporal covariance functions on the interval Ito ti +oo). Note that erf the error function. Also, for the sake of brevity, the terms CkG lo Pki b2 and K2 are defined in Appendix D.2.2
Covariance Function kT
(NT * liT
110 (1) 1:)
4:00
-(4-1/)2
Squared-Fxpooentia (lr)
VolT e
21/
1 (1- erf 210 (20-11-12)
-12pt
Matérn (v/ = + (.) Ehi=o Eha=o Ckikee
Priks (tostsitt))
bounds). As consequence, the approximation errors compensate each other (at least in part) when computing (12). To verify this observation numerically, we compute the relative approximation enror between the criterion (9) and its approximation (12) The results shown in Figure 7 Although the approximation errors in the numerator and the denominator do not entirely compensate each other, the approximation (12) appears to have lost most of its dependency to the dimensionality of the spatial domain and to the spatial lengthscale, making (12) a decent approximation of (9) regardless of or (s. In the main paper, Section 5 corroborates this observation by demonstrating the usefulness of the approximation 12 in practice.
D Convolutions of Usual Covariance Functions
In this appendix, we derive the analytic forms of the convolution of usual covariance functions listed in Tables and which used compute the criterion (12).
D.1 Spatial Covariance Functions
In this subsection, we compute specifically the analytic forms for the convolution of usual spatial covariance functions. We rely on a direct consequence of the convolution theorem, that (k* F-((F' (k)) (x), with F(f) denoting the Fourier transform of and FF (f) the inverse Fourier transform of f.
The Fourier transform F(k) of stationary covariance function l: is called the spectral density of k and is usually denoted Both functions are Fourier duals of each other (see (27 for more details). Furthermore, it is known that if & is isotropic (i.e. it can be written as function of r = ||x:||2), then its spectral density S(s) can be written as a function of s = ||I||2. In that case, the two functions are linked by the pair of transforms (see (27
2x
k(r)
S(s) (
(2Trs)
(58)
- 1
2
S(s)
k (2Trs)r
(59)
where Jg -1 is Bessel function of the first kind and of order d/2 - 1.
As an immediate consequence of 58 (59) and the convolution theorem, we have the following corollary.
20

---

Corollary D.1. Let k be a stationary, isotropic covariance function with spectral density S. Let r =|||||| and || Then,
2T (k* k)())r)
7 -
(60)
where -1 is a Bessel function of the first kind and of order d/2
now derive the analytic forms for the convolutions of usual spatial covariance functions kg.
D.1.1 ExpenentaEl Covariance Function
Lemma D.2. Let kg be Squared--xppoenntnn covariance function (see Table 7 with lengthscale Is > 0. Then,
1
(ks * (x) = -
"3
(61)
Proof. The spectral density of a Squared-Exponential covariance function kg is (see 27
S(s) (2mlz)
(62)
According to Corollary D.I
we have
2T
(ks (x)
r
2T
(2xl2) 100-mmm1
1n (
J (2
r
(2T)
1dd
000
1gs
J (2 ds
(63)
where r ||0:||2 and Jg. 1 a Bessel function of the first kind of order d 1.
is known (see 31 that
ax² (rttJo(Bx)de
8
2²
mm
(2xx)
Therefore,
(2 o+1 12d
(2 -- 2,2,3
(ks * ks)()))
16
r
22)) 4
ige 43
(64)
Replacing r by ||æ||| in 664 concludes the proof.
D.1.2 Matérn Covariance Function
Lemma D.3. Let kg be a Matérn covariance function (see Table 7.. with smoothness parameter > 0 and lengthscale Ig >> 0. Then,
24 - - dy?
12vv
2v-
[Ieella/22 (65)
(kis (x)
||
K2u+
K
T(v))T(2u
1s
1g
where I' is the Gamma function and Ka is modified Bessel function of the second kind of order a.
Proof. The spectral density of a Matérn covariance function kg is (see 27
T (2v)
12v + 4m 22
S(s)
(66)
[(v)T2"
P2
21

---

where I' is the Gamma function.
According to Corollary D.I
have
2
(ky kz))
r (S) - (2Trs)
2x
00 2nd I + (2
2v
-2--d
4m²
- (2trs)s
r
[2(v)
12
22d+ a T2 + (20) 2M
2v
+ 4m 8
(2 do
(1) p2 -1
12
2 T (2v) 20
2v
- 2v
u2
J. (ru) du
(67)
("(
22
where J is Bessel function of the first kind of order variable 2TS.
2 and 67 comes from the change of
is known (see 31) that
(0 ( (a dr
a0-
2n 1) K (aba(ab)
where Kμ-a is modified Bessel function of the second kind of order μ a. Therefore,
(ks (s)(z)= 2 T (20)
(2)))
V2y
22
y 2u+
(rv2v
Kavt
12(v)13"
ls
Wootd-I[(21) + d)
1s
2 - 20+1 π 1 T(v + 5)2
V22
20-
rV22v
2 2 K2ut
p
(68)
((v)TT(2v d)
is
1s
Replacing r by ||0||| in 68 concludes the proof.
D.2 Temporal Covariance Functions
In this section, we derive the analytic expression of the convolutions of the most popular temporal covariance functions ky restricted the interval it t,, to0). Therefore, compute many integrals of the form
kr kr(tt)brrtt -t +|)dt, fo-4
which can be rewritten
kr( (1+ tz)dt.
(69)
1o-
since 2 0 and 1/ L. for all L E (lo 1i +00). The form 669 will be used in every proof of this section.
D.2.1 Squared-Exponentte Covariance Function
Lemma D.4. Let ko be a Then,
covariance functtin (see Table 7 with lengthscale
Is >
Vāly -(4-1,) 24
(1 (1-erf) --erf" 221 4
220-4411 444)
(Jer * ke t000 (t4 12)
(70)
2
2hr
where erf is the error function.
22

---

Proof. Since by is a sennared-xxponentide function, 69
becomes
-(1-1,t1))2
e.
2
6
2
dt
to-
+00
-12 2008-20 -200-100000 10 2525)
2/
dt.
(71)
14-4
It is known (see 31 that
(ax +oa2++2bat)) dr
2
1 π
8-ss
erf
b (Vxx
+
Va
where erf the error function. Therefore,
71 becomes
V n
-14-1p2 21/
(1 -r-f
2l0 - 1, 1j
21r
D.2.2 Matérn Covariance Function
Lemma D.S. Let kr be Matérn covariance function (see Table 77. with smoothness parameter p= + 3.P E N and lengthscale bT > Then,
(le * kr kor)a-t. tto tahoo-te tto- *(t - 11) = Σ
Ckyke
(to.titt)
(72)
R1 1==
where
p!
2
(p + kI) K2)!
(2/2p +
k1 12 k2-1
CRi
(73)
(2p))
by 1k:
k1 k1) Kz)!
4
2p 2p-k1 kg
br
pika) (to tt),
(74)
2v22+1 +I
P(t) = t" ti tÚ + tt) and p(k) the kth derivative of P(L) with respect to
(75)
Proof. The Matérn covariance function has simpler form when its smoothness parameter v is half-integer, that p+ + 2. f E (see 30) ). In that case,
kT(()) = -IFTIL pl
(p+ k2)!
12/2+11
p-hi
(2p)
kp kj)!
k1=0
Therefore,
2v +1
P
P
(t 12)
Σ
C
'7
P(t)
(76)
lr
k2=0
with Ck,k, defined in 73 and P(t) defined in Integrating (6).. get
775
2v +
too
-2v3rt1t
Ch
28
P(t)dt
(77)
1T
A=
t0-ti
thanks to the linearity of the integral.
23

---

It is known (see 31 that
I
6a00
P.R (x)
P(t)e do
a
k=0
where Pm is a polynomial of degree and ((l) is the keth derivative of P... Therefore,
f
11 bp
lr
2v2pt bp
2p ky kka
(t)
P(t)dt-
(78)
2v
-I
(2 2p
Combining 77 and 78 get
P
2a Ti2to
14)
(KT - (ts - ty)
Ckikge
Pke (to, ti,tt)
=== k2=
with Ckyka defined in 73 and Pkika defined in
74
This concludes the proof.
E Extension to Anisotropic Spatial Kernels
In this appendix, we illustrate how Theorem 4.1 could be extended anisotropic spatial kernels by considering an Automatic Relevance Detection (ARD) gquared-Exponentie (SE). It has the following form:
(ar y) = - (x=-)) M (x=-))
-
(79)
where M diag (1....... (2) is diagonal matrix that gathers different lengthscale for each dimension. Observe that the isotropic SE with lengthscale Is is retrieved by setting M = IsI.
Because the ARD SE kernel 79 is anisotropic, the convolution with itself
(ks * ks) y) =
off las 2) =)dz
(80)
cannot be simplified to a one-dimensional integral through a change to polar coordinates, as done in Corollary [.]] The integral becomes more complex, but can still be computed exactly for some kernel such as the ARD SE.
Lemma F.1. Let ké be ARD SE covariance function with parameter M. Then,
(ks kx) y) o det (M) e (a) M ((x-v)
(81)
Proof. For the ARD SE kernel with parameter M diag (4, (.. .4), the convolution (80) is +(== (v--))
(ks * ka y) = of (== M (x-
c { P (he
c = M Pztz M 2 dz
(ce uy of Eiin (mx ax) / 1-1-111011 ... dzd
(#M M 2n)
II
e (mx - (u - -2)/
(82)
km
Integrating the k-thh term in 82 we get
e cul ) = -
Vilge
(zg. 1gg
41g
erf
+
I& + Uk
t00
2
21
Loa 1
Tlpe
41
(83)
24

---

Injecting
83
into
82
we have
(ks ks)( y) e - 1 3 ²8) II VTlue
"
श det (M) e (a+u) ((+t)) "a "att
(84)
π det (M) ((m-v)" M (a--)))))
(85)
where (84) holds because the determinant of a diagonal matrix is the product of its diagonal elements.
As a safety check, observe that Lemma D.2 is special case of Lemma E.1 where M = Ig that is, when ky an isotropic SE kernel.
Relative Quantification of Relevancy
In this appendix, discuss how and its approximation (12) address the dependency the covariance function hyperparameters introduced by (8). For the sake this discussion, we take kş and by two Squared- ExtiareedaEoonntiialll (SE) covariance functions (see Table T. A similar reasoning can be conducted with Matérn covariance functions.
Let us start by rewriting the product of spatial and temporal convolutions C((æ, t), (a', (()) with the formulas provided in Tables 3 for the SE covariance functions. We get
C((ar 1), (()) π 12
-Il-- U3
VT
lre
-t-e
21
erf
2lo 1
2
2/r
1
11e-o'18
-(e-1122
2to - t
π
Isbre
45
217.
1 erf
2
2/T
1 d+1
(x),t)
(86)
2
The dependency on the covariance function hyperparameters X, (s. br appears clearly in (86 Both (s and br are used, not only to scale the spatial distance ||x: æ'112 and the temporal distance
in
C", but also as a scaling constant of the magnitude of the output of the product of convolutions itself.
Because C((x, f), (ar', (')) is involved in every term of (10) and 11] in Theorem 4.1
2 1gly
be factored out of 10) and 11 Overall, both equations have in common the factor 2 * X
Clearly, this shows how the covariance function hyperparameters 0 (A, (7, may control the magnitude of the Wasserstein distances.
To capture the intrinsic relevancy of an observation regardless of the hyperparameters values, one can compute (12), that is the ratio between (IO) and (II). Doing so, the factors which are common the two equations cancel out. Considering the application of Theorem 4.1 with ky and ky being SE covariance functions, the undesirable factor 1 d+1 >²13/r is removed. Clearly, (12 remains function of is and ly, but the hyperparameters only used scale the spatial and temporal distances, that to control correlation lengths. However, the undesirable scaling exposed in (86) no longer exists.
G Removal Budget
In this appendix, discuss why the removal budget of W-DBO (see Algorithm D. denoted b, at given time f, has the form
bo
= 1,
(87)
bit
be(1 + a At/lr
with by the temporal lengthscale (see Assumption 3.2p and a the hyperparameter of W-DBO. Crucially, by and must be expressed in the same unit of time.
25

---

Table 5: Comparison of removal budgets (87) and (88) when doing experiments of different durations
on the Hartmann3d synthetie function. All experiments use the same time domain [0, 1].
Duration D
Second (axis unit)
Lengthscale br (axis unit)
Lengthscale lr
Budget (7) (1+ a) byt
Budget (8)) (1 a a)"
(seconds)
(seconds)
300
1/300
3/5
180
(1 + a )6/3
(1 + a )300
600
1/600
3/5
480
(1 + a )6/3
(1 + a 600
1800
1/1800
3/5
1080
(1 + a)//3
(1 + a) 1800
First, note that the expression of the budget is intuitive because 9 measures ratio, expressed
percentage. Therefore, the budget must accumulate in multiplicative way, leading to the exponential form (87p.
More interestingly, let us discuss the exponent At/br. Arguably, an alternative, more intuitive form of the removal budget would be
bo =1, (itat (1+++ a)
(88)
Although easier understand, the budget (88) presents major problem since it depends arbi- trary choices made by the user. This is illustrated by Table [] where the same synthetic function (Hartmann3d) is optimized under three different durations. When the duration varies, the removal budget (88). which depends on the number of elapsed seconds only, also varies. Conversely, the budget (87) remains the same. This is because the number of temporal lengthscales elapsed during the experiment remains constant, regardless of the experiment duration.
Using the removal budget (88) becomes really troublesome when it comes to making recom- mendation for the hyperparameter a. the analysis in Section [5. had used the budget (88), its recommendation a" would have been a function of the temporal lengthscale, and it would have been valid only for experiments with the same duration (e.g., ten minutes). Any other experiment duration would have required another sensitivity analysis.
Conversely, the recommendation made in Section 3.T. using the budget (87). a single number that is valid regardless of experiment duration. This is much more general insight.
H Empirical Results
H.1 Experimental Settings
In each experiment, the d-dimensional spatial domain is scaled in S' = [0, 1j and the temporal domain (viewed as the (d + 1)th dimension) is normalized in [o, 1]. Additionally, each optimization task lasts 600 seconds (10 minutes).
Unless stated otherwise, each DBO algorithm exploits a Matern- 5/2 kernel as its spatial covariance function. GP-UCB, R-GP-UCB and ET-GP-UCB do not explicitly take into account temporal correlations, while TV-GP-UCB uses own temporal covariance function. Eventually, ABO and W-DBO exploits Matern-3/2 kernel as their temporal covariance function.
Each DBO algorithm begins its optimization task with 15 initial observations, uniformly sampled in S' x [0, 40 - At each iteration (at time t), (i) the noise level as well as the kernel parameters are estimated, and (ii) the GP-UCB acquisition function is optimized to get the next query. The sum of the times taken to perform tasks (i) and (ii) is the response time of the DBO algorithm, denoted by At. Clearly, At is function of the dataset size of the DBO algorithm. Consequently, it varies throughout the optimization, getting larger when the DBO algorithm adds new observation to its dataset, and getting smaller when the DBO algorithm removes at least one point. Once (i) and (ii) are performed, the objective function is immediately sampled (except for ABO which can decide to sample / at specific time in the future) and Gaussian noise with variance equal to 5 % the signal variance added. Then, the next iteration begins at time t At (except for ABO if it decides to sample f later).
1
26

---

Rastrigin5d
Rastrigin5d
60
GP UCB
GP UCB TV-GP-UCB W-DRD
GP UCB
G UC3 IV-GRUCE
ABO
ABO EI-GRUCB
FT-GPAICA
WIBBO
50
Average Regres
10?
40
Datazel
30
10
2
2
100
200
300
400
500
600
Average Response
Duration (5)
Figure 8: (Left) Average response time and average regrets of the DBO solutions during the opti- mization of the Rastrigin synthetic function. (Right) Dataset sizes of the DBO solutions during the optimization of the Rastrigin synthetic function.
For the sake of benchmarking fairness, all the solutions have been implemented using the same popular BO Python library, namely BOTorch [32] (MIT License). To comply with the technical choices (i.e., Python front-end, C++ back-end), the computationalll-heen part of W-DBO (i.e., the evaluation of the formulas in Section 4 have been implemented in C++ and bound to the Python code with PyBind11 33 (BSD License). All experiments have been independently replicated 10 times a laptop equipped with an Intel Core 19.9980HK @ 2.40 GHz with 8 cores (16 threads).
H.2 Benchmarks and Figures
We provide here detailed description of each implemented benchmark and the associated figures. There are two figures associated with each benchmark, showing their average regrets and the size of their datasets throughout the experiment.
In the following, the synthetic benchmarks will be described as functions of a point z in the + 1- dimensional ppatio-tempora domain
More precisely, the point z is explicitly given by (Z1, Zz. t). Also, we will write d' = d + for the sake of brevity.
Rastrigin. The Rastrigin function is d' d'-dimensional, and has the form
f(z) = ad' + Σ 2f a cos (2mzz)
For the numerical evaluation, we set
10. d' and optimized the function on the domain [-4, 4]d". The results are provided in Figure g
Schwefel. The Schwefel function is d' d'ddimensional, and has the form
f(z) 414.988899
sin |2:
For the numerical evaluation, we set d' and we optimized the function on the domain [-500,500]d" The results are provided in Figure 9 This benchmark has also been used to replicate results using the ARD covariance function studied in Appendix E The results provided in Figure 10
Styblinski-Tang
The Syblinski- Tang function is d' d'-dimmesional, and has the form
f(z)
162, + 5zz
1=1
27

---

Schwefel4d
Schwefel4d
1000
GP-JCB
R-GP-UCB TV-CP-UCB
GP-UCB
R-GP-UCH
ABO
ABO
TV-GPAIC3
ET-GP-UCE
ET-GP-UCE
W-DBO
Average Rearet
Dataset
400
200
101
2
2
2
2
2
0
100
200
300
400
500
600
Average Response
Duration
Figure 9: (Left) Average response time and average regrets of the DBO solutions during the opti- mization of the Schwefel synthetic function. (Right) Dataset sizes of the DBO solutions during the optimization of the Schwefel synthetic function.
Schwefel4d ard
Schwefel4d ard
GP-UCB
R-GP-UC
GP-LCB
R-GP-UCB
1000
ABO
ET-GP-UCB
TV-GP-UCB
ABO
ET-GP-UCB
W-DBO
IV-GP-UCB
Average Regget
Datase: Size
10"
0
600
400
103.
20
100
200
400
Average Response Time (s)
Duration (s)
Figure 10: (Left) Average response time and average regrets of the DBO solutions using an ARD SE kernel during the optimization the Schwefel synthetic function. (Right) Dataset sizes of the DBO solutions using an ARD SE. kernel during the optimization of the Schwefel synthetic function.
For the numerical evaluation, we set d' = 4 and we optimized the function on the domain [-5, 55." The results are provided in Figure
Eggholder:
The Eggholder function is 2-dimensional, and has the form
f(z) (22 + 47) sin (v
21 2
+47 sin (Viz1 2/ 47
22 +
For the numerical evaluation, we optimized the function on the domain --12,,12]] The results are provided in Figure [2
Ackley.
The Ackley function is d' d'-dimensional, and has the form
8'
f(z) exp
d'
2
2
exp
d'
cos(c2i)
+ + exp(1).
1=
For the numerical evaluation, we set 20, b- 0.2, 2n, d' = 4 and we optimized the function
on the domain [-32, 32]d', The results are provided in Figure 13 This benchmark has also been used
to replicate our results using the ARD covariance function studied Appendix E The results are provided in Figure [4
28

---

StyblinskiTangee
StyblinskiTang4d
GRUUCB ABO
R-GP-UCE
100
TV-GP-UCB
Datasat Sizp
10a
Average
GP-UCB ABO
TV-GP-UCB
ET-GP-UCB
W-DBO
2
2
100
200
300
40D
500
600
Average Response (s)
Duration
Figure 11: (Left) Average response time and average regrets of the DBO solutions during the optimization of the Styblinski-Tang synthetic function. (Right) Dataset sizes of the DBO solutions during the optimization of the Styblinski-Tang synthetic function.
Eggholder2d
Eggholder2d
600
Average Regret
GP UCB
R GP UCB
Dataset Size
102-
ABO
TV-GP-UCB
ET-GP-UCB
300
101
GP-UCB ABC
R-GP-UCB
TV-GP-UCB W-DBO
200
EI-GP-UCS
2
2
100
200
300
400
500
600
Average Response Time (5)
(s)
Figure 12: (Left) Average response time and average regrets of the DBO solutions during the optimization of the Eggholder synthetic function. (Right) Dataset sizes of the DBO solutions during optimization of the Eggholder synthetic function.
Rosenbrock.
The Rosenbrock function is d' d'-dimensional, and has the form
f(z) 100( 21+1 =22) + (2. 1
For the numerical evaluation, we set d' and we optimized the function on the domain [-1, 1.5]
The results provided in Figure 15
Shekel. The Shekel function 4-dimensional, and has the form
f(z) = Σ
Σ 2g Cj.) +
1=1
For the numerical evaluation, we set
10,
B
1
(1,
3,,,6,,,
TO
4
1
8
6
5
8
6
7
1
8
6
3
1
2
3.6
C
4
1
8
6
5
8
6
7
4
1
8
6
7
3
1
2
3.6
and optimized the function on the domain [0, 10
The results are provided in Figure 16
29

---

Ackley4d
Ackley4d
GP UCB
G UCB TV-GPUCC
ABO
ET-GP-UCB
Average Regret
Dataset Size
101
GP-UCB
R-GP-UCB
ABO
TV-GP-UC3 W-DBO
EI-GP-UCB
20
2
2
100
200
300
400
500
600
Average Response
Duration
Figure 13: (Left) Average response time and average regrets of the DBO solutions during the optimization of the Ackley synthetic function. (Right) Dataset sizes of the DBO solutions during the optimization of the Ackley synthetic function.
Ackley4d ard
Ackley4d ard
GP-UCB
R-GP-UCB
GP-UCB
R-GP-UCB
TV-GP-UCB
ABO
TV-CP-UCa
LT-GP-UCD
W-DBO
ET-CP-UCB
W-DBO
Average Regret
Dataset Size
10?
10a
2
2
100
200
300
400
500
660
Average Response Time
Duration
Figure 14: (Left) Average response time and average regrets of the DBO solutions using an ARD SE Ackley synthetic function. (Right) Dataset sizes the DBO
kernel during the optimization solutions using an ARD SE kernel during the optimization of the Ackley synthetic function.
Hartmann-3. The Hartmann-3 function is 3-dimensional, and has the form
f(z)
ag exp
Σ Anl Pb)2
1=1
For the numerical evaluation, we set a (1.0, 3.2,00
3
10
30
33689
1170
2673
0.1
10
35
4699
4387
7470
A
P. P 10
3
10
30
1091
8732
5547
0.1
10
35
381
5743
8828
and optimized the function the domain [0. 1]3. The results provided in Figure 17
Hartmann-6. The Hartmann-6 function is 6-dimensional, and has the form
4
6 - An - - j=1
f(x) = a; exp -
Σ
30

---

Rosenbrock3d
Rosenbrock3d
GP-UCB
R-GRUCC
GP-UCB
R-GP-UCB
ABO
TV-GP-UCB
ABO
TV-GP.UCB
500
EI-GP-UCH
WDBBO
ET-GP-UCB
W-DBO
Average Regret
Dataset Size
10
100
2a
25
2
100
200
300
430
500
600
Average Response [S]
Duration
Figure 15: (Left) Average response time and average regrets of the DBO solutions during the optimization of the Rosenbrock synthetic function. (Right) Dataset sizes of the DBO solutions during the optimization of the Rosenbrock synthetic function.
Shekel
Shekel
3.00
GP-UCS ABO
R-GP-UCE
TV-GP-UCB
2.75
ET-GP-UCE
W-DBO
2.50
102.
2.25
Everaga
Dataset
2.00
1.75
1.50
GP-UCB
R-GP-UCB TV-GP-UCR
10a
1.25
ET-GP-UCB
20
21
0
100
200
300
100
500
600
Average Response
Duration (s)
Figure 16: (Left) Average response time and average regrets of the DBO solutions during the optimization of the Shekel synthetic function. (Right) Dataset sizes of the DBO solutions during the optimization of the Shekel synthetic function.
For the numerical evaluation, we a (1.0, 3.2,3.0,
10
17
3.50
1.7
1312
1696
5569
124
8283
5886
0.05
10
17
0.1
14
2329
4135
8307
3736
1004
9991
A
P = 10
3
3.5
1.7
10
17
2348
1451
3522
2883
3047
6650
17
0.05
10
0.1
14
4047
8828
8732
5743
1091
381
and we optimized the function on the domain [0, 1]5. The results are provided in Figure 18
Powell.
The Powell function -imeenssionll, and has the form
f(z) Σ (za-ss+ 1024:-2)2 5(245-1 2a;)² + (Z4i-2 2215-1)* 10(24-33 22r)
For the numerical evaluation, we d' and we optimized the function on the domain [-4, 5)d" The results are provided in Figure 19
Temperature. This benchmark comes from the temperature dataset collected from 46 sensors deployed at Intel Research Berkeley. It is famous benchmark, used in other works such as 22 The goal the DBO task is to activate with time. To make the benchmark more interesting, we interpolate the data in space-time. With this interpolation, the algorithms can activate any point space-time, making a 3-dimensional benchmark (2 spatial dimensions for a location in Intel Research Berkeley, I temporal dimension).
sensor with the highest temperature, which will vary
31

---

Hartmann3d
Hartmann3d
GP-UCB
R-GP-UCB
144
ABO
TV-GP-UCE
ET-GP-UCB
W-DBO
Average Hogget
1.0
Dataset Size
10?
0.8
0.6
GP-UCB ABO
0.4
10a
TV-GP-UCB W-DBO
ET-GP-UCB
24
2
2²
100
200
300
400
500
600
Average Response
Duration
Figure 17: (Left) Average response time and average regrets of the DBO solutions during the optimization of the Hartmann-3 synthetic function. (Right) Dataset sizes of the DBO solutions during the optimization of the Hartmann-3 synthetic function.
Hartmann6d
Hartmann6d
1.8
GP-LCB
R B-GP-UCB
TV-GP-UCB
1.6
ET-GP-UCB
W-DBO
Average Regret
1.4-
1.2
10'
Dataset
1.0
GP TV-GP-UCB
ABO
0.8
ET-GP-UCB
W-DB0
0.6
0.2
10
100
200
300
400
500
600
Average Response Time (6)
Duration (s)
Figure 18: (Left) Average response time and average regrets of the DBO solutions during the optimization of the Hartmann-6 synthetic function. (Right) Dataset sizes of the DBO solutions during the optimization of the Hartmann-6 synthetic function.
For the numerical evaluation, used the first day of data. The results provided in Figure 20
WLAN. This benchmark aims at maximizing the throughput of a Wireless Local Area Network (WLAN). 18 moving end-users associated with of fixed nodes and continuously stream a large amount of data. As they move in space, they change the radio environment of the network, which should adapt accordingly that can be tuned for the purpose of reaching the best trade-off between serving all its users and not causing interference for the neighboring nodes.
improve its performance. To do so, each node has power level
The performance of the network is computed as the sum of the Shannon capacities for each pair of node and associated end-users. The Shannon capacity 34 sets theoretical upper bound the throughput wireless communication. We denote it (ii.j), we express it in bits per second (bps). Signal-to-Interferenne It depends on S.j the Signal-to-Ine Interference plus Noise Ratio (SINR) of the communication between
node i and end-user j. as well as on W, the bandwidth of the radio channel (in Hz):
Czy t) = W log (1 + (1y(w, t))-
Then, the objective function is
f(x, t) Σ Σ Ci((,tt) 1=1 JEN.
NG the end-users associated with node i.
32

---

Powell4d
Powell4d
17500
GRAUCB
R.GPACC
GRUCB EI-GPUUCE
R-GP-UCE
15000
ET-GP-UCE
W-DB0
TVVERUCR
TV-GP-UCB
WADBO
Average Regret
12500
Dataset Size
10000
10?
7500
5000
2500
22
22
25
2
100
200
300
400
500
600
Average Response Time (s)
Duration (s)
Figure 19: (Left) Average response time and average regrets of the DBO solutions during the optimization of the Powell synthetic function. (Right) Dataset sizes of the DBO solutions during the optimization of the Powell synthetic function.
Temperature3d
Temperature3d
2.0
GPACB
R.GP-UCH
GP-UCB
R-GP-UC3
ABO
TV-GP-UCB
ABO
TV:GP-UCB
ET-GP-UCB
W-DBC
EI-GP-UCH
W-D80
Average Rogres
166
Dotaset Size
10a
1.4
1.0
0.8
101
0.6
20
2
22
0
100
200
300
400
500
660
Average Resporse Time (s)
(s)
Figure 20: (Left) Average response time and average regrets of the DBO solutions during the Tem- perature real- world experiment. (Right) Dataset sizes of the DBO solutions during the Temperature real-world experiment.
For the numerical evaluation, we optimized the power levels in the domain [10 102.514 For this experiment, the DBO solutions were evaluated with a Matérn-57/ 5/2 for the spatial covariance function and Matérn-1/2 for the temporal covariance function. The results are provided in Figure 21
H.3 Discussion on Empirical Performance
In this section, we discuss the performance achieved by all the DBO solutions on the benchmarks introduced in the previous section.
This baseline, which does take into account temporal correlations, obtains surpris- ingly good performance in this experimental setting (continuous time, hyperparameters estimated on the fly). Its simple behavior (i.e., keep all observations in the dataset until the end of the experiment) hampers its response time, but this drawback is balanced by the fact that it has only three parameters estimate with MLE (i.e., X, ls, o²). Overall, is dominated by R-GP-UCD TV-GP-UCB and W-DBO, but behaves surprisingly well against ABO and ET-GP-UCB (see Figure 5)
GP-UCB.
ABO. ABO performs poorly in this experimental setting. We explain this poor performance by the fact that the hyperparameters (including the spatial and temporal lengthscales) have to be estimated on the fly. Since ABO can decide to postpone its next query to the near future (a fraction of the temporal lengthscale lT away), overestimating lT may cause ABO to wait for a long time before querying f again. This interpretation is supported by the fact that the functions ABO performs the poorest on are the ones with the largest temporal lengthscales lr, e.g., Rosenbrock (see Figure 15) and Powell (see Figure 19 Conversely, ABO obtains competitive performance on functions with
33

---

WLAN5d
WLAN5d
25
GP-UCB
R-GP-UCB TV-GP-UCE
GP-UCB
R-GP-UCE
ABO
TV-GP.UCB
ET-GP-UCE
W-DBO
20-
Aver age Regret
Datasat Sire
102-
104
2-
2
2
2
0
100
200
300
400
500
60C
Average Response
Duration
Figure 21: (Left) Average response time and average regrets of the DBO solutions during the WLAN real-world experiment. (Right) Dataset sizes of the DBO solutions during the WLAN real-world experiment.
smaller temporal lengthscales, c.g. Eggholder (see Figure 12 or Shekel (see Figure 16). These results highlight the lack of robustmess of ABO.
ET-GP-UCB. Like GP-UCB, ET-GP-UCB does not take into account temporal correlations, and deals with stale observations by resctting its dataset each time a condition is met. Our experimental setting exposes the lack of robustness of ET-GP-UCH since its performance is quite poor most benchmarks. This is mainly due to the fact that, because the hyperparameters (including the observational noise level o²) inferred the fly, the MLE explains the variance in the observations with an increasingly large observational noise level o² as time goes by. However, by construction of ET-GP-UCB, the greater o², the less the dataset will be reset. As consequence. on some benchmarks, the event is never triggered (or not triggered enough) and the performance of ET-GP-UCB is close (sometimes worse than) the performance of GP-UCB. For examples, refer Hartmann6d (see Figure (8). Shekel (see Figure [6) or Ackley (see Figure 13.. Some other times, the variance in the observations cannot be explained by an increasingly large observational noise. In these cases, the triggering occurs properly and ET-GP-UCB obtains competitive performance, e.g., with Powell (see Figure
R-GP-UCB. R-GP- UCB deals with stale data by resetting its dataset (like ET-GP-UCB). It indi- rectly accounts for temporal correlations by estimating an hyperparameter c, and the reset is triggered each time the dataset size exceeds N(c) (given in (21). More often than not, its performance is better than GP-UCB, because stale data is frequently removed from the dataset. As a consequence, R-GP-UCB has the lowest average response time of all the DBO solutions. Overall, because of its low response time, R-GP-UCB obtains very good performance on some benchmarks, e.g., Eggholder (see Figure 12p. Hartmann3d (see Figure 17 or Temperature (see Figure 20
TV-GP-UCB. TV-GP-UCB directly accounts for temporal correlations by computing a specific covariance function controlled by an hyperparameter e. Although the temporal covariance function is based on distance between indices instead of a distance between points in time, the DBO solution turns out to be quite robust in our experimental setting. However, its response time is hampered by the irrelevant observations that kept in the dataset. Because of them, TV-GP-UCB has one of the largest response time on many benchmarks, e.g., Rastrigin (see Figure 88. Schwefel (see Figure 9 or Shekel (see Figure 16). Nevertheless, its average performance is significantly better than the other state-of-the-art DBO solutions.
W-DBO. Because of its ability to measure the relevancy of its observations and remove irrelevant observations, W-DBO achieves simultaneously good predictive performance and low response time. Depending on the benchmark, its dataset size follows different patterns. When the objective function evolves smoothly, e.g., Powell (see Figure 19) or Rosenbrock (see Figure 15). W-DBO behaves roughly like GP-UCB and TV-GP-UCB and keeps most of its observations in its dataset (although it manages to identify and delete some irrelevant observations). When the objective function's
34

---

0.8
0.6
Space (X]
0.4
Observations
ArgMex
0.0
0.2
0.4
0.5
0.3
Time (t)
Figure 22: Snapshot from one of the videos showing the optimization conducted by W-DBO. The normalized temporal dimension is shown the x-axis and the normalized spatial dimension is shown on the y-axis. The observations that are in the dataset are depicted as red dots. while the deleted observations depicted as black crosses. The maximal arguments farg max,es far t), t T are depicted with cyan curve. The predictions of W-DBO are shown with a contour plot. Finally, the present time is depicted black vertical line labelled fg-
variations are more pronounced, the dataset size W-DBO experiences sudden drops, as can be seen with Ackley (see Figure 13 Shekel (see Figure [6] Temperature (see Figure 20). This suggests that W-DBO is also able to "reset" its dataset, although in a more refined way as it is able to keep the few observations still relevant for future predictions. Thanks its ability adapt in very different contexts, W-DBO outperforms state-of-the-art DBO solutions by a comfortable margin. This performance gap can be seen its average performance across all benchmarks (see Figure (5) but also on most of the benchmarks themselves, c.g., Schwefel (see Figure 9 Ackley (see Figure Shekel (see Figure [6]. Hartmann-6 (see Figure 18) or Powell (see Figure
H.4 Animated Visualizations
In this section, we describe and discuss the two animated visualizations provided as supplementary material for the paper. These videos show W-DBO optimizing two 2-dimensional synthetic functions. They depict W-DBO's predictions, collected observations and deleted observations into the spatio- temporal domains of the functions.
One of the videos depict the optimization of the Six-Hump Camel function on the domain [-2, 2)2. The SHC function is
I1 x² + rt (-4 412) 12.
SHC(#, 1)
2.1x²
To study how W-DBO reacts to sudden changes in the objective function, the other video depict the optimization of the piecewise function
SHC(x,,)) t) SHC(t,
ift -
f(x,t)
otherwise.
snapshot from the latter can be found Figure 22 illustrates that the benefits brought by
W-DBO are substantial, since the algorithm is able track max (x, over the time while
'The video is accessible at hatpps///abbaduu...m
Hump_ Hump_Camel__ 25.
240 mp4
35

---

simultaneously deleting a significant portion of collected observations. Indeed, many observations are deemed irrelevant, either because (i) they have become stale (there are only a few observations collected at the start of the experiment that have been kept in the dataset) or because (ii) they are redundant with observations that are already in the dataset (many observations are located near the maximal argument, and many of them are deleted soon after being collected).
I Limitations
For the sake of completeness, we explicitly discuss the limitations of W-DBO this appendix. Four limitations were identified:
As for any BO algorithm, W-DBO exploits a GP as surrogate model (see Assumption 3.1 If the objective function f cannot be properly approximated by GP, expect the perior- mance of W-DBO to decline.
As for any BO algorithm, W-DBO conducts GP inference, which causes it to manipulate inverses Gram matrices that scale with the dataset size. Although the main motivation of introducing W-DBO is to reduce the dataset size, the cubic complexity of matrix inversion algorithms can still constitute limitation if too many observations are kept in the dataset.
We also introduce structure for spatio-temporal correlations with Assumption 3.2 Although less restrictive than the one enforced by 21 [22, equivalent to the one in [20] and partially
relaxed in Appendix E this is still a limitation since we expect the performance of W-DBO to worsen the objective function does not meet this assumption.
Finally, W-DBO is not exempt from the effects the sampling frequency. In fact, as for any DBO algorithm, the performance of W-DBO will drop if the function varies too much between observations. As an example, if f evolves so rapidly that two successive observations become basically independent, W-DBO will not be able to infer anything meaningful about the objective function.
NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract states that the paper introduces a criterion that measures the relevancy of the observations and DBO algorithm that exploits the criterion to improve upon state-of-the-arr, These contributions can be found in Sections aaddd
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
The limitations are explicitly discussed in Appendix I and the approximation error of our proposed upper bounds is discussed at length in Appendix C
Justification:
Guidelines:
36

---

The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g. if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: The assumptions explicitly given in the main text (see Section 3), and the
proofs for each result are provided Appendices A. Q and D
Guidelines:
The answer NA means that the paper does include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification; The experimental details are discussed briefly in the main paper (see Section 5) and discussed at length in Appendix H.1 Furthermore, the benchmarks are detailed in Appendix HH.2
Guidelines:
The answer NA means that the paper does not include experiments.
37

---

If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g.. with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors welcome describe the particular way they provide for reproducibilityy In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: A GitHub repository with minimal working example and detailed docu- mentation is provided along with the camera-reate ready version.f
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https //nips.cc/public/eee.
www.ded/eeooeeuummemccoo ssi ssionPolicy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
38

---

Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification: As previously mentioned, all experimental details are discussed in Appendix H A sensitivity study and discussion about the hyperparameter of W-DBO are also provided in Section 3.1 and Appendix respectively.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Standard errors are depicted graphically in Section S and Appendix H.2
Guidelines:
The answer NA means that the paper does include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (c.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of exccution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
The computing resources described in Appendix H.I
Guidelines:
The answer NA means that the paper does not include experiments.
39

---

The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
Justification:
The research conducted in this paper does not violate the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there a special consid- eration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
This paper proposes DBO algorithm that can be applied to many different applied research areas (e.g., robotics, computer networks, biology). These impacts are briefly mentioned in Section 6
Justification:
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML.)
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
40

---

Answer:
[NA]
Justification:
In our opinion, the paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: The paper uses BOTorch and PyBindl1, which are explicitly mentioned in Appendix H.T and whose licenses allow free exploitation for scientific research purposes.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URI
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, pwpprswithhaeecoo thcode. com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-pan se-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors encouraged to reach out the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not introduce new assets yet. Upon acceptance, a GitHub repository with a detailed documentation will be provided.
Guidelines:
The answer NA means that the paper does release assets.
Researchers should communicate the details the dataset/code/moodll as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
41

---

14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The research does involve crowdsourcing nor human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/teview based the requirements of your country institution) were obtained?
Answer: [NA]
Justification:
The research does involve crowdsourcing nor human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
42


Reviews:
[{"review_id": "FZsvDvJS6L", "reviewer": "Reviewer_mGHg", "summary": "This paper addresses the challenge of optimizing a time-varying black-box function using GP-based Dynamic Bayesian Optimization (DBO). Unlike traditional Bayesian Optimization (BO), DBO seeks to handle dynamic functions where the optimum changes over time by incorporating time in the GP model covariance function. With the goal of expediting DBO, this paper proposes a strategy to remove ‘irrelevant’ data from the GP model, specifically by removing points that minimally affect the distribution over future predictions. To this end, the authors introduce a Wasserstein distance-based criterion (and associated bounds/approximations) and propose the W-DBO algorithm. This algorithm dynamically removes irrelevant observations, maintaining good predictive performance and high sampling frequency. Numerical experiments demonstrate W-DBO's superiority over state-of-the-art methods.", "questions": "- Is there ever a downside with the ‘greedy’ of removal of points? E.g., removing the first selected point prevents otherwise being able to remove two points?\n- Can the algorithm remove observations that become relevant again in the future? This is difficult to track, since $\\mathcal{GP}_\\mathcal{D}$ is updated with the removed sample in Algorithm 1, meaning there is no tracking of the original full model with no data removed.", "limitations": "Yes, discussed in appendices.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The introduction of a Wasserstein distance-based criterion to measure the relevancy of observations is an intuitive and interesting that effectively addresses the introduced challenge of minimally altering the GP predictions.\n- The proposed W-DBO algorithm appears to perform well on the extensive numerical experiments, demonstrating the importance of removing irrelevant observations. The algorithm remains computationally efficient, avoiding the prohibitive growth of dataset size over time.\n- The strategy is relatively generalizable; a Wasserstein distance-based approach can be integrated with any BO algorithm to identify candidates for removal from the dataset, enhancing applicability of this paper across various domains.", "weaknesses": "- The motivation for this challenge largely is based on GP-based Bayesian optimization, which might not always be applicable for this problem setting. Specifically, GP-based BO is most relevant when samples are expensive (and are limited as a result). When samples are cheap to obtain, derivative-free optimization methods such as evolutionary strategies are more relevant, but these are not compared. \n- While the numerical experiments are extensive, they are limited to synthetic examples—most of which are not time-dependent (Appendix G.2). The inclusion of some real-world examples and or dynamic optimization case studies could strengthen the motivation and highlight applicability of the proposed method.\n- The scalability of the W-DBO algorithm itself (e.g., checking the proposed metric for all points) in high-dimensional spaces or with very large datasets is not fully addressed, which could be a limitation in some applications.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer mGHg,\n\nThank you for the detailed review. We are discussing below the weaknesses and questions you have raised. Also, please make sure to read the global response as we discuss some of your questions there.\n\n**The scalability of W-DBO (e.g., in high-dimensional spaces and/or with large datasets) is not discussed.**\n\nThe Wasserstein distance is applied in the output space, which is one-dimensional regardless of the dimensionality of the input space, *i.e.*, the objective function domain.\n\nOur framework, which measures the observations relevancy and removes the irrelevant observations, can be seen as a simple post-processing stage running at the end of each iteration. As a consequence, it can be used in conjunction with any BO algorithm. The W-DBO algorithm presented in this paper exploits vanilla GP-UCB, and as such, will struggle in high-dimensional input spaces, precisely because GP-UCB struggles in high-dimensional input spaces. However, our framework could be paired with any state-of-the-art high-dimensional BO algorithm (e.g., [1]) and show good performance when optimizing high-dimensional dynamic black-boxes. This is an interesting future work, thank you for suggesting it.\n\nRegarding very large dataset sizes, W-DBO suffers from the same limitations as any BO algorithm. This is because it also manipulates inverses of Gram matrices, which scale with the dataset size. That is precisely why we argue that removing irrelevant observations gives W-DBO a decisive advantage upon other DBO solutions. Nevertheless, upon acceptance, we will discuss this limitation explicitly in the Appendix H of the camera-ready version.\n\n**Is there ever a downside with the ‘greedy’ of removal of points?**\n\nOur framework could easily be extended to the power set $2^\\mathcal{D}$ of the dataset $\\mathcal{D}$, in order to measure the relevancy of a subset of observations, instead of only quantifying the relevancy of a single observation.\n\nThe greedy approach causes W-DBO to overestimate the impact of the observation removals on the GP posterior. To rapidly get an intuition explaining this, let us describe a simple example. Consider a dataset that contains an irrelevant pair of observations $\\mathcal{P} = \\\\{\\mathbf x_i, \\mathbf x_j\\\\}$, in the sense that $W_2\\left(\\mathcal{GP}\\_\\mathcal{D}, \\mathcal{GP}\\_{\\mathcal{D} \\setminus \\mathcal{P}}\\right) = 0$. Furthermore, let us assume that $W_2\\left(\\mathcal{GP}\\_\\mathcal{D}, \\mathcal{GP}\\_{\\mathcal{D} \\setminus \\\\{\\mathbf x_i\\\\}}\\right) = W_2\\left(\\mathcal{GP}\\_{\\mathcal{D} \\setminus \\\\{\\mathbf x_i\\\\}}, \\mathcal{GP}\\_{\\mathcal{D} \\setminus \\mathcal{P}}\\right) = \\epsilon > 0$. If W-DBO were able to directly capture the irrelevancy of the couple $\\mathcal{P}$, it would have removed it without consuming any of its removal budget. However, the greedy removal procedure as described in Algorithm 1 causes W-DBO to remove first $\\mathbf x_i$ and next $\\mathbf x_j$. Consequently, it will consume $2\\epsilon$ from its removal budget and, in that sense, it will overestimate the impact of removing $\\mathbf x_i$ and $\\mathbf x_j$.\n\nAlthough working with the power set $2^\\mathcal{D}$ of the dataset $\\mathcal{D}$ is clearly more advantageous to avoid such budget depletion caused by greedy removals, the combinatorial explosion of $|2^\\mathcal{D}|$ makes this strategy prohibitive in practice, even for moderately-sized datasets. That is why it is not advertised in the paper. Upon acceptance, we will include this remark in Section 4 of the camera-ready version. Thank you for pointing it out.\n\n**Can the algorithm remove observations that become relevant again in the future?**\n\nFor usual stationary, decreasing kernels (e.g., Squared-Exponential or Matérn) and because of Assumption 3.2, the covariance between a function value at a time $t$, *i.e.*, $f(\\mathbf x, t)$, and a function value at a future time $t’$, *i.e.*, $f(\\mathbf x’, t’)$ with $t < t’$, can only decrease as $t’$ gets further and further away from $t$. Consequently, the relevancy of the observing $f(\\mathbf x, t)$ will only decrease as time goes by.\n\nAn interesting case that is not discussed in the paper, is when the time kernel $k_T$ is a periodic kernel, e.g., $k_T(t, t’) = \\exp\\left(-\\frac{2 \\sin^2\\left(\\pi |t - t’| / p\\right)}{l^2}\\right)$. The objective function is thus periodic along its time dimension. In that case, the concept of stale observation vanishes completely since the observation of $f(\\mathbf x, t)$ will always be useful to predict $f(\\mathbf x, t’)$, even when $t’ \\gg t$. The Wasserstein distance between the two posterior GPs with this time-periodic kernel can be shown to diverge to $+\\infty$, and so does our approximation, as it indeed should.\n\n**References**\n\n[1] Bardou, A., Thiran, P., & Begin, T. Relaxing the Additivity Constraints in Decentralized No-Regret High-Dimensional Bayesian Optimization. In The Twelfth International Conference on Learning Representations. 2024.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "XygcU8OSAi", "reviewer": "Reviewer_JgwK", "summary": "The paper proposes a new algorithm for dynamic Bayesian optimization (DBO). To develop this new algorithm, the authors first derive a Wasserstein distance-based criterion that is a way of measuring how relevant a given collected data point is during optimization. Since dynamic functions change over time, each observation collected during DBO becomes less relevant over time, until eventually the computational cost of continuing to consider the observation outweighs any benefit it provides. When this becomes the case it therefore makes sense to remove these older observations from the dataset, especially since we care a lot about reducing computational cost since having high sampling frequency is especially important in DBO. In order to define a principled way to remove these irrelevant observations during DBO, the authors use their Wasserstein distance-based criterion to measure the relevancy of collected data points and decide whether to remove them from the dataset during the course of optimization. This strategy leads to the author’s novel DBO algorithm (W-DBO), which the authors show to perform better than relevant baselines with a convincing set of experimental results.", "questions": "The authors state that computational biology is a field that makes “heavy use of DBO”. I am curious what scenarios in comp bio require the use of DBO rather than just traditional BO? In what scenarios do we have biological black-box functions that are changing over time? It might be useful to have some additional citations and/or discussion of this in the paper.", "limitations": "Yes.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Originality: The author’s proposed Wasserstein distance-based criterion and resultand W-DBO algorithm is clearly a novel approach to solve a relevant problem in DBO. \nQuality: The paper is very well-written. Additionally, the figures and tables are all of good quality - they are both easy to parse and do a nice job of displaying relevant results.  \nClarity: The paper is clear and easy to follow from start to finish. The figures and tables are clear and easy to read. The paper is also clearly motivated and it’s easy to understand what the authors did and why. \nSignificance: It is obvious to me that the problem the authors seek to solve here (DBO algorithms keeping around increasingly irrelevant data points despite the need to maintain high sample efficiency over time) is important and relevant to the community. The author’s W-DBO algorithm provides a reasonable and intuitive solution to this problem. \nConvincing results: The experimental results do indeed show that W-DBO outperforms other methods on a large number of BO tasks.", "weaknesses": "Typo (Line 305): “to the best of your knowledge” should instead be “to the best of our knowledge”. \n\nOne ablation I would’ve like to see in this paper: The paper makes it clear that leveraging quantification of relevancy of observations can be used to improve DBO performance with their W-DBO algorithm. However, it is less clear to me from the experimental results that using the author’s proposed Wasserstein distance-based criterion method for relevancy quantification is necessarily better than another strategy for quantifying relevancy. While the Wasserstein distance-based criterion method the authors propose is intuitive and clearly works well in this setting, it would be interesting to see a direct comparison against other simpler (or even ad-hoc) methods for attempting to quantify how relevant a given observation is. For example, rather than using the Wasserstein distance-based criterion, one could just assume some constant rate of decline in quantitative relevancy over time. I do not think that this additional experiment is necessary for this paper to be accepted, but I do think it would strengthen the paper by showing the importance of using the author's proposed Wasserstein distance-based criterion method specifically for relevancy quantification in this setting.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer JgwK,\n\nThank you for the detailed review. We are discussing below the weaknesses and questions you have raised. Also, please make sure to read the global response as we discuss some of your questions there.\n\n**Direct comparison against other simpler (or even ad-hoc) methods for relevancy quantification (e.g., constant rate of decline over time).**\n\nTV-GP-UCB assumes a constant decreasing rate of the covariance function over time. Neglecting spatial relevancy, this precisely translates into a constant decreasing rate of the relevancy of an observation over time. Our approach (W-DBO) differs from TV-GP-UCB in three ways:\n\n(i) W-DBO can use an arbitrary time kernel instead of a fixed kernel for TV-GP-UCB (see [1]),\n\n(ii) W-DBO has a more sophisticated quantification of relevancy (*i.e.*, the Wasserstein distance) that involves the relevancy both in the time dimension (staleness) and in the spatial domain,\n\n(iii) W-DBO is able to remove observations when they are deemed irrelevant.\n\nTogether, (i), (ii) and (iii) explain the difference in the performance of TV-GP-UCB and the performance of W-DBO.\n\nOne could use the very same core idea of our framework (a distance function on stochastic processes) and come up with other definitions of relevancy. As an example, one could define the relevancy of an observation using the KL-divergence on the two GP posteriors. However, this would require a very different analysis to find an approximation of this distance, which would lead more to a new paper than to an ablation study.\n\n**The authors state that computational biology is a field that makes “heavy use of DBO”.**\n\nIn Section 1, we cite computational biology as an application of BO, not of DBO. However, it is true that in our concluding remarks of Section 6, we cite computational biology as a field that could apply DBO. This is a typo that will be removed from the camera-ready version upon acceptance. Thank you for pointing it out.\n\n**References**\n\n[1] Ilija Bogunovic, Jonathan Scarlett, and Volkan Cevher. Time-varying gaussian process bandit optimization. In Artificial Intelligence and Statistics, pages 314–323. PMLR, 2016.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "hkaTmzCwbQ", "reviewer": "Reviewer_iNzn", "summary": "The authors propose W-DBO, a statistical distance-based criterion for removing \"stale\" observations in the dynamic BO setting. Observations are removed based on their impact on the GP globally, as measured by an approximate integrated Wasserstein distance, for which the authors prove the approximation quality. Results show that w-DBO outperforms the competition in terms of regret over time.", "questions": "Questions are apparent from the \"Weaknesses\" section.", "limitations": "Limitations have been addressed.", "rating": 3, "confidence": 5, "soundness": 1, "presentation": 2, "contribution": 2, "strengths": "__Relevant problem__: The optimization of time-varying functions appears highly relevant, and seems underexplored.\n\n__Intuitive solution__: Measuring staleness by the impact on the GP globally is a good idea.\n\n__Theory__: There is a good amount of theory included on the proposed algorithm.", "weaknesses": "Unfortunately, I have a substantial concern about the evaluation, which I believe requires the authors to (ideally) release code of their experimental setup and unfortunately, re-run experiments with a more conventional setup.\n\n__Benchmarking__: The authors do not appear to use an ARD kernel in their work. As far as I can tell, the method is not restricted to a single lengthscale for all dimensions, so is there a reason for this unconventional choice? \n\nNotably, using a non-ARD kernel is an issue of fairness in benchmarking. Since almost all test functions (Rastrigin, Schwefel, Shekel, Ackley, Styblinski-Tang and to a lesser degree Rosenbrock), that are chosen are symmetric, this choice is rather convenient. Specifically, as soon as one lengthscale is estimated correctly, all of them will be, which is _very_ beneficial for performance. For one, ABO estimates all of the lengthscales, and will as such be as a large disadvantage. The authors mention that standard GP-UCB is non-ARD as well (App. G3), which could very well be the explanation for its good performance on the symmetric functions. To emphasize this even more, the obvious non-symmetric function (Hartmann-3) is one one where W-DBO and GP-UCB perform worse. (Sidenote: Hartmann-6 is fairly uniform in its lengthscales.)\n\nIn short: Can the authors please run their experiments with an ARD kernel on the symmetric test functions? Right now, the non-ARD provides a alrge advantage, specifically in relation to the algorithms that run ARD (ABO at least). I do not believe running a non-ARD kernel is warranted, since it is not a design choice that is relevant to the dynamic setting and provides an outsized and unrealistic advantage.\n\n__Relevance of \"stale\" observations__: It is not apparent to me that \"stale\" observations are a large problem if there is a time-varying component to the kernel which naturally decreases its relevance with time. As such, the algorithm seems to strictly target high-throughput settings due to the computational saving, which is more restrictive setting than dynamic BO generally.\n\n__Evidence for removing stale observations__: While the criterion makes sense, I fail to see evidence that the removed observations are indeed \"stale\". What is good evidence for this? Simply plotting downstream regret performance (especially when the benchmarking has other flaws) is not a good metric for staleness.\n\nAn Illustrative example of which observations are being removed are a must. However, I also encourage the authors to assess what a good metric for staleness is, and present that in their results.\n\n__Benchmarking, part 2__: The authors propose a dynamic BO algorithm, but the synthetic benchmarks appear to all be non-dynamic (there is no time-dependence in any of the benchmarks in App. G2). I would have expected these to have a time-varying component in order for the algorithm to be relevant. If the benchmarks indeed do not have a time-varying component, shouldn't standard GP-UCB be the gold standard for performance?\n\n__Benchmarking, part 3__: Please add standard errors to the tables.\n\n__Content in Appendices__: There is a lot of important content that has been moved to the supplementary material, most importantly the variable definitions in 4.1 (App. A & B), intuition for the Wasserstein ratio in 3.2 (App. E) . I believe the authors should try and make room for these important pieces in the main paper, possibly by shortening the background section or moving some of it to SM.\n\n__Normalization by empty GP__: The idea is novel and appears intuitive, but it is not thoroughly explained insofar as why the metric is appropriate. I recommend the authors to elaborate and add a 1D illustrative example, possibly compared to simply using the Wasserstein distance.\n\nMinor: \n- \"BO inference\" is used numerous times in place of \"GP inference', which is the more precise description of what is occuring.\n- L160: \"uses\" -> \"used\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer iNzn,\n\nThank you for the detailed review. We are discussing below the weaknesses and questions you have raised. Also, please make sure to read the global response as we discuss some of your questions there.\n\n**On the usage of an ARD kernel.**\n\nWe make no mention of setting an ARD kernel for ABO, nor for any other DBO algorithms. Appendix G.1 lists the kernels used for each solution. Let us detail precisely the number of hyperparameters that each solution has to infer:\n\n_GP-UCB and ET-GP-UCB_.\n* Spatial cov.: Matérn-5/2\n* Temporal cov.: N/A\n* Number of parameters: 3 (scale $\\lambda$, spatial lengthscale $l_S$, noise level $\\sigma^2_0$)\n\n_R-GP-UCB and TV-GP-UCB_.\n* Spatial cov.: Matérn-5/2\n* Temporal cov.: fixed kernel with decay parameter $\\epsilon$\n* Number of parameters: 4 (scale $\\lambda$, spatial lengthscale $l_S$, decay $\\epsilon$, noise level $\\sigma^2_0$)\n\n_ABO and W-DBO_. \n* Spatial cov.: Matérn-5/2\n* Temporal cov.: Matérn-3/2\n* Number of parameters: 4 (scale $\\lambda$, spatial lengthscale $l_S$, temporal lengthscale $l_T$, noise level $\\sigma^2_0$)\n\nAlthough we acknowledge the usefulness of anisotropic kernels in practice, we would like to address (i) the issue of “conventionality” of ARD kernels in benchmarking, and (ii) the performance of W-DBO with an ARD kernel.\n\n(i) Among the DBO papers (*i.e.*, [1, 2, 3]), the vast majority of the benchmarking was done with a non-ARD kernel. More specifically, [1] explicitly uses a non-ARD kernel, [2] uses an ARD kernel on 1 experiment out of 7 (cf. Tables 4-10 in [2]), and [3] makes no mention of an ARD kernel. This makes us consider that using a non-ARD kernel for benchmarking is fairly standard, at least in the DBO literature. We do not understand therefore why the comprehensive benchmarking would be unconventional and/or unfair towards some solutions. Upon acceptance, we propose to explicitly list the kernel hyperparameters for each DBO solution in Appendix G.1 to prevent any confusion about our benchmarking.\n\n(ii) Let us discuss how our framework can be applied to anisotropic kernels. Most of our analysis remains unchanged, but we need additional convolution formulas because the ones provided in Tables 3 and 4 only hold under Assumption 3.2. We now illustrate this point by providing additional results with the anisotropic SE kernel. Consider a $d \\times d$ diagonal matrix $\\Sigma = \\text{diag}\\left(l_1, \\cdots, l_d\\right)$ that gathers the $d$ lengthscales in its diagonal. It can be shown that the anisotropic convolution formula is $(k_S * k_S)(\\mathbf x) = \\pi^{d / 2} \\det\\left(\\Sigma\\right) e^{-\\mathbf x^\\top \\Sigma^{-2} \\mathbf x / 4}$. Observe that the formula reduces to the isotropic formula described in Table 3 when $l_1 = \\cdots = l_d = l_S$. Additional experiments with the ARD SE kernel are listed on the PDF, which had to be uploaded in the global rebuttal (where they are discussed too).\n\n**I fail to see that the removed observations are stale. Provide an illustrative example of which observations are removed and a metric for staleness.**\n\nThe concept of staleness is used in the title, as well as in Sections 1 and 2 to introduce our work because it is a term commonly used in the literature, but we replace it with the broader concept of “relevancy” from Section 3 onwards. Relevancy can be temporal (an observation that is too distant in time to be relevant for GP inference, a.k.a staleness) but also spatial (an observation that is too close from another to bring substantial information for GP inference). In other words, staleness coincides with temporal relevancy, whereas W-DBO accounts for both temporal and spatial relevancy. We provide a metric for the relevancy of an observation $\\mathbf x_i$, namely the Wasserstein distance between the posterior GP conditioned on the entire dataset $\\mathcal{D}$ and the posterior GP conditioned on $\\tilde{\\mathcal{D}} = \\mathcal{D} \\setminus \\\\{\\mathbf x_i\\\\}$, defined on the cartesian product of the spatial domain $\\mathcal{S}$ and the future $[t_0, +\\infty]$, with $t_0$ the present time.\n\nWe use this metric because it is intuitive and it captures properly the concept of observation relevancy (e.g., it is 0 when removing an observation does not change the GP posterior in any way, and it can be arbitrarily large otherwise, depending on how the removal of an observation affects the posterior). We remove observations according to a normalized version of this metric, to ensure that the removed observations are indeed “stale” or, more broadly speaking, “irrelevant”.\n\nWe also provide an example illustrating which observations are being removed. Animated visualizations are provided in the supplementary material, and are discussed in Appendix G.4.\n\n**Please add standard errors to the tables.**\n\nDue to lack of space, we did not report the standard errors directly in Table 2. However, they are provided graphically whenever applicable, namely in Figures 1-3 and 6-17.\n\nTo further address this weakness while respecting the NeurIPS template, we propose to replicate Table 2 in a dedicated appendix of the camera-ready version, where we will also provide the standard errors.\n\n**There is a lot of important content that has been moved to the SM.**\n\nUpon acceptance, we will use the additional page allowed to move some important content from the SM to the main paper, and by possibly also reducing the Background section to make room for this material.\n\n**Typos.**\n\nThank you for pointing out these typos. Upon acceptance, we will fix them in the camera-ready version.\n\n**References**\n\n[1] I. Bogunovic, J. Scarlett, and V. Cevher. Time-varying gaussian process bandit optimization. In AISTATS, pages 314–323. PMLR, 2016.\n\n[2] P. Brunzema, A. von Rohr, F. Solowjow, and S. Trimpe. Event-triggered time-varying bayesian optimization. arXiv preprint arXiv:2208.10790, 2022.\n\n[3] Nyikosa, F. M., Osborne, M. A., & Roberts, S. J. (2018). Bayesian optimization for dynamic problems. arXiv preprint arXiv:1803.03432.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "EVJgB2sMV5", "reviewer": "Reviewer_NVX3", "summary": "The paper proposes a novel algorithm for dynamic Bayesian optimization (DBO). DBO defers from traditional BayesOpt as the black-box function to be optimized is changing in time, and the goal of the optimization procedure is to keep track of the optimum across a continuous time index. The continuous time component adds an interesting dimension to the problem: the black-box sampling frequency becomes important to be able to query the function frequently enough to keep track of the optimum. This means there is a necessity for the algorithm and acquisition function solve to be as efficient as possible.\n\nOne important important factor in the acquisition function solve time is the training of the GP surrogate model. In particular, the training is of order $O(N^3)$ on the number of training data-points. In a dynamical system where the function evolves with time, certain training-points can become redundant, meaning we can make the optimization more efficient if we discard them. The paper proposes a method to choose which points are worth discarding.\n\nThe main intuition behind the method, is that a point which does not change the *future* predictive distribution of the GP (with respect to time) is a point that is no longer needed in the data-set. Therefore the authors propose a criterion that looks at the predictive distribution under the full data-set $D$ and measures the 2-Wasserstein distance to a GP trained on $D \\backslash (x_i, y_i)$ where $(x_i, y_i)$ is training point that could potentially be removed. To make the measure interpretable (in terms of the magnitude of the criterion), we can normalize by looking at the 2-Wasserstein distance between the full data-set and the prior GP, giving the final criterion. Unfortunately, actually calculating the criterion is too computationally expensive, especially when trying to reduce costs to have a high sampling frequency. The paper then shows that both 2-Wasserstein distances can be upper-bounded, and proposes using the upper-bounds to replace the actual distance, creating an approximate criterion. \n\nIn the appendix, the magnitude of the error of the approximation is investigated in practice. While it is shown to be non-negligible, the division in the criterion helps balance out the error in part, and most importantly, the empirical results of using the criterion are strong. Once the criterion is established, the paper introduces the algorithm: before every acquisition function solve, points are removed from the data-set in a greedy manner, i.e., the point with the least relevancy is compared against a budget, and if the budget allows it, then the points is removed and the budget is decreased, this is repeated until the budget is small enough. The budget is defined by the amount of relative error allowed by removing data-points. \n\nNumerical results show that a trade-off between removing enough information to be computationally tractable, while keeping enough data to identify the optimum is required, and through empirics the authors are able to provide a recommended hyper-parameter value. Then the algorithms are tested on 10 synthetic functions and 2 real-world experiments. The proposed algorithm is the strongest in 10 of the 12 benchmarks.", "questions": "See weaknesses", "limitations": "Limitations are discussed at length in Appendix C and H.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- (Originality) The paper proposes a novel algorithm in an understudied area of Bayesian optimization. The paper provides a strong theoretical backing of the method and a good empirical study. Some important literature connections seem to have been missed though. \n\n- (Quality) All claims in the paper seem backed-up and well justified. Details of the method are well reported, and important parts of the algorithm, such as sensitivity to hyper-parameters and approximation error are investigated.\n\n- (Clarity) The paper is exceptionally written, and the method clear to understand. \n\n- (Significance) As mentioned in the paper, Dynamical Bayesian optimization is an understudied area, however I can think of many possible applications of it. I believe the method to be potentially significant in many applications.", "weaknesses": "- Perhaps the biggest weakness is a lack of comparison, and overall lack of discussion of sparse Gaussian processes and their relationship to online learning. Indeed, methods for selecting inducing points share similarities with the proposed algorithm, e.g. by greedily making the data-set smaller based on information criterions (see Section 9 of [1]) and even with application to online learning and BO [2, 3, 4] (the final reference is very recent, so it was impossible for the authors to include, but it looks relevant). Such methods seem to have the potential of performing well in this setting, and if possible it should be compared against them, but at least they should be discussed.\n\n- The proposed approximation of the 2-Wasserstein distance criterion is potentially loose, however, the authors recognize this and the empirical algorithmic performance is strong.\n\n- A lot of the synthetic testing is carried out in non-dynamical benchmarks. While the importance of the method is that it is able to maintain a small data-set and still carry out BO well, I believe a larger pool of experiments that includes the main application domain is important.\n\n[1] Quinonero-Candela, Joaquin, and Carl Edward Rasmussen. \"A unifying view of sparse approximate Gaussian process regression.\" The Journal of Machine Learning Research 6 (2005): 1939-1959.\n\n[2] Galy-Fajou, Théo, and Manfred Opper. \"Adaptive inducing points selection for gaussian processes.\" arXiv preprint arXiv:2107.10066 (2021).\n\n[3] Moss, Henry B., Sebastian W. Ober, and Victor Picheny. \"Inducing point allocation for sparse Gaussian processes in high-throughput Bayesian optimisation.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2023.\n\n[4] Maus, Natalie, et al. \"Approximation-Aware Bayesian Optimization.\" arXiv preprint arXiv:2406.04308 (2024).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Dear Reviewer NVX3,\n\nThank you for the detailed review. We are discussing below the weaknesses and questions you have raised. Also, please make sure to read the global response as we discuss some of your questions there.\n\n**Lack of discussion of sparse Gaussian processes.**\n\nThank you for pointing out this very interesting connection between the sparse GP literature and our work.\n\nAlthough these works are definitely related to our problem (because they also seek to discard some observations while still preserving the quality of the GP inference), applying them to a dynamic setting would require some non-trivial modifications. More precisely,\n\n(i) In their current form, these works introduce solutions that take place into a spatial domain only. Consequently, they try to approximate an exact GP over a bounded domain, and assume that they can place inducing inputs arbitrarily in this very same domain. Conversely, in the dynamic setting, the domain is constantly changing and the domain where the inducing inputs can be located is restricted to the cartesian product of the spatial domain and the past $[0, t_0]$, with $t_0$ the present time. Moreover, the exact GP should be approximated on a different, complementary, unbounded domain, that is the cartesian product between the spatial domain and the future $[t_0, +\\infty]$.\n\n(ii) Most sparse GP solutions have an hyperparameter that controls the number of inducing inputs to place. Unfortunately, this hyperparameter should be set specifically for each objective function. Conversely, W-DBO can adapt its dataset size depending on the nature of the objective function, without the need to adjust its hyperparameter.\n\nStill, we acknowledge the connections described by the reviewer. Upon acceptance, we will definitely discuss the mentioned papers in Section 2 of the camera-ready version.\n\n**The proposed approximation of the 2-Wasserstein distance criterion is potentially loose.**\n\nWe plan to study more the approximation of the ratio of Wasserstein distances in future works, in order to better understand it theoretically.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
Dear Reviewers,

We thank you all for your detailed reviews. In this global response, we address some weaknesses and questions that were raised in more than one review. We also discuss the additional figures shown in the PDF uploaded with the global rebuttal.

**This is a high-throughput problem, BO is not the ideal candidate for this setting.**

Most evaluations of DBO algorithms are performed in a setting that assumes a constant time step between two consecutive iterations (e.g., see [1, 2]). The underlying assumption is that the objective function is expensive to evaluate, making the GP inference time complexity negligible.

In Section 2, we argue that this is an unreasonable assumption to make for two reasons, regardless of how expensive an evaluation of the objective function is:

(i) Unlike static BO, the GP inference time in DBO limits how much two consecutive observations can be correlated. This in turn impacts the quality of the GP inference.

(ii) Considering optimization tasks with very large (e.g., infinite) time horizons, the GP inference time will eventually become more expensive than the evaluation of the objective function itself, if nothing is done to prevent the dataset size from diverging.

Now, if the objective function is expensive to evaluate, then a sample-efficient strategy should still be used, and BO is indeed the gold standard in this area. In other words, our problem setting still calls for sample-efficient optimization techniques (e.g., BO) because the objective function is expensive to evaluate. Our problem setting is also more general than the one considered in [1, 2], since it relaxes the assumption of a fixed time step between two consecutive iterations.

**The benchmarks are synthetic and not time-dependent.**

Among 12 benchmarks, 10 are synthetic while two are real-world experiments (WLAN and Temperature).

Although the dynamic nature of the two real-world experiments is quite explicit, we agree that we did not stress enough the fact that all synthetic examples are also time-dependent. Time is simply taken as the last dimension (*i.e.*, the $d$-th dimension) of each  $d$-dimensional synthetic function, while the first $(d-1)$ dimensions form the spatial domain. This explains why GP-UCB is not the best-performing strategy on the synthetic benchmarks.

This is stated in Appendix G.1, but not clearly enough in the main text. We will make this point clear in the camera-ready version, by mentioning at the beginning of Section 5 that the temporal domain is always the $d$th dimension of each $d$-dimensional objective function. To further avoid confusion, we will also rewrite Appendix G.2 to make the time variable $t$ explicitly appear in the definition of each synthetic function.

**Additional experiments with an ARD kernel**

In the PDF uploaded with the rebuttal, you can find additional experiments run with the anisotropic SE kernel (see Figure 1). The results are quite similar to the isotropic case, although the average response time and the average regret of each solution are larger. Both of these observations can be explained by the fact that each DBO solution has to infer more hyperparameters (more precisely, three spatial lengthscales instead of one) at each iteration.

We plan to add these results to a dedicated Appendix in the camera-ready version. Thank you for pointing out that it would be beneficial to extend our work to the anisotropic case.

**1D example motivating the normalization by an empty GP.**

A 1D example could indeed help motivating the normalization by an empty GP. You can find it in the PDF uploaded with the rebuttal. We discuss it below.

Figure 2 illustrates that the very same non-normalized Wasserstein distance ($W_2(\mathcal{GP}\_\mathcal{D}, \mathcal{GP}\_{\tilde{\mathcal{D}}}) = 0.46$ in both cases) can lead to different results on the posteriors, depending on their covariance parameters (here, only the lengthscale is considered). The left part of Figure 2 depicts two widely different posteriors, while the right part of Figure 2 depicts two similar posteriors. A good metric of relevancy should capture this difference.

Figure 3 depicts four different scenarios illustrating that the normalized Wasserstein distance is able to capture this difference. In fact, when the metric is low (*i.e.*, close to $0$), it always implies that the two posteriors are similar, regardless of the covariance hyperparameters (*i.e.*, the lengthscale in this example). Similarly, when the metric is large (*i.e.*, close to $1$), the two posteriors are very different, regardless of the covariance hyperparameters.

Upon acceptance, we plan to add this example to Appendix E of the camera-ready version. Thank you for pointing out that a 1-dimensional example would be beneficial.

**References**

[1] Ilija Bogunovic, Jonathan Scarlett, and Volkan Cevher. Time-varying gaussian process bandit optimization. In Artificial Intelligence and Statistics, pages 314–323. PMLR, 2016.

[2] Paul Brunzema, Alexander von Rohr, Friedrich Solowjow, and Sebastian Trimpe. Event-triggered time-varying bayesian optimization. arXiv preprint arXiv:2208.10790, 2022.

Meta-review: All except one of the reviewers are positive about this work. The reviewer with a negative rating acknowledges a number of the work's strengths, such as solving a relevant problem, proposing an intuitive solution, and having good theory, but opted for a very low rating in spite of this due to concerns with benchmarking.

After fairly extensive discussion with other reviewers, no-one else found these concerns sufficiently significant to change their scores, and the paper's other reviewers have all opted in favor of acceptance. Between the two camps, I am inclined to follow the majority's recommendation.
Decision: Accept (poster)