This Too Shall Pass: Removing Stale Observations in Dynamic Bayesian Optimization

Authors: anthony  bardou
patrick  thiran
giovanni  ranieri
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 10:29:28

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.

Main paper: 
 introduction: 
 Many real-world problems require the optimization of a costly-to-evaluate objective function f : S ⊆ R d → R with an unknown closed form (i.e., either the closed form expression of f exists but remains unknown to the user, or it does not exist). Such a setting occurs frequently, and examples can be found in hyperparameters tuning [1] , networking [2, 3] , robotics [4] or computational biology [5] . In such applications, f can be seen as a black-box and cannot be optimized by usual first-order approaches. Bayesian Optimization (BO) is an effective framework for black-box optimization. Its core idea is to leverage a surrogate model, usually a Gaussian Process (GP), to query f at specific inputs. By doing so, a BO algorithm is able to simultaneously discover and optimize the objective function.Since its inception, BO has proven to be very effective at optimizing black-boxes in a variety of contexts, such as high-dimensional input spaces [6, 7, 8] , batch mode [9] or multi-objective optimization [10] . However, few works study BO in dynamic contexts (i.e., with a time-varying objective function), despite its critical importance. Indeed, dynamic black-box optimization problems arise whenever an optimization task is conducted within an environment that comprises exogenous factors that may vary with time and significantly impact the objective function. Dynamic black-boxes are found in network management [11] , unmanned aerial vehicles tasks [12] , hyperparameter tuning in online deep learning [13] , online clustering [14] or crossing waypoints location in air routes [15] .In a dynamic context, f : S × T → R is a time-varying, black-box, costly-to-evaluate objective function with spatial domain S ⊆ R d and temporal domain T ⊆ R. Unlike common black-box objective functions that take as input a point x ∈ S in a spatial domain only, the function f takes as input a point (x, t) ∈ S × T in space and time. Since the problem of optimizing f is 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 an extra dimension to S. It changes the nature of the optimization problem on at least three aspects: (i) at present time t 0 , a DBO algorithm can query any arbitrary point x ∈ S but cannot query past (i.e., t < t 0 ) nor future 1 (i.e., t > t 0 ) points in time, (ii) as time (only) moves forward, a previously collected observation ((x, t), f (x, t)) 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 it 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.Interestingly, (ii) and (iii) imply that each observation eventually becomes stale and, as such, a 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 are twofold. First, we propose a 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 an online fashion 2 .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 a DBO algorithm that exploits the criterion to improve upon state-of-the-art. These contributions can be found in Sections 3, 4 and 5.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in 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. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
background: 
 Let us start by describing the BO framework with a GP prior, as introduced by [16] . Given an objective function f : S ⊆ R d → R, BO assumes that f is a GP (0, k(x, x ′ )). For any x ∈ S, and given a dataset of observations D = {(x i , y i )} i∈ 1,n , where x i ∈ S is a previously queried input with (noisy) function value y i = f (x i ) + ϵ, ϵ ∼ N 0, σ 2 , the posterior distribution of f (x) is N (µ(x), σ 2 (x)), whereEQUATIONEQUATIONwith X = (x 1 , • • • , x n ), y = (y 1 , • • • , y n ), ∆ = k(X, X) + σ 2 I and k(X , Y) = (k (x i , x j )) xi∈X ,xj ∈Y .To find x n+1 , the input to query at the (n + 1)th iteration, a BO algorithm exploits an acquisition function φ : S → R. The acquisition function φ quantifies the benefits of querying the input x 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 [18] or Probability of Improvement [19] . Formally, the BO algorithm determines its next queried input by finding x n+1 = arg max x∈S φ(x).BO extends naturally to dynamic problems, by adapting the covariance function k to properly capture temporal correlations (discussed later in this section). The resulting inference formulas are very similar to (1) and (2) . Given a dataset of observations D = {((x i , t i ) , y i )} i∈ 1,n , where (x i , t i ) ∈ S × T is a previously queried input with the (noisy) function value y i = f (x i , t i ) + ϵ, the posterior distribution of f (x, t) is N µ(x, t), σ 2 (x, t) for any (x, t) ∈ S × T , with µ(x, t) = k ⊤ ((x, t) , X)∆ -1 y,σ 2 (x, t) = k((x, t) , (x, t)) -k ⊤ ((x, t) , X)∆ -1 k((x, t) , X)with X = ((x 1 , t 1 ) , • • • , (x n , t n )), y = (y 1 , • • • , y n ), ∆ = k(X, X) + σ 2 I and k(X , Y) = (k ((x i , t i ) , (x j , t j ))) (xi,ti)∈X ,(xj ,tj )∈Y .Exploiting the usual acquisition functions in a dynamic context is also straightforward. Since a DBO algorithm can only query an input at the current running time t 0 , the next queried input is simply (x n+1 , t 0 ), with x n+1 = arg max x∈S φ(x, t 0 ). Some DBO algorithms (e.g., [20] ) extend the querying horizon to the near future, that is, from t 0 to a time interval [t 0 , t 0 + δ t ]. In that case, the next queried input is (x n+1 , t n+1 ) = arg max (x,t)∈S×[t0,t0+δt] φ(x, t).BO is an active field of research, but relatively few works address DBO, despite the natural extension of BO to dynamic problems described above. We conclude this section by reviewing them. In [21] , the objective function is allowed to evolve in time according to a simple Markov model, controlled by a hyperparameter ϵ ∈ [0, 1]. On the one hand, the authors propose R-GP-UCB, which handles data staleness by resetting the dataset every N (ϵ) iterations. On the other hand, the authors also propose TV-GP-UCB that incorporates data staleness by weighing the covariance of two queries q i = (x i , t i ) and q j = (x j , t j ) by (1 -ϵ) |i-j|/2 . 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 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, it also overlooks a critical effect of the response times of their algorithms. In fact, the response time of a BO algorithm heavily depends on its dataset size n, since BO inference is in O(n 3 ). Although it is reasonable to ignore this for classical BO because the objective function f 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 [21, 23, 24] ) recommend to keep all the collected observations in their datasets, whereas in practice, their response times would asymptotically become prohibitive. Other algorithms (e.g., 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((x, t), (x′ , t ′ )) = k S (||x -x ′ || 2 )k T (|t -t ′ |)to accurately model complex spatio-temporal correlations and samples the objective function only when deemed necessary. Although this reduces the size of ABO'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 an observation can be found in the sparse GPs literature (e.g., see [25, 26] ). However, they require non-trivial adjustments to account for the particular nature of the time dimension. As far as we know, there is no 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 to 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.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: The limitations are explicitly discussed in Appendix I and the approximation error of our proposed upper bounds is discussed at length in Appendix C.Guidelines:• 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 a 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 a specific model and empirical evaluation, it may be necessary to either make it 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 submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
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 a surrogate model for f (see [16] ).Assumption 3.1. f is a GP (0, k((x, t) , (x ′ , t ′ )), whose mean is 0 (without loss of generality) and whose covariance function is denoted by k : S × T × S × T → R + .In order to accurately model complex spatio-temporal dynamics, we make the same assumption on the decomposition and isotropy in time and space of the covariance function k as in [20] . Assumption 3.2.EQUATIONwith λ > 0, k S : R + → [0, 1] and k T : R + → [0, 1] two positive correlation functions, parameterized by lengthscales l S > 0 and l T > 0, respectively. The factor λ > 0 scales the product of the two correlation functions and hence, controls the magnitude of the covariance function k. The lengthscales l S and l T 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 k T , it does not accurately measure the relevancy of an observation. The next section addresses this question.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 i = ((x i , t i ), y i ) by measuring the impact that the removal of o i has on the GP posterior.Let GP D be the GP conditioned on the datasetD = {((x i , t i ), y i )} i∈ 1,n , with (x i , t i ) ∈ S × T and y i = f (x i , t i ) + ϵ, ϵ ∼ N 0, σ 2 .Without loss of generality, let us measure the impact of removing ((x 1 , t 1 ), y 1 ) from the dataset on the GP posterior. Clearly, the measure should be defined on the domain of future predictions at time t 0 , denoted by F t0 , which must include the whole space S and only the future time interval [t 0 , +∞):EQUATIONWe compare a GP conditioned on the whole dataset, denoted by GP D , with a GP conditioned on D, the dataset without (x 1 , t 1 , y 1 ), denoted by GP D . For an arbitrary point (x, t) ∈ F t0 , GP D provides a posterior distribution N D (x, t) = N µ D (x, t), σ 2 D (x, t) , and so does GP D with N D (x, t) = N µ D (x, t), σ 2 D (x, t) . We compare these two distributions by using the 2-Wasserstein distance [28] , given byEQUATIONA natural extension of the 2-Wasserstein distance from a point (x, t) ∈ F t0 to the domainEQUATION)Observe that ( 8) is a criterion that effectively captures the impact of removing the observation o 1 = ((x 1 , t 1 ), y 1 ) from the dataset on the GP posterior. However, as discussed in Appendix F, the covariance function hyperparameters θ = (λ, l S , l T ) control the magnitude of ( 8). This is illustrated by Figure 1 , which depicts two couples of GP posteriors that achieve the same value (8) . Depending on the lengthscale magnitude, the posteriors may be quite different or, conversely, very similar. As a result, (8) cannot be directly used as a gauge of observation relevancy.To remove this ambiguity, we normalize (8) by W 2 (GP D , GP ∅ ) (i.e., the 2-Wasserstein distance between the GP conditioned on D and the prior GP), and we obtain the ratioEQUATIONIntuitively, W 2 (GP D , GP ∅ ) 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 a direct consequence, , bottom) row depicts couples of posteriors that yield a small (resp., large) ratio (9) . The left (resp., right) column depicts couples of posteriors controlled by a small (resp., large) lengthscale. The prior GP mean µ ∅ = 0 is shown as a black dashed line, and the Wasserstein distance between the posterior and the prior as a gray shaded area. ( 9) 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 ( 9) is small, one can infer that the observation o 1 can be safely removed from the dataset since it will have virtually no impact on the posterior. Conversely, when ( 9) is large, one can infer that removing o 1 would alter the posterior too much, and conclude that it is a relevant observation that must remain in the dataset. The exploitation of ( 9) is discussed in more details in Section 4.2.The criterion ( 9) is useful for a 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.4 Using the Criterion in Practice 
computational tractability: 
 In [29] , the authors provide an algorithm to approximate the 2-Wasserstein distance between two GPs up to an 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-iteration cost as small as possible to ensure a high sampling frequency. In this section, we put this issue to rest by deriving an explicit approximation of (9) . 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 is to 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 W 2 (GP D , GP ∅ ). Theorem 4.1. Let t 0 be the present time and D = {((x i , t i ), y i )} i∈ 1,n be a dataset of observations made before t 0 . Let D = {((x i , t i ), y i )} i∈ 2,n be the dataset without the first observation and F t0 = S × [t 0 , +∞) be the domain of future predictions. Then, an upper bound forEQUATIONwhereC(X , Y) = (k S * k S )(x j -x i ) • (k T * k T ) +∞ t0-ti (t j -t i ) (xi,ti)∈X (xj ,tj )∈Y, where (f * g) denotes the convolution between f and g, and (f * g) b a denotes the convolution between f and g restricted to the interval [a, b]. The terms a, b, c, E and M are explicited in Appendices A and B.Moreover, an upper bound forEQUATION)This theorem provides the analytic form of an upper bound for the Wasserstein distance between GP D and GP D and the Wasserstein distance between GP D and GP ∅ on the domain of future predictions F t0 . Using it, we can compute an approximation R of the relative criterion (9) , that isEQUATION)In Appendix C, we study the error between the criterion (9) and the approximation (12) . In essence, we bound the approximation error caused by estimating the integrals in (8) by a self-convolution of the covariance functions k S and k T (i.e., k S * k S and k T * k T ). Furthermore, we provide numerical evidence that the approximation errors in the numerator and the denominator of (12) compensate each other at least in part, making (12) a decent approximation for (9) .In practice, the upper bounds (10) and (11) can only be computed efficiently if the convolutions of the covariance functions can themselves be computed efficiently. The analytic forms for the convolution Matérn (ν, l)2 1-ν Γ(ν) √ 2ν||x||2 l S ν K ν √ 2ν||x||2 l Sof two usual covariance functions listed in Table 1 , namely Squared-Exponential (SE) and Matérn [30] , are provided in Appendix D 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 a 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, 4 in Appendix D and Theorem 4.1 show that the approximation of (9) given by ( 12) can be computed efficiently in an online setting. In an effort to generalize our results to a class of covariance functions that extends beyond Assumption 3.2, we also discuss how to compute the selfconvolution of an anisotropic spatial SE covariance function in Appendix E. We now leverage (12) to propose a DBO algorithm able to pinpoint and remove irrelevant observations in its dataset.The metric (9) and its approximation (12) can be seen as a relative error (or drift), expressed as a percentage, that separates GP D and GP D . Indeed, the Wasserstein distance W (GP D , GP D ) is scaled by the Wasserstein distance W (GP D , GP ∅ ), that is, the distance between GP D and the prior.In other words, (9) and its approximation (12) measure the relative drift from GP D to GP D 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 a consequence, removing multiple observations could, in the worst case, make GP D drift exponentially fast from GP D . To keep this exponential drift under control, one can use a removal budget b(t) = (1 + α) t that allows a maximal relative drift from GP D of α per unit of time (e.g., if α = 0.1, the allowed maximal drift is 10 % per unit of time). The cost of removing an observation is given by (12) .Algorithm 1 describes W-DBO, a DBO algorithm exploiting (12) to remove irrelevant observations on the fly. As described above, the removal budget is controlled by a single hyperparameter α and grows exponentially as time goes by (see line 24) . At each iteration, ( 12) is used to compute the relevancy of each observation in the dataset (see lines [10] [11] [12] [13] . The relevancy of the least relevant observation is then compared to the removal budget, and the observation is removed if the budget allows it (see lines [14] [15] [16] [17] . This process is repeated until all the budget is consumed. Such a greedy observation removal policy causes W-DBO to overestimate the impact of removing multiple observations 3 . We discuss and motivate the expression of the removal budget in Appendix G. The sensitivity analysis conducted in Section 5.1 supports this removal budget, by showing that the same value of the hyperparameter α is valid for a 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 1 to take as inputs any GP model GP and any acquisition function φ. Similarly, observe that lines 5-8 in Algorithm 1 describe the usual BO optimization loop, without any modification. 
numerical results: 
 In this section, we 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 a 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 Observe y t = f (x t , t) + ϵ, ϵ ∼ N (0, σ 2 )7:D = D ∪ {(x t , t, y t )} 8:Condition GP on D and get the MLE parameters (λ, l S , l T , σ2 ) for all (x i , t i , y i ) ∈ D [in parallel] do 11:D = D \ {(x i , t i , y i )} 12:Compute Ri using (12) 13:end for b = b(1 + α) (t-t ′ )/l T 25: end while parameters and the noise level are estimated on the fly. Please refer to Appendix H.1 for further experimental details, and to Appendix H.2 for a detailed description of the dynamic benchmarks.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 the algorithm. This evaluation protocol reveals that a trade-off must be achieved between having an accurate model of the objective function (which requires a large dataset) and being able to track the optimal argument of the function as it evolves (which requires a high sampling frequency, thus a moderately-sized 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 in bold text. The performance of algorithms whose confidence intervals overlap the best performing algorithm's confidence interval is underlined. To study this trade-off, we compare the performance of W-DBO with several values of its hyperparameter α, as 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 a trade-off between the size of the dataset and the sampling frequency can significantly improve the performance of W-DBO. Clearly, the sweet spot is reached for α = 1 4 . This hyperparameter value is used to evaluate W-DBO in the next section.The competing baselines are the relevant algorithms reported in Section 2, namely R-GP-UCB and TV-GP-UCB [21] , 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 H.2, and the performance of each DBO solution is discussed at length in Appendix H.3. Furthermore, the provided supplementary animated visualizations are discussed in Appendix H.4.In this section. we only depict the performance of the DBO solutions on the Ackley synthetic function 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 its deep hole at the center of its domain. Observe that most DBO solutions miss that hole, as their average regrets skyrocket between 200 and 2 . 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 the DBO solutions is shown in black. 400 seconds. In contrast, W-DBO manages to rapidly exploit the hole at the center of the function domain, thereby maintaining a low average regret.This performance gap can be explained by studying the dataset size of W-DBO (see the right side of Figure 4 ). 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 a summary of the performance of the DBO solutions across all our benchmarks, please refer to the last row of Table 2 , and to the visual summary in Figure 5 . 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-UCB appears significantly better than the aforementioned algorithms. Remarkably, W-DBO shows significantly better performance than TV-GP-UCB and outperforms the other DBO solutions by a comfortable margin. In fact, it obtains the lowest average regret on almost every benchmark.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes]Justification: A GitHub repository with a minimal working example and a detailed documentation is provided along with the camera-ready version.f Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) 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 (e. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
conclusion: 
 The ability to remove irrelevant observations from the dataset of a DBO algorithm is 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 in an online setting and (iii) a DBO algorithm, W-DBO, that exploits this approximation. We have evaluated W-DBO against the state-of-the-art of DBO on a variety of benchmarks comprising synthetic functions and real-world experiments. The evaluation was conducted in 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-art of DBO by a comfortable 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 erasure of the whole dataset. By doing so, W-DBO simultaneously ensures a high sampling frequency and a very good predictive performance.In addition to its impact on DBO itself, we believe that W-DBO can have a significant impact on the fields that make heavy use of DBO (e.g., computer networks, robotics). As a future work, we plan on exploring these applications of W-DBO. Furthermore, we plan on better understanding the excellent performance of W-DBO by addressing the difficult problem of deriving a regret bound that holds in a continuous time setting and incorporates the effect of the sampling frequency of the DBO algorithm as well as the deletion of irrelevant observations.A Wasserstein Distance at a Point in F t 0 At the core of ( 8) lies (7) . In this appendix, we provide explicit expressions for (7) . Let us start by proving the following lemma. Lemma A.1.EQUATIONwith ∆ D = k(D, D) + σ 2 I and E, F , G, H defined asEQUATIONEQUATIONEQUATIONEQUATIONwherek (X , Y) = (k ((x i , t i ) , (x j , t j ))) (xi,ti)∈X (xj ,tj )∈Y and ∆ D = k( D, D) + σ 2 I.Proof. The proof is trivial using the inverse of a block matrix:EQUATIONNote that A and D must be invertible.We can use (18) to write ∆ -1 D as a function of ∆ -1 D , sinceEQUATIONNote that, in (19) , λ + σ 2 and ∆ D are invertible. Therefore, (18) can be applied to (19) , and this yields the desired result.To compute (7) , we provide the following results. Proposition A.2.EQUATIONEQUATIONwith c = G + H ⊤ and M = F -∆ -1 D .We now prove Proposition A.2.Proof. According to (3)EQUATIONApplying the same definition for µ D (x, t) we haveEQUATIONwhere (23) follows from Lemma A.1.Finally, we haveEQUATIONwitha = Ey 1 + G ỹ and b ⊤ = H ⊤ y 1 + ỹ⊤ F -∆ -1 D .This concludes the proof.Finally, we prove Proposition A.3.Proof. As (σ D (x, t) -σ D (x, t)) 2 is hard to integrate, to get (8) , we upper bound it byEQUATIONEQUATIONwhere (26) follows from σ D (x, t) ≤ σ D (x, t). Now, (4) yields that for X = D,EQUATIONand forEQUATIONEQUATIONwhere (29) follows from Lemma A.1. Developing the dot product in (30), we haveEQUATIONCombining ( 28) and (31), we getEQUATIONwhere c = G + H ⊤ and M = F -∆ -1 D . Combining ( 27) and ( 32) concludes the proof.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: As previously mentioned, all experimental details are discussed in Appendix H. A sensitivity study and a discussion about the hyperparameter of W-DBO are also provided in Section 5.1 and Appendix G, respectively.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: Standard errors are depicted graphically in Section 5 and Appendix H.2.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at 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 or 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 a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: The computing resources are described in Appendix H.1.Guidelines:• The answer NA means that the paper does not include experiments.• 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 of compute required for each of the individual experimental runs as well as 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). Justification: This paper proposes a 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. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (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 as 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 to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? 
 
 Appendix: In this appendix, we extend the computation of the Wasserstein distance between two posterior Gaussian distributions at an arbitrary point in F t0 (see Propositions A.2 and A.3) to the Wasserstein distance between GP D and GP D on F t0 . We also provide an upper bound for the Wasserstein distance between a posterior GP conditioned on D (i.e., GP D ) and the prior GP (i.e., GP ∅ ).Let us start by proving the following lemma. Lemma B.1. Let t 0 be the present time and D = {((x i , t i ), y i )} i∈ 1,n be a dataset of observations before t 0 . Let D = {((x i , t i ), y i )} i∈ 2,n and F t0 = S×[t 0 , +∞) be the domain of future predictions. Then, for any pair of observations (x i , t i , •), (x j , t j , •) from D we have Proof. For the sake of brevity, let us denote by I the LHS of (33) . According to Assumption 3.2, Fubini's theorem and since F t0 = S × [t 0 , +∞), we haveEQUATIONEQUATION)Let us focus on 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 boundEQUATIONwhere (35) holds since k S is positive and (36) comes from the change of variable u = xx i .Obviously, this upper bound remains interesting because k S is usually an exponentially decreasing function (see Table 1 ). We discuss this point in more details in Appendix C.Observe that (36) is actually the convolution (k S * k S )(x j -x i ). A similar change of variable can be made regarding the time integral, with v = t -t i :EQUATIONNote that (37) is also a convolution, but restricted to the interval [t 0 -t i , +∞). We denote this convolution on a restricted interval 34), ( 36) and (37) yields Lemma B.1.(k T * k T ) +∞ t0-ti (t j -t i ). Combining (We now prove the following lemma. Lemma B.2. Let t 0 be the present time and D = {((x i , t i ), y i )} i∈ 1,n be a dataset of observations before t 0 . Let D = {((x i , t i ), y i )} i∈ 2,n and F t0 = S × [t 0 , +∞) the domain of future predictions. Then, we haveEQUATION)with a = Ey 1 + G ỹ and b ⊤ = H ⊤ y 1 + ỹ⊤ F -∆ -1 D , C(X , Y) = (k S * k S )(x j -x i )(k T * k T ) +∞ t0-ti (t j -t i ) (xi,ti)∈X (xj ,tj )∈Y, 1 the conformable vector of ones, (f * g) the convolution between f and g and (f * g) b a the convolution between f and g restricted to the interval [a, b].A.2, S ∞ t0 (µ D (x, t) -µ D (x, t)) 2 dxdt = S ∞ t0 ak((x, t), (x 1 , t 1 )) + bk (x, t), D 2 dxdt (39) with a = Ey 1 + G ỹ and b ⊤ = H ⊤ y 1 + ỹ⊤ F -∆ -1 D .Expanding the square, we get three different integrals, denoted A, B and C withA = S ∞ t0 a 2 k 2 ((x, t), (x 1 , t 1 ))dxdt, B = S ∞ t0 2abk((x, t), (x 1 , t 1 ))k((x, t), D)dxdt, C = S ∞ t0 bk (x, t), D k ⊤ (x, t), D b ⊤ dxdt.Using Lemma B.1, computing an upper bound for A, B and C is immediate. In fact,EQUATIONEQUATION)where ⊙ is the Hadamard product and 1 the conformable vector of ones.Adding (40), ( 41) and (42) together concludes our proof.To get the first part of Theorem 4.1, we prove the following lemma.Lemma B.3. Let t 0 be the present time and D = {((x i , t i ), y i )} i∈ 1,n be a dataset of observations before t 0 . Let D = {((x i , t i ), y i )} i∈ 2,n and F t0 = S × [t 0 , +∞) the domain of future predictions. Then, we have Proof. This proof is conceptually identical to the previous one. By Proposition A.3,EQUATIONwith c = G + H ⊤ , M = F -∆ -1 D , C(X , Y) = (k S * k S )(x j -x i )(k T * k T ) +∞ t0-ti (t j -t i ) (xi,EQUATIONwherec = G + H ⊤ and M = F -∆ -1 D .By linearity of the integral, the RHS of (44) can be split into three different integrals A, B and C whereA = S ∞ t0 Ek((x, t), (x 1 , t 1 )) 2 dxdt B = S ∞ t0 k((x, t), (x 1 , t 1 ))ck((x, t), D)dxdt C = S ∞ t0 k ⊤ ((x, t), D)M k((x, t), D)dxdtOnce again, Lemma B.1 can be used to compute an upper bound for A, B and C. In fact,EQUATIONEQUATIONEQUATION)where ⊙ is the Hadamard product and 1 is the conformable vector of ones.Adding (45), ( 46) and (47) together concludes the proof. Proof. Recall that, according to GP ∅ , f (x, t) ∼ N (0, λ) for any point (x, t) ∈ F t0 . Consequently,EQUATIONwith a = ∆ -1 y and C(X , Y) = (k S * k S )(x j -x i )(k T * k T ) +∞ t0-ti (t j -t i ) (xi,EQUATION)These two integrals in (49) can be computed with the same techniques as above. For the mean integral, we haveEQUATION)Regarding the variance integral in (49), we haveEQUATIONwhere (51) holds because √ λσ D (x, t) ≥ σ 2 D (x,In Appendix B, we provided a computationally tractable upper bound of W 2 (GP D , GP D ) and W 2 (GP D , GP ∅ ). In this appendix, we provide the expression of the corresponding approximation error and we study its magnitude with the Squared-Exponential (SE) covariance function.First, and without loss of generality, assume S = [0, 1] d . In Appendix B, we have to approximate the Wasserstein distance because the integration over S in (34) is difficult to compute. The upper bound proposed in ( 35) isEQUATION)Recall that the upper bounded quantity is a product of functions which are decreasing exponentially (see Table 1 ). As a consequence, their product decreases exponentially as well, so that extending the integration from S = [0, 1] d to R d has a bounded impact on the result. Clearly, a first absolute approximation error for ( 54) is(R\[0,1]) d k S (||x -x i || 2 )k S (||x -x j || 2 )dx.Because the upper bounded quantity is a product of two correlation functions k S 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:EQUATION)Obtaining a 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 to l S . Conversely, the right term in (55) is a decreasing function with respect to l S . This observation allows us to derive the spatial lengthscale for which (55) is maximal.Proposition C.1. Let (x i , x j ) ∈ S 2 , with S = [0, 1] d . Let k S be a SE kernel with lengthscale l S .Then,EQUATION)with W 0 the principal branch of the Lambert function.Proof. Because the two terms in (55) are respectively increasing and decreasing with respect to the spatial lengthscale l S , (55) is maximal when both terms are equal. Therefore, from (55), we have the following relation where (57) uses a result derived in Appendix D and reported in Table 3 . Solving for l S concludes the proof.EQUATIONLet us illustrate (55) and Proposition C.1 by plotting the relative approximation error (55) with respect to the spatial lengthscale l S . The integral over S is computed with a 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 a spatial lengthscale l * S given by (56), as anticipated above. For l S < l * S , 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 l S > l * S . The bottom row of Figure 6 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 l S 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 Table 3 : Analytic forms for the convolution of usual spatial covariance functions. Γ is the Gamma function, J α is a Bessel function of the first kind of order α, K α is a modified Bessel function of the second kind of order α.Covariance Function k S (k S * k S )(x) Squared-Exponential (l S ) π d 2 l d S e -||x|| 2 2 4l 2 S Matérn (ν, l S ) 2 d 2 -2ν+1 π d 2 Γ(ν+ d 2 ) 2 Γ(ν) 2 Γ(2ν+d) √ 2ν l S 2ν-d 2 ||x|| 2ν+ d 2 2 K 2ν+ d 2 ||x||2√ 2ν l S Table 4 : Analytic forms for the convolution of usual temporal covariance functions on the interval [t 0 -t i , +∞). Note that erf is the error function. Also, for the sake of brevity, the terms C k1k2 , P k1k2 and Q k1k2 are defined in Appendix D.2.2.Covariance Function k T (k T * k T ) +∞ t0-ti (t j -t i ) Squared-Exponential (l T ) √ πl T 2 e -(t i -t j ) 2 2l 2 T 1 -erf 2t0-ti-tj 2l T Matérn (ν = p + 1 2 , l T ) p k1=0 p k2=0 C k1k2 e - √ 2p+1(2t 0 -t i -t j ) l T P k1k2 (t 0 , t i , t j )bounds). As a consequence, the approximation errors compensate each other (at least in part) when computing (12) . To verify this observation numerically, we compute the relative approximation error between the criterion (9) and its approximation (12) . The results are 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 d or l S . In the main paper, Section 5 corroborates this observation by demonstrating the usefulness of the approximation (12) in practice.In this appendix, we derive the analytic forms of the convolution of usual covariance functions listed in Tables 3 and 4 , which are used to compute the criterion (12) .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 is(k * k)(x) = F -1 F 2 (k) (x), with F(f ) denoting the Fourier transform of f and F -1 (f ) the inverse Fourier transform of f .The Fourier transform F(k) of a stationary covariance function k is called the spectral density of k, and is usually denoted S. Both functions are Fourier duals of each other (see [27] for more details). Furthermore, it is known that if k is isotropic (i.e. it can be written as a function of r = ||x|| 2 ), then its spectral density S(s) can be written as a function of s = ||s|| 2 . In that case, the two functions are linked by the pair of transforms (see [27] )EQUATION)where J d 2 -1 is a 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.Corollary D.1. Let k be a stationary, isotropic covariance function with spectral density S. Let r = ||x|| 2 and s = ||s|| 2 . Then,(k * k)(r) = 2π r d 2 -1 ∞ 0 S 2 (s)J d 2 -1 (2πrs)s d 2 ds (60)where J d 2 -1 is a Bessel function of the first kind and of order d/2 -1. We now derive the analytic forms for the convolutions of usual spatial covariance functions k S .Lemma D.2. Let k S be a Squared-Exponential covariance function (see Table 1 ), with lengthscale l S > 0. Then,EQUATIONProof. The spectral density of a Squared-Exponential covariance function k S is (see [27] )EQUATION)According to Corollary D.1, we haveEQUATIONwhere r = ||x|| 2 and J d 2 -1 is a Bessel function of the first kind of order d 2 -1. It is known (see [31] ) that∞ 0 e -αx 2 x ν+1 J ν (βx)dx = β ν (2α) ν+1 e -β 2 4α .Therefore,EQUATIONReplacing r by ||x|| 2 in (64) concludes the proof.Lemma D.3. Let k S be a Matérn covariance function (see Table 1 ), with smoothness parameter ν > 0 and lengthscale l S > 0. Then,EQUATIONwhere Γ is the Gamma function and K α is a modified Bessel function of the second kind of order α.Proof. The spectral density of a Matérn covariance function k S is (see [27] )S(s) = 2 d π d 2 Γ(ν + d 2 )(2ν) ν Γ(ν)l 2ν S 2ν l 2 S + 4π 2 s 2 -ν-d 2 (66)where Γ is the Gamma function.According to Corollary D.1, we haveEQUATION)where J d 2 -1 is a Bessel function of the first kind of order d 2 -1, and (67) comes from the of variable u = 2πs.It is known (see [31] ) that∞ 0 (a 2 + x 2 ) -(µ+1) J α (bx)x α+1 dx = a α-µ b µ 2 µ Γ(µ + 1) K µ-α (ab)where K µ-α is a modified Bessel function of the second kind of order µ -α.Therefore,EQUATIONReplacing r by ||x|| 2 in (68) concludes the proof.In this section, we derive the analytic expression of the convolutions of the most popular temporal covariance functions k T restricted to the interval [t 0 -t i , +∞). Therefore, we compute many integrals of the form+∞ t0-ti k T (|t|)k T (|t j -t i -t|)dt,which can be rewritten asEQUATION)since t ≥ 0 and t j -t i -t ≤ 0 for all t ∈ [t 0 -t i , +∞). The form (69) will be used in every proof of this section.Lemma D.4. Let k T be a Squared-Exponential covariance function (see Table 1 ), with lengthscale l S > 0. Then,EQUATION)where erf is the error function.Proof. Since k T is a Squared-Exponential function, (69) becomesEQUATIONIt is known (see [31] ) thate -(ax 2 +2bx+c) dx = 1 2 π a e b 2 -ac a erf √ ax + b √ awhere erf is the error function.Therefore, (71) becomes√ πl T 2 e -(t i -t j ) 2 2l 2 T 1 -erf 2t 0 -t i -t j 2l T .Lemma D.5. Let k T be a Matérn covariance function (see Table 1 ), with smoothness parameter ν = p + 1 2 , p ∈ N and lengthscale l T > 0. Then,EQUATIONwhereEQUATIONEQUATIONEQUATIONand P (k) the kth derivative of P (t) with respect to t.Proof. The Matérn covariance function has a simpler form when its smoothness parameter ν is a half-integer, that is ν = p + 1 2 , p ∈ N (see [30] ). In that case,k T (t) = e - √ 2p+1t l T p! (2p)! p k1=0 (p + k 1 )! k 1 !(p -k 1 )! 2 √ 2p + 1t l T p-k1.Therefore,EQUATIONwith C k1k2 defined in (73) and P (t) defined in (75).Integrating (76), we get2 √ 2p + 1 l T p k1=0 p k2=0 C k1k2 e - √ 2p+1(t i -t j ) l T +∞ t0-ti e -2 √ 2p+1t l T P (t)dt (77)thanks to the linearity of the integral.It is known (see [31] ) thatP (x)e ax dx = e ax a m k=0 (-1) k P (k) m (x) a kwhere P m is a polynomial of degree m and P (k) m is the kth derivative of P m . Therefore,EQUATIONCombining ( 77) and (78) we get(k T * k T ) +∞ t0-ti (t i -t j ) = p k1=0 p k2=0 C k1k2 e - √ 2p+1(2t 0 -t i -t j ) l T P k1k2 (t 0 , t i , t j )with C k1k2 defined in (73) and P k1k2 defined in (74).This concludes the proof.In this appendix, we illustrate how Theorem 4.1 could be extended to anisotropic spatial kernels by considering an Automatic Relevance Detection (ARD) Squared-Exponential (SE). It has the following form:EQUATIONwhereM = diag (l 1 , • • • , l d) is a diagonal matrix that gathers a different lengthscale for each dimension. Observe that the isotropic SE with lengthscale l S is retrieved by setting M = l S I.Because the ARD SE kernel (79) is anisotropic, the convolution with itselfEQUATIONcannot be simplified to a one-dimensional integral through a change to polar coordinates, as done in Corollary D.1. The integral becomes more complex, but can still be computed exactly for some kernel such as the ARD SE. Lemma E.1. Let k S be an ARD SE covariance function with parameter M . Then,EQUATIONProof. For the ARD SE kernel with parameterEQUATION)Integrating the k-th term in (82), we getEQUATION)Injecting ( 83) into (82), we haveEQUATIONwhere (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 a special case of Lemma E.1 where M = l S I, that is, when k S is an isotropic SE kernel.In this appendix, we discuss how ( 9) and its approximation ( 12) address the dependency on the covariance function hyperparameters introduced by ( 8). For the sake of this discussion, we take k S and k T as two Squared-Exponential (SE) covariance functions (see Table 1 ). A similar reasoning can be conducted with Matérn covariance functions.Let us start by rewriting the product of spatial and temporal convolutions C((x, t), (x ′ , t ′ )) with the formulas provided in Tables 3 and 4 for the SE covariance functions. We getEQUATION)The dependency on the covariance function hyperparameters λ, l S , l T appears clearly in (86). Both l S and l T are used, not only to scale the spatial distance ||x -x ′ || 2 and the temporal distance |t -t ′ | in C * , but also as a scaling constant of the magnitude of the output of the product of convolutions itself. Because C((x, t), (x ′ , t ′ )) is involved in every term of ( 10) and ( 11) in Theorem 4.1, 1 2 π d+1 2 l d S l T can be factored out of (10) and (11) . Overall, both equations have in common the factor 1 2 π d+1 2 λ 2 l d S l T . Clearly, this shows how the covariance function hyperparameters θ = (λ, l S , l T ) 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 ( 10) and (11) . Doing so, the factors which are common to the two equations cancel out. Considering the application of Theorem 4.1 with k S and k T being SE covariance functions, the undesirable factor 1 2 π d+1 2 λ 2 l d S l T is removed. Clearly, ( 12) remains a function of l S and l T , but the hyperparameters are only used to scale the spatial and temporal distances, that is to control correlation lengths. However, the undesirable scaling exposed in (86) no longer exists.In this appendix, we discuss why the removal budget of W-DBO (see Algorithm 1), denoted b t at a given time t, has the formb 0 = 1, b t+∆t = b t (1 + α) ∆t/l T (87)with l T the temporal lengthscale (see Assumption 3.2) and α the hyperparameter of W-DBO. Crucially, l T and t must be expressed in the same unit of time. First, note that the expression of the budget is intuitive because (9) measures a ratio, expressed as a percentage. Therefore, the budget must accumulate in a multiplicative way, leading to the exponential form (87).More interestingly, let us discuss the exponent ∆t/l T . Arguably, an alternative, more intuitive form of the removal budget would beEQUATION)Although easier to understand, the budget (88) presents a major problem since it depends on arbitrary choices made by the user. This is illustrated by Table 5 , 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 a recommendation for the hyperparameter α. If the analysis in Section 5.1 had used the budget (88), its recommendation α * 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 5.1, using the budget (87), is a single number that is valid regardless of experiment duration. This is a much more general insight.In each experiment, the d-dimensional spatial domain is scaled in S ′ = [0, 1] d and the temporal domain (viewed as the (d + 1)th dimension) is normalized in [0, 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 its own temporal covariance function. Eventually, ABO and W-DBO exploits a Matern-3/2 kernel as their temporal covariance function.Each DBO algorithm begins its optimization task with 15 initial observations, uniformly sampled in S ′ × 0, 1 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 ∆t. Clearly, ∆t is a function of the dataset size of the DBO algorithm. Consequently, it varies throughout the optimization, getting larger when the DBO algorithm adds a 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 f at a specific time in the future) and a Gaussian noise with variance equal to 5 % of the signal variance is added. Then, the next iteration begins at time t + ∆t (except for ABO if it decides to sample f later). 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 computationally-heavy 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 on a laptop equipped with an Intel Core i9-9980HK @ 2.40 GHz with 8 cores (16 threads).We provide here a 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 d + 1dimensional spatio-temporal domain S × T . More precisely, the point z is explicitly given by z = (x 1 , • • • , x d , t). Also, we will write d ′ = d + 1 for the sake of brevity.Rastrigin. The Rastrigin function is d ′ -dimensional, and has the formf (z) = ad ′ + d ′ i=1 z 2 i -a cos (2πz i ) .For the numerical evaluation, we set a = 10, d ′ = 5 and we optimized the function on the domain [-4, 4] d ′ . The results are provided in Figure 8 .Schwefel. The Schwefel function is d ′ -dimensional, and has the formf (z) = 418.9829d ′ - d ′ i=1 z i sin |z i | .For the numerical evaluation, we set d ′ = 4 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 our results using the ARD covariance function studied in Appendix E. The results are provided in Figure 10 .Styblinski-Tang. The Syblinski-Tang function is d ′ -dimensional, and has the form Eggholder. The Eggholder function is 2-dimensional, and has the formf (z) = 1 2 d ′ i=1 z 4 i -16z 2 i + 5z i .f (z) = -(z 2 + 47) sin |z 2 + z 1 2 + 47| -z 1 sin |z 1 -z 2 -47| .For the numerical evaluation, we optimized the function on the domain [-512, 512] 2 . The results are provided in Figure 12 .Ackley. The Ackley function is d ′ -dimensional, and has the formf (z) = -a exp   -b 1 d ′ d ′ i=1 z 2 i   -exp   1 d ′ d ′ i=1 cos(cz i )   + a + exp(1).For the numerical evaluation, we set a = 20, b = 0.2, c = 2π, 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 in Appendix E. The results are provided in Figure 14 . Rosenbrock. The Rosenbrock function is d ′ -dimensional, and has the formf (z) = d ′ -1 i=1 100(z i+1 -z 2 i ) 2 + (z i -1) 2 .For the numerical evaluation, we set d ′ = 3 and we optimized the function on the domain [-1, 1.5] d ′ . The results are provided in Figure 15 .Shekel. The Shekel function is 4-dimensional, and has the formf (z) = - m i=1   4 j=1 (z j -C ji ) 2 + β i   -1 .For the numerical evaluation, we set m = 10, β = Hartmann-3. The Hartmann-3 function is 3-dimensional, and has the formf (z) = - 4 i=1 α i exp   - 3 j=1 A ij (z j -P ij ) 2   .For the numerical evaluation, we set α = (1.0, 1.2, 3.0, 3. and we optimized the function on the domain [0, 1] 3 . The results are provided in Figure 17 .Hartmann-6. The Hartmann-6 function is 6-dimensional, and has the form and we optimized the function on the domain [0, 1] 6 . The results are provided in Figure 18 .f (z) = - 4 i=1 α i exp   - 6 j=1 A ij (z j -P ij ) 2   .Powell. The Powell function is d ′ -dimensional, and has the formf (z) = d ′ /4 i=1 (z 4i-3 + 10z 4i-2 ) 2 + 5(z 4i-1 -z 4i ) 2 + (z 4i-2 -2z 4i-1 ) 4 + 10(z 4i-3 -z 4i ) 4 .For the numerical evaluation, we set d ′ = 4 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 a famous benchmark, used in other works such as [21, 22] . The goal of the DBO task is to activate the sensor with the highest temperature, which will vary with time. To make the benchmark more interesting, we interpolate the data in space-time. With this interpolation, the algorithms can activate any point in space-time, making it a 3-dimensional benchmark (2 spatial dimensions for a location in Intel Research Berkeley, 1 temporal dimension). For the numerical evaluation, we used the first day of data. The results are provided in Figure 20 .WLAN. This benchmark aims at maximizing the throughput of a Wireless Local Area Network (WLAN). 18 moving end-users are associated with one of 4 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 to improve its performance. To do so, each node has a power level 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.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 a theoretical upper bound on the throughput of a wireless communication. We denote it C(i, j), we express it in bits per second (bps). It depends on S ij the Signal-to-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):C ij (x, t) = W log 2 (1 + S ij (x, t)).Then, the objective function isf (x, t) = 4 i=1 j∈Ni C ij (x, t),with N i the end-users associated with node i. For the numerical evaluation, we optimized the power levels x in the domain [10 0.1 , 10 2.5 ] 4 . For this experiment, the DBO solutions were evaluated with a Matérn-5/2 for the spatial covariance function and a Matérn-1/2 for the temporal covariance function. The results are provided in Figure 21 .In this section, we discuss the performance achieved by all the DBO solutions on the benchmarks introduced in the previous section. GP-UCB. This baseline, which does not take into account temporal correlations, obtains surprisingly 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 to estimate with MLE (i.e., λ, l S , σ 2 ). Overall, it is dominated by R-GP-UCB, TV-GP-UCB and W-DBO, but behaves surprisingly well against ABO and ET-GP-UCB (see Figure 5 ).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 l T away), overestimating l T 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 l T , e.g., Rosenbrock (see Figure 15 ) and Powell (see Figure 19 ). Conversely, ABO obtains competitive performance on functions with smaller temporal lengthscales, e.g. Eggholder (see Figure 12 ) or Shekel (see Figure 16 ). These results highlight the lack of robustness of ABO.ET-GP-UCB. Like GP-UCB, ET-GP-UCB does not take into account temporal correlations, and deals with stale observations by resetting its dataset each time a condition is met. Our experimental setting exposes the lack of robustness of ET-GP-UCB, since its performance is quite poor on most benchmarks. This is mainly due to the fact that, because the hyperparameters (including the observational noise level σ 2 ) are inferred on the fly, the MLE explains the variance in the observations with an increasingly large observational noise level σ 2 as time goes by. However, by construction of ET-GP-UCB, the greater σ 2 , the less the dataset will be reset. As a consequence, on some benchmarks, the event is never triggered (or not triggered enough) and the performance of ET-GP-UCB is close to (sometimes worse than) the performance of GP-UCB. For examples, refer to Hartmann6d (see Figure 18 ), Shekel (see Figure 16 ) 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 19 ). R-GP-UCB. R-GP-UCB deals with stale data by resetting its dataset (like ET-GP-UCB). It indirectly accounts for temporal correlations by estimating an hyperparameter ϵ, and the reset is triggered each time the dataset size exceeds N (ϵ) (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 12 ), 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 ϵ. Although the temporal covariance function is based on a 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 are 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 8 ), 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 to remove irrelevant observations, W-DBO achieves simultaneously good predictive performance and a 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 variations are more pronounced, the dataset size of W-DBO experiences sudden drops, as can be seen with Ackley (see Figure 13 ), Shekel (see Figure 16 ) or 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 to its ability to adapt in very different contexts, W-DBO outperforms state-of-the-art DBO solutions by a comfortable margin. This performance gap can be seen in its average performance across all benchmarks (see Figure 5 ), but also on most of the benchmarks themselves, e.g., Schwefel (see Figure 9 ), Ackley (see Figure 13 ), Shekel (see Figure 16 ), Hartmann-6 (see Figure 18 ) or Powell (see Figure 19 ).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 spatiotemporal domains of the functions.One of the videos depict the optimization of the Six-Hump Camel function 4 on the domain [-2, 2] 2 . The SHC function is SHC(x, t) = 4 -2.1x 2 + x 4 3 x 2 + xt + -4 + 4t 2 t 2 .To study how W-DBO reacts to sudden changes in the objective function, the other video depict the optimization of the piecewise functionf (x, t) = SHC(x, t) if t < -1 2 , SHC(t, x) otherwise.A snapshot from the latter can be found in Figure 22 . It illustrates that the benefits brought by W-DBO are substantial, since the algorithm is able to track max x∈S f (x, t) over the time t while 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).For the sake of completeness, we explicitly discuss the limitations of W-DBO in this appendix. Four limitations were identified:• As for any BO algorithm, W-DBO exploits a GP as a surrogate model (see Assumption 3.1).If the objective function f cannot be properly approximated by a GP, we expect the performance of W-DBO to decline.• As for any BO algorithm, W-DBO conducts GP inference, which causes it to manipulate inverses of 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 a limitation if too many observations are kept in the dataset.• We also introduce a 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 if the objective function does not meet this assumption.• Finally, W-DBO is not exempt from the effects of 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 ChecklistAt least, not immediately.Its Python documented implementation can be found at https://github.com/WDBO-ALGORITHM/wdbo_ algo. A PyPI package can be quickly installed with the command pip install wdbo-algo.To prevent this, the criterion (12) could be computed on every element of 2 |D| at each iteration. Unfortunately, this policy does not scale well with the dataset size |D|.The video is accessible at https://abardou.github.io/assets/vid/PermSix-Hump_Camel_25.0_ 240.mp4 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: - 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. 
- 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.
- 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.
questions: 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?
- 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: limitations: Yes, discussed in appendices.
paper_summary: paper_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.
strength_weakness: strength_weakness: Typo (Line 305): “to the best of your knowledge” should instead be “to the best of our knowledge”. 

One 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.
questions: 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: limitations: Yes.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

__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? 

Notably, 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.)

In 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.

__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.

__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.

An 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.

__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?

__Benchmarking, part 3__: Please add standard errors to the tables.

__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.

__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.

Minor: 
- "BO inference" is used numerous times in place of "GP inference', which is the more precise description of what is occuring.
- L160: "uses" -> "used"
questions: questions: Questions are apparent from the "Weaknesses" section.
limitations: limitations: Limitations have been addressed.
paper_summary: paper_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.

One 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.

The 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. 

In 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. 

Numerical 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.
strength_weakness: strength_weakness: - 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.

- The proposed approximation of the 2-Wasserstein distance criterion is potentially loose, however, the authors recognize this and the empirical algorithmic performance is strong.

- 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.

[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.

[2] Galy-Fajou, Théo, and Manfred Opper. "Adaptive inducing points selection for gaussian processes." arXiv preprint arXiv:2107.10066 (2021).

[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.

[4] Maus, Natalie, et al. "Approximation-Aware Bayesian Optimization." arXiv preprint arXiv:2406.04308 (2024).
questions: questions: See weaknesses
limitations: limitations: Limitations are discussed at length in Appendix C and H.

Comments: title: Rebuttal Acknowledgement
Response to Further Comments (2/2)
Response to Further Comments (1/2)
Further comments
More Clarifications
Response to Rebuttal
Response to Rebuttal
comment: I would like to thank authors for their response and addressing the points I raised in my review. I am happy to keep my assessment of their work the same.
> Regret is not enough to assess the quality of the point removal metric

In our previous response, we misinterpreted the reviewer’s point about the removal metric. Allow us to address this concern with a more adequate answer.

As far as we understand, the reviewer asks for evidence that the removed points are indeed irrelevant. In the paper, we clearly define the relevancy of an observation $\mathbf x_i$ as the Wasserstein distance between two GP posteriors, one conditioned on the whole dataset $\mathcal{D}$ and the other conditioned on $\tilde{\mathcal{D}} = \mathcal{D} \setminus \\{\mathbf x_i\\}$. All four reviewers agree that this definition is intuitive and reasonable.

We have further illustrated this metric of relevancy in the rebuttal. In fact, in the PDF uploaded with the rebuttal, Figures 2 and 3 bring evidence that the metric is indeed aligned with our intuition. Concretely, it is low (respectively, large) when the two GP posteriors are similar (resp., dissimilar).

A corollary to this definition is that observation relevancy is a continuous non-negative variable; it is zero when the GP posteriors are identical and can be arbitrarily large depending on how the removal of an observation changes the posterior. However, removing an observation from a dataset is a binary decision, since it boils down to labeling each observation as “kept in the dataset” or “removed from the dataset”. To do so, W-DBO (described in Algorithm 1) consumes a removal budget (discussed in Section 4 and Appendix F) by greedily removing the observation with the **smallest value of relevancy** in its dataset, according to the definition of relevancy introduced in the paper.

That is why we were confused when the reviewer asked us for a proof that the observations removed by W-DBO were indeed irrelevant. Because the metric used by W-DBO follows directly from the very definition of observation relevancy, the observations removed by W-DBO are **by definition** the ones minimizing relevancy.

We do not claim that the removal policy advertised in the paper (*i.e.*, the point removal metric, the removal budget and the greedy removal of observations) is the optimal policy to clean a dataset from unwanted observations. Our claim and main result is that a DBO algorithm that uses this policy (that is, W-DBO) outperforms, by a significant margin, state-of-the-art DBO algorithms that either propose to (i) not remove any observation (TV-GP-UCB), (ii) sample the objective function less frequently (ABO) or (iii) reset the dataset according to some condition (R-GP-UCB, ET-GP-UCB). This is discussed and demonstrated at length in Section 5 and in Appendix G.2 and G.3. Note that regret is then the correct metric to support this claim.

> Assess the query removal in isolation (e.g., predictive performance before and after removal compared to previous methods)

Such an assessment has been provided in the PDF uploaded with the rebuttal. Figures 2 and 3 show that our relevancy metric about an observation takes a small (respectively, large) value when the GP posteriors obtained before and after the observation is removed are similar (resp., dissimilar). Because the predictive performance is a distance function applied to the GP posterior and to the objective function, if the two GP posteriors are similar (i.e., the relevancy metric is low), then their predictive performance will be similar as well.

We can make this more explicit by plotting the objective function when we will include Figures 2 and 3 in the camera-ready version.

**In conclusion**

We hope our answers have put the last concerns of the reviewer about potential confounding factors at rest. Indeed, on the one hand, the better performance of W-DBO cannot be explained by implementation quality. On the other hand, we clarified that the main goal of the paper is to introduce a new DBO algorithm that outperforms state-of-the-art DBO algorithms. No claim is made on the optimality of the relevancy metric nor on the observation removal policy, except that it is this metric that gives a decisive advantage to W-DBO.

We remain available to answer further questions from the reviewer before the discussion period ends.

**References**

[1] Balandat, M., Karrer, B., Jiang, D., Daulton, S., Letham, B., Wilson, A. G., & Bakshy, E. (2020). BoTorch: A framework for efficient Monte-Carlo Bayesian optimization. Advances in neural information processing systems, 33, 21524-21538.

[2] Eriksson, D., Pearce, M., Gardner, J., Turner, R. D., & Poloczek, M. (2019). Scalable global optimization via local Bayesian optimization. Advances in neural information processing systems, 32.
Dear Reviewer iNzn,

Thank you for your once again prompt feedback. We further address your points below.

> This metric is heavily impacted by the implementation quality of each algorithm and runtime-impacting hyperparameters such as acquisition optimization budget

Indeed, and that is why we have been very careful to implement all algorithms with the same and well-established BO library (namely BOTorch [1]), as stated in Appendix G.1. The runtime-impacting tasks (e.g., acquisition function optimization, hyperparameters optimization) are all implemented with the same methods (namely `botorch.optim.optimize_acqf` and gradient descent on marginal log-likelihood, respectively) and with the same budget. We can emphasize this more in Appendix G.1 of the camera-ready version.

> Scalable methods such as TurBO [...] should have been part of the initial evaluation

TurBO [2] is a trust region (TR)-based BO algorithm designed to optimize a static objective function. Its core idea is to maintain TRs on a static domain, centered on the best solutions found so far, and iteratively update them until they shrink towards the global maximum of the static objective function.

TurBO sets aside observations but does not discard them in the same way as DBO algorithms (e.g., R-GP-UCB, ET-GP-UCB, W-DBO) do. Indeed, DBO algorithms remove a discarded observation from memory, whereas TurBO simply does not exploit it, but may exploit it at a later stage if one of its TRs is expanded or is reset. More importantly, maintaining TRs in a dynamic environment, which includes both a dynamic domain for the queries (queries in the past cannot be made anymore, queries in the future cannot be made yet) and a global maximum that is constantly moving, requires to develop a substantially different TR-based BO algorithm than TurBO.

Other solutions that could be interesting in the dynamic setting include sparse GPs, as pointed out by Reviewer NVX3. Unfortunately, these solutions were designed to approximate exact GPs on a static objective function, and would also require non-trivial work to be adapted to a dynamic setting. As Reviewer NVX3 puts it, “the sparse GP methodology is relevant but cannot be applied in the current version of the papers and non-trivial work would need to be carried out”.

Following Reviewer NVX3’s remark, we will mention and discuss sparse GPs in the camera-ready version. We will also add TR-based BO algorithms to this discussion.

> Plots with # iterations in the X-axis need to be included

The number of iterations of each algorithm, for each experiment, can be deduced from its average response time given on the x-axes of Figures 6-17: given an experiment duration $d$ and an average response time $r$, the number of iterations is $d / r$. Note that, as stated in Appendix G.1, $d = 600$ seconds in our experimental setting.

In general, W-DBO has a lower response time (*i.e.*, it makes more queries in a given amount of time) than TV-GP-UCB and ABO, but has a larger response time (*i.e.*, it makes fewer queries in a given amount of time) than R-GP-UCB.

Figure 2 shows the running average regret w.r.t. time for the Ackley experiment, and W-DBO’s average regret rapidly becomes significantly lower than the average regret of the other solutions. This value remains lower during the whole experiment, showing the sample efficiency of W-DBO.

Although we are not sure if adding plots with the number of iterations on the x-axis would provide much more additional information, we can add them in a dedicated appendix of the camera-ready version.
__Benchmarking:__ I'm inclined to agree. I'm surprised that the setup in DBO deviates so substantially from the rest of the field, but I digress.

__Evidence for removal:__ Thanks to the authors for clarifying. Unfortunately, I had missed the fact that runtime (and not the number of iterations) was on the x-axis.  Since this metric is heavily impacted by the implementation quality of each algorithm and runtime-impacting hyperparameters such as acquisition optimization budget.

If runtime and throughput is the only relevant metric, scalable methods such as TurBO, which naturally discards observations, should have been part of the initial evaluation. If sample-efficiency is still important, plots with #iterations in the X-axis need to be included, at least in the appendix. 

As I hope the authors see, there are still substantial confounding factors that makes it very difficult to assess if the removal itself is actually good. Regret over time is a downstream metric, and so that by itself is used to assess the quality of the point removal metric, there can be just about no confounding factors. Otherwise, experiments that assess the query removal _in isolation_ (e.g., predictive performance before and after removal compared to previous methods) should certainly be included.
Dear Reviewer iNzn,

Thank you for your feedback on our response. We further address your points below.

**Benchmarking**

We agree with the reviewer that this discussion boils down to an alleged unfairness issue in the benchmarking. We think that this issue has been put to rest in three ways: (i) we have shown that most benchmarking in the *DBO* literature was conducted with a non-ARD kernel, (ii) we have made it clear that, for each experiment in the paper, each DBO solution uses a non-ARD kernel and (iii) we have run additional experiments with an ARD kernel to make sure that W-DBO was still the best-performing solution with ARD kernels.

**Evidence for removal**

An instance of W-DBO that does not remove any point would have a behavior very close to the behavior of TV-GP-UCB. In fact, when it does not remove any point, W-DBO is simply TV-GP-UCB with a different temporal kernel. This suggested instance of W-DBO would therefore be quite redundant with TV-GP-UCB, and would not exhibit superior performance for the same reason TV-GP-UCB does not.

As time goes by, TV-GP-UCB and similar solutions that do not remove irrelevant observations (e.g., ABO and the instance of W-DBO suggested by the reviewer) experience a larger and larger GP inference time because irrelevant observations accumulate in their datasets. This increasingly large inference time prevents these solutions from querying the objective function as often as W-DBO, which regularly removes observations that are deemed irrelevant according to the definition in the paper (which all four reviewers find intuitive and reasonable). Because we benchmark the solutions in a continuous time setting, failing to query the objective function often enough ultimately hinders the optimization of the dynamic objective function by limiting how much two consecutive observations can be correlated. In other words, W-DBO outperforms TV-GP-UCB and ABO precisely because it removes some observations. At the same time, recall that W-DBO also outperforms R-GP-UCB and ET-GP-UCB, two solutions that also remove observations but in a less refined way.

This leads us to make two conclusions: (i) a version of W-DBO that does not remove any irrelevant observation would **not** be the gold standard suggested by the reviewer and (ii) the very fact that W-DBO exhibits better performance than all the other DBO solutions (ABO and TV-GP-UCB on the one hand, R-GP-UCB and ET-GP-UCB on the other hand) is a strong evidence that the observations that it removes are indeed irrelevant, or at least, irrelevant enough so that removing them is a better strategy than keeping them in the dataset.

Finally, we acknowledge that our paper opens multiple interesting research questions about observation relevancy, and does not provide all the answers to these questions. In our opinion, this is more a strength than a weakness. W-DBO is an original answer to the DBO problem that outperforms the existing solutions while opening various research questions, which is why we believe our paper is of interest to the BO community.

We hope our answers have lifted any remaining misunderstanding.
Thanks to the authors for their response.

__Benchmarking:__ Thanks to the authors for the clarification. Ultimately, while ARD is clearly the convention in BO (I cannot speak for DBO specifically, so I thank the authors for clarifying) the issue is ultimately one of fairness in evaluation. As mentioned previously, ARD kernels will be disadvantaged on symmetric functions. I was under the impression that ABO ran an ARD kernel (their paper seems to suggest it), but if the authors say that all methods run a non-ARD kernel, I will take their word for it.


__Evidence for removal:__ There has not been evidence that suggests that the removal criterion is actually efficient. One would expect there to be a comparison between not removing at all (should be ideal) and removing stale observations (should be close to ideal if the method is potent). As stated previously, regret performance as a measure of this is not convincing enough, in my view. The ablation plot on $\alpha$ is a start, but running the D-BO without removing any observations (and an otherwise identical setup to W-DBO) as a gold standard is, in my view, a must have for all experiments.

__Wasserstein Distance figure:__ I think this is a good addition. I would also consider adding where the removed observation was located.
______
Given unconvincing experimental evidence that the proposed method does what it seeks to do (namely remove unimportant observations), otherwise unconvincing experimental results, and a presentation that could use some polish, I will maintain my rating.
Thanks for the response. In particular, I see I have misunderstood how the experimental case studies were defined (see authors' global response), and so this stated weakness is largely alleviated. I have raised my score correspondingly, and look forward to the authors' stated clarifications in the camera-ready (or otherwise future) version.
Thanks for your response. I agree that the sparse GP methodology is relevant but cannot be applied in the current version of the papers and non-trivial work would need to be carried out. Adding a discussion of these challenges and the connections with the current problem would be sufficient. 

Thanks further for the clarification about the synthetic benchmarks being dynamic, I do agree it should be made clearer.

I remain positive about the work, and keep my recommendation of acceptance.