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

Anthony Bardou
IC, EPFL
Lausanne, Switzerland
<EMAIL_ADDRESS>
Patrick Thiran
IC, EPFL
Lausanne, Switzerland
<EMAIL_ADDRESS>
Giovanni Ranieri
IC, EPFL
Lausanne, Switzerland
<EMAIL_ADDRESS>
###### 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.

##  1 Introduction

Many real-world problems require the optimization of a costly-to-evaluate
objective function $f:\mathcal{S}\subseteq\mathbb{R}^{d}\to\mathbb{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. 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 traditional 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:\mathcal{S}\times\mathcal{T}\to\mathbb{R}$ is a time-
varying, black-box, costly-to-evaluate objective function with spatial domain
$\mathcal{S}\subseteq\mathbb{R}^{d}$ and temporal domain
$\mathcal{T}\subseteq\mathbb{R}$. Unlike common black-box objective functions
that take as input a point $\bm{x}\in\mathcal{S}$ in a spatial domain only,
the function $f$ takes as input a point
$(\bm{x},t)\in\mathcal{S}\times\mathcal{T}$ in space-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 $\mathcal{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 $\bm{x}\in\mathcal{S}$ but cannot query past (i.e.,
$t<t_{0}$) nor future111At least, not immediately. (i.e., $t>t_{0}$) points in
time, (ii) as time (only) moves forward, a previously collected observation
$\left((\bm{x},t),f(\bm{x},t)\right)$ 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 algorithm’s sampling frequency 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 turns into 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 Background

Let us start by describing BO with a GP prior, as introduced by [16]. Given an
objective function $f:\mathcal{S}\subseteq\mathbb{R}^{d}\to\mathbb{R}$, BO
assumes that $f$ is a $\mathcal{GP}\left(0,k(\bm{x},\bm{x}^{\prime})\right)$.
For any $\bm{x}\in\mathcal{S}$, and given a dataset of observations
$\mathcal{D}=\left\\{\left(\bm{x}_{i},y_{i}\right)\right\\}_{i\in\llbracket
1,n\rrbracket}$, where $\bm{x}_{i}\in\mathcal{S}$ is a previously queried
input with (noisy) function value $y_{i}=f(\bm{x}_{i})+\epsilon$,
$\epsilon\sim\mathcal{N}\left(0,\sigma^{2}\right)$, the posterior distribution
of $f(\bm{x})$ is $\mathcal{N}(\mu(\bm{x}),\sigma^{2}(\bm{x}))$, where

$\displaystyle\mu(\bm{x})$
$\displaystyle=\bm{k}^{\top}(\bm{x},\bm{X})\bm{\Delta}^{-1}\bm{y},$ (1)
$\displaystyle\sigma^{2}(\bm{x})$
$\displaystyle=k(\bm{x},\bm{x})-\bm{k}^{\top}(\bm{x},\bm{X})\bm{\Delta}^{-1}\bm{k}(\bm{x},\bm{X})$
(2)

with $\bm{X}=\left(\bm{x}_{1},\cdots,\bm{x}_{n}\right)$,
$\bm{y}=\left(y_{1},\cdots,y_{n}\right)$,
$\bm{\Delta}=\bm{k}(\bm{X},\bm{X})+\sigma^{2}\bm{I}$ and
$\bm{k}(\mathcal{X},\mathcal{Y})=\left(k\left(\bm{x}_{i},\bm{x}_{j}\right)\right)_{\bm{x}_{i}\in\mathcal{X},\bm{x}_{j}\in\mathcal{Y}}$.

To find $\bm{x}_{n+1}$, the input to query at the ($n+1$)th iteration, a BO
algorithm exploits an acquisition function $\varphi:\mathcal{S}\to\mathbb{R}$.
The acquisition function $\varphi$ quantifies the benefits of querying the
input $\bm{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
$\bm{x}_{n+1}=\operatorname*{arg\,max}_{\bm{x}\in\mathcal{S}}\varphi(\bm{x})$.

BO inference extends naturally to dynamic problems, by simply 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
$\mathcal{D}=\left\\{\left(\left(\bm{x}_{i},t_{i}\right),y_{i}\right)\right\\}_{i\in\llbracket
1,n\rrbracket}$, where
$\left(\bm{x}_{i},t_{i}\right)\in\mathcal{S}\times\mathcal{T}$ is a previously
queried input with the (noisy) function value
$y_{i}=f(\bm{x}_{i},t_{i})+\epsilon$,
$\epsilon\sim\mathcal{N}\left(0,\sigma^{2}\right)$, the posterior distribution
of $f(\bm{x},t)$ is
$\mathcal{N}\left(\mu(\bm{x},t),\sigma^{2}(\bm{x},t)\right)$ for any
$(\bm{x},t)\in\mathcal{S}\times\mathcal{T}$, with

$\displaystyle\mu(\bm{x},t)$
$\displaystyle=\bm{k}^{\top}(\left(\bm{x},t\right),\bm{X})\bm{\Delta}^{-1}\bm{y},$
(3) $\displaystyle\sigma^{2}(\bm{x},t)$
$\displaystyle=k(\left(\bm{x},t\right),\left(\bm{x},t\right))-\bm{k}^{\top}(\left(\bm{x},t\right),\bm{X})\bm{\Delta}^{-1}\bm{k}(\left(\bm{x},t\right),\bm{X})$
(4)

with
$\bm{X}=\left(\left(\bm{x}_{1},t_{1}\right),\cdots,\left(\bm{x}_{n},t_{n}\right)\right)$,
$\bm{y}=\left(y_{1},\cdots,y_{n}\right)$,
$\bm{\Delta}=\bm{k}(\bm{X},\bm{X})+\sigma^{2}\bm{I}$ and
$\bm{k}(\mathcal{X},\mathcal{Y})=\left(k\left(\left(\bm{x}_{i},t_{i}\right),\left(\bm{x}_{j},t_{j}\right)\right)\right)_{\left(\bm{x}_{i},t_{i}\right)\in\mathcal{X},\left(\bm{x}_{j},t_{j}\right)\in\mathcal{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
$\left(\bm{x}_{n+1},t_{0}\right)$, with
$\bm{x}_{n+1}=\operatorname*{arg\,max}_{\bm{x}\in\mathcal{S}}\varphi(\bm{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}+\delta_{t}]$.
In that case, the next queried input is
$\left(\bm{x}_{n+1},t_{n+1}\right)=\operatorname*{arg\,max}_{(\bm{x},t)\in\mathcal{S}\times[t_{0},t_{0}+\delta_{t}]}\varphi(\bm{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 $\epsilon\in[0,1]$. On the
one hand, the authors propose R-GP-UCB, which handles data staleness by
resetting the dataset every $N(\epsilon)$ iterations. On the other hand, the
authors also propose TV-GP-UCB that incorporates data staleness by weighing
the covariance of two queries $\bm{q}_{i}=(\bm{x}_{i},t_{i})$ and
$\bm{q}_{j}=(\bm{x}_{j},t_{j})$ by $(1-\epsilon)^{|i-j|/2}$. In [22], the
authors use the same model with an event-triggered erasure of the entire
dataset to deal with stale data. 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 work under the assumption 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 a DBO algorithm by
allowing the use of proof techniques similar to the ones used for classical
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, since BO inference is in $\mathcal{O}(n^{3})$, with $n$
the dataset size. Although it is reasonable to ignore this for classical BO
tasks 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 collected observations in their datasets,
whereas in practice, their response times would asymptotically become
prohibitive, especially for DBO tasks with infinite time horizon. 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 static during the optimization. However, not estimating
hyperparameters on the fly is unrealistic, 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((\bm{x},t),(\bm{x}^{\prime},t^{\prime}))=k_{S}(||\bm{x}-\bm{x}^{\prime}||_{2})k_{T}(|t-t^{\prime}|)$
to accurately model complex spatio-temporal correlations and sample the
objective function only when deemed necessary. Although this reduces the size
of ABO’s dataset, the authors do not propose a way to remove stale
observations, only to add less frequently new observations. Therefore, using
ABO will still become prohibitive in the long run.

As far as we know, there is no method in the DBO literature able to quantify
the relevancy of an observation in its dataset 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.

##  3 A Wasserstein Distance-Based Criterion

###  3.1 Core Assumptions

To address the DBO problem under suitable smoothness conditions, let us make
the usual assumption of BO, using a Gaussian Process (GP) as a surrogate model
for $f$ (see [16]).

######  Assumption 3.1.

$f$ is a
$\mathcal{GP}\left(0,k(\left(\bm{x},t\right),\left(\bm{x}^{\prime},t^{\prime}\right)\right)$,
whose mean is $0$ (without loss of generality) and whose covariance function
is denoted by
$k:\mathcal{S}\times\mathcal{T}\times\mathcal{S}\times\mathcal{T}\to\mathbb{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.

$k((\bm{x},t),(\bm{x}^{\prime},t^{\prime}))=\lambda
k_{S}(||\bm{x}-\bm{x}^{\prime}||_{2},l_{S})k_{T}(|t-t^{\prime}|,l_{T}),$ (5)

with $\lambda>0$, $k_{S}:\mathbb{R}_{+}\to[0,1]$ and
$k_{T}:\mathbb{R}_{+}\to[0,1]$ two positive correlation functions,
parameterized by lengthscales $l_{S}>0$ and $l_{T}>0$, respectively. The
factor $\lambda>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 [25] 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.

###  3.2 Measuring the Relevancy of an Observation

By definition, when an irrelevant observation gets removed from the dataset,
the GP future predictions experience hardly any change. Therefore, we propose
to measure the relevancy of an observation
$\bm{o}_{i}=((\bm{x}_{i},t_{i}),y_{i})$ by measuring the impact that the
removal of $\bm{o}_{i}$ has on the GP future predictions.

Let $\mathcal{GP}_{\mathcal{D}}$ be the GP conditioned on the dataset
$\mathcal{D}=\left\\{\left((\bm{x}_{i},t_{i}),y_{i}\right)\right\\}_{i\in\llbracket
1,n\rrbracket}$, with $(\bm{x}_{i},t_{i})\in\mathcal{S}\times\mathcal{T}$ and
$y_{i}=f(\bm{x}_{i},t_{i})+\epsilon,\epsilon\sim\mathcal{N}\left(0,\sigma^{2}\right)$.
Without loss of generality, let us measure the impact of removing
$((\bm{x}_{1},t_{1}),y_{1})$ from the dataset on the GP future predictions.
Clearly, the domain of future predictions at time $t_{0}$, denoted by
$\mathcal{F}_{t_{0}}$, must include the whole space $\mathcal{S}$ and only the
future time interval $[t_{0},+\infty)$:

$\mathcal{F}_{t_{0}}=\mathcal{S}\times[t_{0},+\infty).$ (6)

We compare a GP conditioned on the whole dataset, denoted by
$\mathcal{GP}_{\mathcal{D}}$, with a GP conditioned on $\tilde{\mathcal{D}}$,
the dataset without $(\bm{x}_{1},t_{1},y_{1})$, denoted by
$\mathcal{GP}_{\tilde{\mathcal{D}}}$. For an arbitrary point
$(\bm{x},t)\in\mathcal{F}_{t_{0}}$, $\mathcal{GP}_{\mathcal{D}}$ provides a
posterior distribution
$\mathcal{N}_{\mathcal{D}}(\bm{x},t)=\mathcal{N}\left(\mu_{\mathcal{D}}(\bm{x},t),\sigma^{2}_{\mathcal{D}}(\bm{x},t)\right)$,
and so does $\mathcal{GP}_{\tilde{\mathcal{D}}}$ with
$\mathcal{N}_{\tilde{\mathcal{D}}}(\bm{x},t)=\mathcal{N}\left(\mu_{\tilde{\mathcal{D}}}(\bm{x},t),\sigma^{2}_{\tilde{\mathcal{D}}}(\bm{x},t)\right)$.
We compare these two distributions by using the 2-Wasserstein distance [26],
given by

$W_{2}\left(\mathcal{N}_{\mathcal{D}}(\bm{x},t),\mathcal{N}_{\tilde{\mathcal{D}}}(\bm{x},t)\right)=\left(\left(\mu_{\mathcal{D}}(\bm{x},t)-\mu_{\tilde{\mathcal{D}}}(\bm{x},t)\right)^{2}+\left(\sigma_{\mathcal{D}}(\bm{x},t)-\sigma_{\tilde{\mathcal{D}}}(\bm{x},t)\right)^{2}\right)^{\frac{1}{2}}.$
(7)

A natural extension of the 2-Wasserstein distance from a point
$(\bm{x},t)\in\mathcal{F}_{t_{0}}$ to the entire domain $\mathcal{F}_{t_{0}}$
is

$W_{2}\left(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\tilde{\mathcal{D}}}\right)=\left(\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}W_{2}^{2}\left(\mathcal{N}_{\mathcal{D}}(\bm{x},t),\mathcal{N}_{\tilde{\mathcal{D}}}(\bm{x},t)\right)d\bm{x}dt\right)^{\frac{1}{2}}.$
(8)

Observe that (8) is a criterion that effectively captures the impact of
removing the observation $\bm{o}_{1}=((\bm{x}_{1},t_{1}),y_{1})$ from the
dataset on the domain of future predictions $\mathcal{F}_{t_{0}}$. However, as
illustrated in Appendix E, the covariance function hyperparameters
$\bm{\theta}=(\lambda,l_{S},l_{T})$ control the magnitude of (8). As a result,
(8) cannot be directly uses as a gauge of the relevancy of an observation. In
fact, the very same quantity computed by (8) might indicate a much relevant
observation when the hyperparameters in $\bm{\theta}$ are small enough, but
also an irrelevant observation when the hyperparameters in $\bm{\theta}$ are
large enough.

To remove this ambiguity, we normalize (8) by
$W_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset})$ (i.e., the
2-Wasserstein distance between the GP conditioned on $\mathcal{D}$ and the
prior GP), and we obtain the ratio

$R(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\tilde{\mathcal{D}}})=\frac{W_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\tilde{\mathcal{D}}})}{W_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset})}.$
(9)

Intuitively, $W_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset})$
measures the impact of resetting the whole dataset $\mathcal{D}$ on the future
predictions and serves as a baseline that puts into perspective the distance
measured by (8). Technically, taking the ratio (9) successfully cancels out
the the influence of the covariance function hyperparameters on the magnitude
of the Wasserstein distances, as further discussed in Appendix E.

Exploiting this criterion is straightforward. When (9) is low, one can infer
that $\bm{o}_{1}$ can be safely removed from the dataset since it will have
virtually no impact on future predictions. Conversely, when (9) is large, one
can infer that removing $\bm{o}_{1}$ would alter the future predictions 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

###  4.1 Computational Tractability

In [27], 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 the computational cost of an iteration must crucially be kept as low as
possible to ensure a high sampling frequency. In this section, we put this
issue to rest, by deriving explicit formulas for approximating (9). These
formulas are computationally cheap enough to be exploited on the fly by a DBO
algorithm.

In Appendix A, we show that (7) can be computed efficiently. Next, in Appendix
B, we apply these results to obtain a closed-form expression of an upper bound
of (8). The key observation for deriving this result is to approximate the
integrals in (8) with a product of convolutions of the covariance functions
with themselves in space and time. The same trick can be used for
approximating $W_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset})$.

######  Theorem 4.1.

Let $t_{0}$ be the present time and
$\mathcal{D}=\left\\{((\bm{x}_{i},t_{i}),y_{i})\right\\}_{i\in\llbracket
1,n\rrbracket}$ be a dataset of observations made before $t_{0}$. Let
$\tilde{\mathcal{D}}=\left\\{((\bm{x}_{i},t_{i}),y_{i})\right\\}_{i\in\llbracket
2,n\rrbracket}$ and $\mathcal{F}_{t_{0}}=\mathcal{S}\times[t_{0},+\infty)$ be
the domain of future predictions. Then, an upper bound for
$W_{2}^{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\tilde{\mathcal{D}}})$ on
$\mathcal{F}_{t_{0}}$ is

$\begin{split}\hat{W}_{2}^{2}\left(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\tilde{\mathcal{D}}}\right)=\lambda^{2}(a^{2}+\bm{E})C((\bm{x}_{1},t_{1}),&(\bm{x}_{1},t_{1}))+\lambda^{2}(2a\bm{b}+\bm{c})C((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}})\\\
&+\lambda^{2}\bm{1}^{\top}\left(\left(\bm{b}\bm{b}^{\top}+\bm{M}\right)\odot
C(\tilde{\mathcal{D}},\tilde{\mathcal{D}})\right)\bm{1}\end{split}$ (10)

with
$C(\mathcal{X},\mathcal{Y})=\left((k_{S}*k_{S})(\bm{x}_{j}-\bm{x}_{i})\cdot(k_{T}*k_{T})_{t_{0}-t_{i}}^{+\infty}(t_{j}-t_{i})\right)_{\begin{subarray}{c}(\bm{x}_{i},t_{i})\in\mathcal{X}\\\
(\bm{x}_{j},t_{j})\in\mathcal{Y}\end{subarray}}$, $\bm{1}$ the conformable
vector of ones, $\odot$ the Hadamard product, $(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]$. The terms $a$, $\bm{b}$, $\bm{c}$, $\bm{E}$ and
$\bm{M}$ are explicited in Appendices A and B.

Moreover, an upper bound for
$W_{2}^{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset})$ on
$\mathcal{F}_{t_{0}}$ is

$\hat{W}_{2}^{2}\left(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset}\right)=\lambda^{2}\left(\bm{y}^{\top}\bm{\Delta}^{-\top}C\left(\mathcal{D},\mathcal{D}\right)\bm{\Delta}^{-1}\bm{y}+\bm{1}^{\top}\left(\bm{\Delta}^{-1}\odot
C\left(\mathcal{D},\mathcal{D}\right)\right)\bm{1}\right).$ (11)

This theorem provides the analytic form of an upper bound for the Wasserstein
distance between $\mathcal{GP}_{\mathcal{D}}$ and
$\mathcal{GP}_{\tilde{\mathcal{D}}}$ and the Wasserstein distance between
$\mathcal{GP}_{\mathcal{D}}$ and $\mathcal{GP}_{\emptyset}$ on the domain of
future predictions $\mathcal{F}_{t_{0}}$. Using it, we can compute an
approximation $\hat{R}$ of the relative criterion (9), that is

$\hat{R}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\tilde{\mathcal{D}}})=\frac{\hat{W}_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\tilde{\mathcal{D}}})}{\hat{W}_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset})}.$
(12)

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 product of convolutions of the covariance
functions $k_{S}$ and $k_{T}$ with themselves. 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).

Table 1: Usual covariance functions. $\Gamma$ is the Gamma function and $K_{\nu}$ is a modified Bessel function of the second kind of order $\nu$. Covariance Function | $k(\bm{x})$
---|---
Squared-Exponential ($l$) | $e^{-\frac{||\bm{x}||_{2}^{2}}{2l^{2}}}$
Matérn ($\nu,l$) | $\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}||\bm{x}||_{2}}{l_{S}}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}||\bm{x}||_{2}}{l_{S}}\right)$

In practice, the upper bounds (10) and (11) can only be computed efficiently
if the convolutions of the spatial and the temporal covariance functions can
themselves be computed efficiently. The analytic forms for the convolutions of
two usual covariance functions listed in Table 1, namely the Squared-
Exponential covariance function and the Matérn covariance function [28], are
provided in Appendix D together with Tables 3 and 4 that list the convolutions
for the spatial and the temporal covariance functions with themselves,
respectively. The details for their computation are also provided in this
appendix. In a nutshell, the convolution formulas are obtained first in the
Fourier domain by computing the square of the spectral densities of the
covariance functions, and next by computing the inverse Fourier transform of
the result.

Together, Tables 3 and 4 in Appendix D as well as Theorem 4.1 show that the
approximation of (9) given by (12) can be computed efficiently in an online
setting. We now leverage (12) to propose a DBO algorithm able to pinpoint and
remove irrelevant observations in its dataset.

###  4.2 W-DBO

The metric (9) and its approximation (12) can be seen as a relative error (or
drift), expressed as a percentage, that separates $\mathcal{GP}_{\mathcal{D}}$
and $\mathcal{GP}_{\tilde{\mathcal{D}}}$. Indeed, the Wasserstein distance
$W\left(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\tilde{\mathcal{D}}}\right)$
is scaled by the Wasserstein distance
$W\left(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset}\right)$, that is,
the distance between $\mathcal{GP}_{\mathcal{D}}$ and the prior. In other
words, (9) and its approximation (12) measure the relative drift from
$\mathcal{GP}_{\mathcal{D}}$ to $\mathcal{GP}_{\tilde{\mathcal{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
$\mathcal{GP}_{\tilde{\mathcal{D}}}$ drift exponentially fast from
$\mathcal{GP}_{\mathcal{D}}$. To keep this exponential drift under control,
one can use a removal budget $b(t)=(1+\alpha)^{t}$ that allows a maximal
relative drift from $\mathcal{GP}_{\mathcal{D}}$ of $\alpha$ per unit of time
(e.g., if $\alpha=0.1$, the allowed maximal drift is 10 % per unit of time).
Obviously, the cost of removing an observation is given by (12).

Algorithm 1 W-DBO

1: Input: $\mathcal{GP}$, acquisition function $\varphi$, hyperparameter
$\alpha$, clock $\mathcal{C}$, initial observations $\mathcal{D}$.

2: Get current time $t$ from clock $\mathcal{C}$

3: $b=1$

4: while true do

5: Find $\bm{x}_{t}=\operatorname*{arg\,max}_{\bm{x}}\varphi(\bm{x},t)$

6: Observe $y_{t}=f(\bm{x}_{t},t)+\epsilon$,
$\epsilon\sim\mathcal{N}(0,\sigma^{2})$

7: $\mathcal{D}=\mathcal{D}\cup\left\\{(\bm{x}_{t},t,y_{t})\right\\}$

8: Condition $\mathcal{GP}$ on $\mathcal{D}$ and get the MLE parameters
$(\lambda,l_{S},l_{T},\hat{\sigma}^{2})$

9: while $b>1$ do

10: for all $(\bm{x}_{i},t_{i},y_{i})\in\mathcal{D}$ [in parallel] do

11:
$\tilde{\mathcal{D}}=\mathcal{D}\setminus\left\\{(\bm{x}_{i},t_{i},y_{i})\right\\}$

12: Compute $\hat{R}_{i}$ using (12)

13: end for

14: $i^{*}=\operatorname*{arg\,min}_{i\in\llbracket
1,|\mathcal{D}|\rrbracket}\hat{R}_{i}$

15: if $b>1+\hat{R}_{i^{*}}$ then

16:
$\mathcal{D}=\mathcal{D}\setminus\left\\{(\bm{x}_{i^{*}},t_{i^{*}},y_{i^{*}})\right\\}$

17: $b=b/(1+\hat{R}_{i^{*}})$

18: else

19: break

20: end if

21: end while

22: $t^{\prime}=t$

23: Get current time $t$ from clock $\mathcal{C}$

24: $b=b(1+\alpha)^{(t-t^{\prime})/l_{T}}$

25: end while

Algorithm 1 describes W-DBO, a DBO algorithm that shows how to exploit (12) to
remove irrelevant observations on the fly. As described above, the removal
budget is controlled by a single hyperparameter $\alpha$ 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-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-17). This process is repeated until all the budget is consumed.
Observe that, in the exponent used for the removal budget, the duration
$t-t^{\prime}$ is scaled by the temporal lengthscale $l_{T}$. We discuss and
motivate this specific expression of the removal budget in Appendix F. The
sensitivity analysis conducted in Section 5.1 supports this form of the
removal budget $b$, by showing that the same value of the hyperparameter
$\alpha$ 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 $\mathcal{GP}$ and any acquisition function $\varphi$.
Similarly, observe that lines 5-8 in Algorithm 1 describe the usual BO
optimization loop, without any modification.

##  5 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 parameters and the noise level are estimated on the
fly. Please refer to Appendix G.1 for further experimental details.

###  5.1 Sensitivity Analysis

We start by studying the impact of the W-DBO hyperparameter $\alpha$ 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).

Figure 1: (Left) Sensitivity analysis on the Eggholder function. (Right)
Aggregation of sensitivity analyses of W-DBO made on 10 synthetic functions
and a real-world experiment. For aggregation purposes, the average regrets in
each experiment have been normalized between 0 (lowest average regret) and 1
(largest average regret). The average performance of W-DBO over all the
experiments is shown in black. Standard errors are depicted with colored bars
(left) and shaded areas (right).

To study this trade-off, we compare the performance of W-DBO with several
values of its hyperparameter $\alpha$, as illustrated by the left side of
Figure 1. We apply this protocol on 11 different benchmarks (described in
Appendix G.2). The aggregated results (see the right side of Figure 1) 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 $\alpha=\frac{1}{4}$. This hyperparameter value is
used to evaluate W-DBO in the next section.

###  5.2 Comparison with Baselines

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. Experiment ($d+1$) | GP-UCB | R-GP-UCB | TV-GP-UCB | ET-GP-UCB | ABO | W-DBO
---|---|---|---|---|---|---
Rastrigin ($5$) | $\bm{17.81}$ | $53.67$ | $26.50$ | $\underline{19.54}$ | $36.16$ | $\underline{18.54}$
Schwefel ($4$) | $469.10$ | $954.03$ | $520.40$ | $428.97$ | $662.32$ | $\bm{290.34}$
StyblinskiTang ($4$) | $18.83$ | $45.82$ | $\underline{15.74}$ | $22.16$ | $58.40$ | $\bm{13.04}$
Eggholder ($2$) | $542.53$ | $273.60$ | $287.01$ | $559.61$ | $256.92$ | $\bm{225.68}$
Ackley ($4$) | $4.10$ | $4.45$ | $3.27$ | $3.96$ | $3.63$ | $\bm{2.24}$
Rosenbrock ($3$) | $31.37$ | $25.99$ | $17.55$ | $28.79$ | $171.04$ | $\bm{3.81}$
Shekel ($4$) | $2.56$ | $2.21$ | $2.03$ | $2.70$ | $2.06$ | $\bm{1.72}$
Hartmann3 ($3$) | $1.17$ | $\bm{0.26}$ | $0.82$ | $1.06$ | $0.55$ | $0.35$
Hartmann6 ($6$) | $1.33$ | $1.25$ | $0.44$ | $1.46$ | $0.61$ | $\bm{0.32}$
Powell ($4$) | $1992.1$ | $1167.6$ | $1223.4$ | $534.2$ | $9888.6$ | $\bm{428.1}$
Temperature ($3$) | $1.02$ | $\underline{0.69}$ | $1.36$ | $1.25$ | $1.21$ | $\bm{0.68}$
WLAN ($5$) | $1.46$ | $4.84$ | $1.33$ | $4.98$ | $12.94$ | $\bm{1.19}$
Avg. Perf. | $0.48$ | $0.47$ | $0.29$ | $0.54$ | $0.62$ | $\bm{0.01}$

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
G.2, and the performance of each DBO solution is discussed at length in
Appendix G.3. Furthermore, the provided supplementary animated visualizations
are discussed in Appendix G.4.

Figure 2: (Left) Average regrets of the DBO solutions during the optimization
of the Ackley synthetic function. (Right) Dataset sizes of the DBO solutions
during the optimization of the Ackley function.

In this section. we only depict the performance of the DBO solutions on the
Ackley synthetic function in Figure 2, 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 large hole at the center of its
domain. Observe that most DBO solutions miss that hole, as their average
regrets skyrocket between 200 and 400 seconds. On the other hand, W-DBO
manages to quickly exploit the hole at the center of the function’s domain,
thus 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 2). 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
experiences the flatter outer region of the Ackley function again.
Consequently, its dataset size increases again.

Figure 3: Visual summary of the results reported in Table 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.

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

##  6 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 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 your 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.,
computational biology, 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.

## References

  * [1] James Bergstra, Daniel Yamins, and David Cox.  Making a science of model search: Hyperparameter optimization in hundreds of dimensions for vision architectures.  In International conference on machine learning, pages 115–123. PMLR, 2013.
  * [2] Gregory Hornby, Al Globus, Derek Linden, and Jason Lohn.  Automated antenna design with evolutionary algorithms.  In American Institute of Aeronautics and Astronautics, 2006.
  * [3] Anthony Bardou and Thomas Begin.  Inspire: Distributed bayesian optimization for improving spatial reuse in dense wlans.  In Proceedings of the 25th International ACM Conference on Modeling Analysis and Simulation of Wireless and Mobile Systems, pages 133–142, 2022.
  * [4] Daniel Lizotte, Tao Wang, Michael Bowling, and Dale Schuurmans.  Automatic gait optimization with gaussian process regression.  In Proceedings of the 20th International Joint Conference on Artifical Intelligence, IJCAI’07, page 944–949, San Francisco, CA, USA, 2007. Morgan Kaufmann Publishers Inc.
  * [5] Javier González, Joseph Longworth, David C. James, and Neil D. Lawrence.  Bayesian optimization for synthetic gene design.  In NIPS Workshop on Bayesian Optimization in Academia and Industry, 2014.
  * [6] David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek.  Scalable global optimization via local bayesian optimization.  In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [7] Johannes Kirschner, Mojmir Mutny, Nicole Hiller, Rasmus Ischebeck, and Andreas Krause.  Adaptive and safe bayesian optimization in high dimensions via one-dimensional subspaces.  In International Conference on Machine Learning, pages 3429–3438. PMLR, 2019.
  * [8] Anthony Bardou, Patrick Thiran, and Thomas Begin.  Relaxing the additivity constraints in decentralized no-regret high-dimensional bayesian optimization.  arXiv preprint arXiv:2305.19838, 2023.
  * [9] Javier González, Zhenwen Dai, Philipp Hennig, and Neil Lawrence.  Batch bayesian optimization via local penalization.  In Artificial intelligence and statistics, pages 648–657. PMLR, 2016.
  * [10] Samuel Daulton, David Eriksson, Maximilian Balandat, and Eytan Bakshy.  Multi-objective bayesian optimization over high-dimensional search spaces.  In Uncertainty in Artificial Intelligence, pages 507–517. PMLR, 2022.
  * [11] Seokhyun Kim, Kimin Lee, Yeonkeun Kim, Jinwoo Shin, Seungwon Shin, and Song Chong.  Dynamic control for on-demand interference-managed wlan infrastructures.  IEEE/ACM Transactions on Networking, 28(1):84–97, 2019.
  * [12] Aurelio G Melo, Milena F Pinto, Andre LM Marcato, Leonardo M Honório, and Fabrício O Coelho.  Dynamic optimization and heuristics based online coverage path planning in 3d environment for uavs.  Sensors, 21(4):1108, 2021.
  * [13] Doyen Sahoo, Quang Pham, Jing Lu, and Steven CH Hoi.  Online deep learning: learning deep neural networks on the fly.  In Proceedings of the 27th International Joint Conference on Artificial Intelligence, pages 2660–2666, 2018.
  * [14] Charu C Aggarwal, Jiawei Han, Jianyong Wang, and Philip S Yu.  A framework for projected clustering of high dimensional data streams.  In Proceedings of the Thirtieth international conference on Very large data bases-Volume 30, pages 852–863, 2004.
  * [15] Xiao Mingming, Zhang Jun, Cai Kaiquan, Cao Xianbin, and Tang Ke.  Cooperative co-evolution with weighted random grouping for large-scale crossing waypoints locating in air route network.  In 2011 IEEE 23rd International Conference on Tools with Artificial Intelligence, pages 215–222. IEEE, 2011.
  * [16] Christopher Williams and Carl Rasmussen.  Gaussian processes for regression.  Advances in neural information processing systems, 8, 1995.
  * [17] Niranjan Srinivas, Andreas Krause, Sham M. Kakade, and Matthias W. Seeger.  Information-theoretic regret bounds for gaussian process optimization in the bandit setting.  IEEE Transactions on Information Theory, 58(5):3250–3265, 2012.
  * [18] Jonas Mockus.  Application of bayesian approach to numerical methods of global and stochastic optimization.  Journal of Global Optimization, 4:347–365, 1994.
  * [19] Donald R. Jones, Matthias Schonlau, and William J. Welch.  Efficient global optimization of expensive black-box functions.  Journal of Global optimization, 13(4):455–492, 1998.
  * [20] Favour M Nyikosa, Michael A Osborne, and Stephen J Roberts.  Bayesian optimization for dynamic problems.  arXiv preprint arXiv:1803.03432, 2018.
  * [21] Ilija Bogunovic, Jonathan Scarlett, and Volkan Cevher.  Time-varying gaussian process bandit optimization.  In Artificial Intelligence and Statistics, pages 314–323. PMLR, 2016.
  * [22] Paul Brunzema, Alexander von Rohr, Friedrich Solowjow, and Sebastian Trimpe.  Event-triggered time-varying bayesian optimization.  arXiv preprint arXiv:2208.10790, 2022.
  * [23] Xingyu Zhou and Ness Shroff.  No-regret algorithms for time-varying bayesian optimization.  In 2021 55th Annual Conference on Information Sciences and Systems (CISS), pages 1–6. IEEE, 2021.
  * [24] Yuntian Deng, Xingyu Zhou, Baekjin Kim, Ambuj Tewari, Abhishek Gupta, and Ness Shroff.  Weighted gaussian process bandits for non-stationary environments.  In International Conference on Artificial Intelligence and Statistics, pages 6909–6932. PMLR, 2022.
  * [25] Christopher KI Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * [26] Leonid V Kantorovich.  Mathematical methods of organizing and planning production.  Management science, 6(4):366–422, 1960.
  * [27] Anton Mallasto and Aasa Feragen.  Learning from uncertain curves: The 2-wasserstein metric for gaussian processes.  Advances in Neural Information Processing Systems, 30, 2017.
  * [28] Emilio Porcu, Moreno Bevilacqua, Robert Schaback, and Chris J Oates.  The mat$\backslash$’ern model: A journey through statistics, numerical analysis and machine learning.  arXiv preprint arXiv:2303.02759, 2023.
  * [29] Izrail Solomonovich Gradshteyn and Iosif Moiseevich Ryzhik.  Table of integrals, series, and products.  Academic press, 2014.
  * [30] Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton, Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy.  BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization.  In Advances in Neural Information Processing Systems 33, 2020.
  * [31] Wenzel Jakob, Jason Rhinelander, and Dean Moldovan.  pybind11 – seamless operability between c++11 and python, 2017.  https://github.com/pybind/pybind11.
  * [32] JHB Kemperman.  On the shannon capacity of an arbitrary channel.  In Indagationes Mathematicae (Proceedings), volume 77, pages 101–115. North-Holland, 1974.

##  Appendix A Wasserstein Distance at a Point in $\mathcal{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.

$\bm{\Delta}_{\mathcal{D}}^{-1}=\begin{pmatrix}\bm{E}&\bm{G}\\\
\bm{H}&\bm{F}\end{pmatrix}$ (13)

with
$\bm{\Delta}_{\mathcal{D}}=\bm{k}(\mathcal{D},\mathcal{D})+\sigma^{2}\bm{I}$
and $\bm{E},\bm{F},\bm{G},\bm{H}$ defined as

$\displaystyle\bm{E}$
$\displaystyle=\left(\lambda+\sigma^{2}-\bm{k}^{\top}\left((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}}\right)\bm{\Delta}_{\tilde{\mathcal{D}}}^{-1}\bm{k}\left((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}}\right)\right)^{-1},$
(14) $\displaystyle\bm{F}$
$\displaystyle=\left(\bm{\Delta}_{\tilde{\mathcal{D}}}-\frac{1}{\lambda+\sigma^{2}}\bm{k}\left((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}}\right)\bm{k}^{\top}\left((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}}\right)\right)^{-1},$
(15) $\displaystyle\bm{G}$
$\displaystyle=-\bm{E}\bm{k}^{\top}\left((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}}\right)\bm{\Delta}_{\tilde{\mathcal{D}}}^{-1},$
(16) $\displaystyle\bm{H}$
$\displaystyle=-\frac{1}{\lambda+\sigma^{2}}\bm{F}\bm{k}\left((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}}\right),$
(17)

where
$\bm{k}\left(\mathcal{X},\mathcal{Y}\right)=\left(k\left(\left(\bm{x}_{i},t_{i}\right),\left(\bm{x}_{j},t_{j}\right)\right)\right)_{\begin{subarray}{c}(\bm{x}_{i},t_{i})\in\mathcal{X}\\\
(\bm{x}_{j},t_{j})\in\mathcal{Y}\end{subarray}}$ and
$\bm{\Delta}_{\tilde{\mathcal{D}}}=\bm{k}(\tilde{\mathcal{D}},\tilde{\mathcal{D}})+\sigma^{2}\bm{I}$.

###### Proof.

The proof is trivial using the inverse of a block matrix:

$\displaystyle\begin{pmatrix}\bm{A}&\bm{B}\\\ \bm{C}&\bm{D}\end{pmatrix}^{-1}$
$\displaystyle=\begin{pmatrix}\left(\bm{A}-\bm{B}\bm{D}^{-1}\bm{C}\right)^{-1}&\bm{0}\\\
\bm{0}&\left(\bm{D}-\bm{C}\bm{A}^{-1}\bm{B}\right)^{-1}\end{pmatrix}\begin{pmatrix}\bm{I}&-\bm{B}\bm{D}^{-1}\\\
-\bm{C}\bm{A}^{-1}&\bm{I}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\left(\bm{A}-\bm{B}\bm{D}^{-1}\bm{C}\right)^{-1}&-\left(\bm{A}-\bm{B}\bm{D}^{-1}\bm{C}\right)^{-1}\bm{B}\bm{D}^{-1}\\\
-\left(\bm{D}-\bm{C}\bm{A}^{-1}\bm{B}\right)^{-1}\bm{C}\bm{A}^{-1}&\left(\bm{D}-\bm{C}\bm{A}^{-1}\bm{B}\right)^{-1}\end{pmatrix}.$
(18)

Note that $\bm{A}$ and $\bm{D}$ must be invertible.

We can use (18) to write $\bm{\Delta}_{\mathcal{D}}^{-1}$ as a function of
$\bm{\Delta}_{\tilde{\mathcal{D}}}^{-1}$, since

$\displaystyle\bm{\Delta}_{\mathcal{D}}^{-1}$
$\displaystyle=\begin{pmatrix}k((\bm{x}_{1},t_{1}),(\bm{x}_{1},t_{1}))+\sigma^{2}&\bm{k}^{\top}((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}})\\\
\bm{k}((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}})&\bm{\Delta}_{\tilde{\mathcal{D}}}\end{pmatrix}^{-1}$
$\displaystyle=\begin{pmatrix}\lambda+\sigma^{2}&\bm{k}^{\top}((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}})\\\
\bm{k}((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}})&\bm{\Delta}_{\tilde{\mathcal{D}}}\end{pmatrix}^{-1}.$
(19)

Note that, in (19), $\lambda+\sigma^{2}$ and
$\bm{\Delta}_{\tilde{\mathcal{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.

$\mu_{\mathcal{D}}(\bm{x},t)-\mu_{\tilde{\mathcal{D}}}(\bm{x},t)=ak((\bm{x},t),(\bm{x}_{1},t_{1}))+\bm{b}\bm{k}\left((\bm{x},t),\tilde{\mathcal{D}}\right),$
(20)

with $a=\bm{E}y_{1}+\bm{G}\tilde{\bm{y}}$,
$\bm{b}^{\top}=\bm{H}^{\top}y_{1}+\tilde{\bm{y}}^{\top}\left(\bm{F}-\bm{\Delta}_{\tilde{\mathcal{D}}}^{-1}\right)$
and $\tilde{\bm{y}}=(y_{2},\cdots,y_{n})$.

######  Proposition A.3.

$\begin{split}\left(\sigma_{\mathcal{D}}(\bm{x},t)-\sigma_{\tilde{\mathcal{D}}}(\bm{x},t)\right)^{2}\leq\bm{E}k((\bm{x},t),(\bm{x}_{1},t_{1}))^{2}&+k((\bm{x},t),(\bm{x}_{1},t_{1}))\bm{c}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})\\\
&+\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\bm{M}\bm{k}((\bm{x},t),\tilde{\mathcal{D}}),\end{split}$
(21)

with $\bm{c}=\bm{G}+\bm{H}^{\top}$ and
$\bm{M}=\bm{F}-\bm{\Delta}^{-1}_{\tilde{\mathcal{D}}}$.

We now prove Proposition A.2.

###### Proof.

According to (3)

$\mu_{\tilde{\mathcal{D}}}(\bm{x},t)=\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\bm{\Delta}_{\tilde{\mathcal{D}}}^{-1}\tilde{\bm{y}}.$
(22)

Applying the same definition for $\mu_{\mathcal{D}}(\bm{x},t)$ we have

$\displaystyle\mu_{\mathcal{D}}(\bm{x},t)$
$\displaystyle=\left(k((\bm{x},t),(\bm{x}_{1},t_{1})),k^{\top}((\bm{x},t),\tilde{\mathcal{D}})\right)\bm{\Delta}_{\mathcal{D}}^{-1}\bm{y}$
$\displaystyle=\left(k((\bm{x},t),(\bm{x}_{1},t_{1})),k^{\top}((\bm{x},t),\tilde{\mathcal{D}})\right)\begin{pmatrix}\bm{E}&\bm{G}\\\
\bm{H}&\bm{F}\end{pmatrix}\bm{y}$ (23)
$\displaystyle=\left(k((\bm{x},t),(\bm{x}_{1},t_{1})),k^{\top}((\bm{x},t),\tilde{\mathcal{D}})\right)\begin{pmatrix}\bm{E}y_{1}+\bm{G}\tilde{\bm{y}}\\\
\bm{H}y_{1}+\bm{F}\tilde{\bm{y}}\end{pmatrix}$
$\displaystyle=k((\bm{x},t),(\bm{x}_{1},t_{1}))\left(\bm{E}y_{1}+\bm{G}\tilde{\bm{y}}\right)+k^{\top}((\bm{x},t),\tilde{\mathcal{D}})\left(\bm{H}y_{1}+\bm{F}\tilde{\bm{y}}\right),$
(24)

where (23) follows from Lemma A.1.

Finally, we have

$\begin{split}\mu_{\mathcal{D}}(\bm{x},t)-\mu_{\tilde{\mathcal{D}}}(\bm{x},t)=k((\bm{x},t),(\bm{x}_{1},t_{1}))\left(\bm{E}y_{1}+\bm{G}\tilde{\bm{y}}\right)&+\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\left(\bm{H}y_{1}+\bm{F}\tilde{\bm{y}}\right)\\\
&-\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\bm{\Delta}_{\tilde{\mathcal{D}}}^{-1}\tilde{\bm{y}}\end{split}$

which can be reduced to

$\mu_{\mathcal{D}}(\bm{x},t)-\mu_{\tilde{\mathcal{D}}}(\bm{x},t)=ak((\bm{x},t),(\bm{x}_{1},t_{1}))+\bm{b}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})$
(25)

with $a=\bm{E}y_{1}+\bm{G}\tilde{\bm{y}}$ and
$\bm{b}^{\top}=\bm{H}^{\top}y_{1}+\tilde{\bm{y}}^{\top}\left(\bm{F}-\bm{\Delta}_{\tilde{\mathcal{D}}}^{-1}\right)$.
This concludes the proof. ∎

Finally, we prove Proposition A.3.

###### Proof.

As
$(\sigma_{\mathcal{D}}(\bm{x},t)-\sigma_{\tilde{\mathcal{D}}}(\bm{x},t))^{2}$
is hard to integrate, to get (8), we upper bound it by

$\displaystyle(\sigma_{\mathcal{D}}(\bm{x},t)-\sigma_{\tilde{\mathcal{D}}}(\bm{x},t))^{2}$
$\displaystyle=\sigma_{\mathcal{D}}^{2}(\bm{x},t)+\sigma^{2}_{\tilde{\mathcal{D}}}(\bm{x},t)-2\sigma_{\mathcal{D}}(\bm{x},t)\sigma_{\tilde{\mathcal{D}}}(\bm{x},t)$
$\displaystyle\leq\sigma_{\mathcal{D}}^{2}(\bm{x},t)+\sigma^{2}_{\tilde{\mathcal{D}}}(\bm{x},t)-2\sigma^{2}_{\mathcal{D}}(\bm{x},t)$
(26)
$\displaystyle=\sigma^{2}_{\tilde{\mathcal{D}}}(\bm{x},t)-\sigma^{2}_{\mathcal{D}}(\bm{x},t)$
(27)

where (26) follows from
$\sigma_{\mathcal{D}}(\bm{x},t)\leq\sigma_{\tilde{\mathcal{D}}}(\bm{x},t)$.

Now, (4) yields that for $\bm{X}=\tilde{\mathcal{D}}$,

$\sigma^{2}_{\tilde{\mathcal{D}}}(\bm{x},t)=\lambda-\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\bm{\Delta}_{\tilde{\mathcal{D}}}^{-1}\bm{k}((\bm{x},t),\tilde{\mathcal{D}}).$
(28)

and for $\bm{X}=\mathcal{D}=\tilde{\mathcal{D}}\cup\\{(\bm{x}_{1},t_{1})\\}$

$\displaystyle\sigma^{2}_{\mathcal{D}}(\bm{x},t)$
$\displaystyle=\lambda-\left(k((\bm{x},t),(\bm{x}_{1},t_{1})),\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\right)\bm{\Delta}_{\mathcal{D}}^{-1}\begin{pmatrix}k((\bm{x},t),(\bm{x}_{1},t_{1}))\\\
\bm{k}((\bm{x},t),\tilde{\mathcal{D}})\end{pmatrix}$
$\displaystyle=\lambda-\left(k((\bm{x},t),(\bm{x}_{1},t_{1})),\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\right)\begin{pmatrix}\bm{E}&\bm{G}\\\
\bm{H}&\bm{F}\end{pmatrix}\begin{pmatrix}k((\bm{x},t),(\bm{x}_{1},t_{1}))\\\
\bm{k}((\bm{x},t),\tilde{\mathcal{D}})\end{pmatrix}$ (29)
$\displaystyle=\lambda-\left(k((\bm{x},t),(\bm{x}_{1},t_{1})),\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\right)\begin{pmatrix}\bm{E}k((\bm{x},t),(\bm{x}_{1},t_{1}))+\bm{G}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})\\\
\bm{H}k((\bm{x},t),(\bm{x}_{1},t_{1}))+\bm{F}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})\end{pmatrix}$
(30)

where (29) follows from Lemma A.1. Developing the dot product in (30), we have

$\begin{split}\sigma^{2}_{\mathcal{D}}(\bm{x},t)=\lambda&-\bm{E}k^{2}((\bm{x},t),(\bm{x}_{1},t_{1}))-k((\bm{x},t),(\bm{x}_{1},t_{1}))\bm{G}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})\\\
&-k((\bm{x},t),(\bm{x}_{1},t_{1}))\bm{H}^{\top}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})-\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\bm{F}\bm{k}((\bm{x},t),\tilde{\mathcal{D}}).\end{split}$
(31)

Combining (28) and (31), we get

$\begin{split}\sigma^{2}_{\tilde{\mathcal{D}}}(\bm{x},t)-\sigma^{2}_{\mathcal{D}}(\bm{x},t)=\bm{E}k((\bm{x},t),(\bm{x}_{1},t_{1}))^{2}&+k((\bm{x},t),(\bm{x}_{1},t_{1}))\bm{c}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})\\\
&+\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\bm{M}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})\end{split}$
(32)

where $\bm{c}=\bm{G}+\bm{H}^{\top}$ and
$\bm{M}=\bm{F}-\bm{\Delta}^{-1}_{\tilde{\mathcal{D}}}$.

Combining (27) and (32) concludes the proof. ∎

##  Appendix B Wasserstein Distance on $\mathcal{F}_{t_{0}}$

In this appendix, we extend the computation of the Wasserstein distance
between two posterior Gaussian distributions at an arbitrary point in
$\mathcal{F}_{t_{0}}$ (see Propositions A.2 and A.3) to the Wasserstein
distance between $\mathcal{GP}_{\mathcal{D}}$ and
$\mathcal{GP}_{\tilde{\mathcal{D}}}$ on $\mathcal{F}_{t_{0}}$. We also provide
an upper bound for the Wasserstein distance between a posterior GP conditioned
on $\mathcal{D}$ (i.e., $\mathcal{GP}_{\mathcal{D}}$) and the prior GP (i.e.,
$\mathcal{GP}_{\emptyset}$).

Let us start by proving the following lemma.

######  Lemma B.1.

Let $t_{0}$ be the present time and
$\mathcal{D}=\left\\{((\bm{x}_{i},t_{i}),y_{i})\right\\}_{i\in\llbracket
1,n\rrbracket}$ be a dataset of observations before $t_{0}$. Let
$\tilde{\mathcal{D}}=\left\\{((\bm{x}_{i},t_{i}),y_{i})\right\\}_{i\in\llbracket
2,n\rrbracket}$ and $\mathcal{F}_{t_{0}}=\mathcal{S}\times[t_{0},+\infty)$ be
the domain of future predictions. Then, for any pair of observations
$(\bm{x}_{i},t_{i},\cdot),(\bm{x}_{j},t_{j},\cdot)$ from $\mathcal{D}$ we have

$\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}k\left((\bm{x},t),(\bm{x}_{i},t_{i})\right)k\left((\bm{x},t),(\bm{x}_{j},t_{j})\right)d\bm{x}dt\leq\lambda^{2}(k_{S}*k_{S})(\bm{x}_{j}-\bm{x}_{i})(k_{T}*k_{T})_{t_{0}-t_{i}}^{+\infty}(t_{j}-t_{i})$
(33)

with $(f*g)$ the convolution between $f$ and $g$ and $(f*g)_{a}^{b}$ the
convolution between $f$ and $g$ restricted to the interval $[a,b]$.

###### 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
$\mathcal{F}_{t_{0}}=\mathcal{S}\times[t_{0},+\infty)$, we have

$I=\lambda^{2}\oint_{\mathcal{\mathcal{S}}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x}\int_{t_{0}}^{+\infty}k_{T}(|t-t_{i}|)k_{T}(|t-t_{j}|)dt.$
(34)

Let us focus on the integral on the spatial domain $\mathcal{S}$. Depending on
the expression of the covariance function, this integral can be quite
difficult to compute. Let us turn this expression into a more pleasant upper
bound

$\displaystyle\oint_{\mathcal{\mathcal{S}}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x}$
$\displaystyle\leq\oint_{\mathbb{R}^{d}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x}$
(35)
$\displaystyle=\oint_{\mathbb{R}^{d}}k_{S}(||\bm{u}||_{2})k_{S}(||\bm{x}_{j}-\bm{x}_{i}-\bm{u}||_{2})d\bm{u}$
(36)

where (35) holds since $k_{S}$ is positive and (36) comes from the change of
variable $\bm{u}=\bm{x}-\bm{x}_{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})(\bm{x}_{j}-\bm{x}_{i})$. A similar change of variable can be
made regarding the time integral, with $v=t-t_{i}$:

$\int_{t_{0}}^{+\infty}k_{T}(|t-t_{i}|)k_{T}(|t-t_{j}|)dt=\int_{t_{0}-t_{i}}^{+\infty}k_{T}(|v|)k_{T}(|t_{j}-t_{i}-v|)dv.$
(37)

Note that (37) is also a convolution, but restricted to the interval
$[t_{0}-t_{i},+\infty)$. We denote this convolution on a restricted interval
$(k_{T}*k_{T})_{t_{0}-t_{i}}^{+\infty}(t_{j}-t_{i})$.

Combining (34), (36) and (37) yields Lemma B.1. ∎

We now prove the following lemma.

######  Lemma B.2.

Let $t_{0}$ be the present time and
$\mathcal{D}=\left\\{((\bm{x}_{i},t_{i}),y_{i})\right\\}_{i\in\llbracket
1,n\rrbracket}$ be a dataset of observations before $t_{0}$. Let
$\tilde{\mathcal{D}}=\left\\{((\bm{x}_{i},t_{i}),y_{i})\right\\}_{i\in\llbracket
2,n\rrbracket}$ and $\mathcal{F}_{t_{0}}=\mathcal{S}\times[t_{0},+\infty)$ the
domain of future predictions. Then, we have

$\begin{split}\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}(\mu_{\mathcal{D}}(\bm{x},t)-\mu_{\tilde{\mathcal{D}}}(\bm{x},t))^{2}d\bm{x}dt\leq\lambda^{2}a^{2}C((\bm{x}_{1},t_{1}),(\bm{x}_{1},&t_{1}))+2\lambda^{2}a\bm{b}C((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}})\\\
&+\lambda^{2}\bm{1}^{\top}\left(\bm{b}\bm{b}^{\top}\odot
C(\tilde{\mathcal{D}},\tilde{\mathcal{D}})\right)\bm{1},\end{split}$ (38)

with $a=\bm{E}y_{1}+\bm{G}\tilde{\bm{y}}$ and
$\bm{b}^{\top}=\bm{H}^{\top}y_{1}+\tilde{\bm{y}}^{\top}\left(\bm{F}-\bm{\Delta}_{\tilde{\mathcal{D}}}^{-1}\right)$,
$C(\mathcal{X},\mathcal{Y})=\left((k_{S}*k_{S})(\bm{x}_{j}-\bm{x}_{i})(k_{T}*k_{T})_{t_{0}-t_{i}}^{+\infty}(t_{j}-t_{i})\right)_{\begin{subarray}{c}(\bm{x}_{i},t_{i})\in\mathcal{X}\\\
(\bm{x}_{j},t_{j})\in\mathcal{Y}\end{subarray}}$, $\bm{1}$ the conformable
vector of ones, $\odot$ the Hadamard product, $(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]$.

###### Proof.

By Proposition A.2,

$\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}(\mu_{\mathcal{D}}(\bm{x},t)-\mu_{\tilde{\mathcal{D}}}(\bm{x},t))^{2}d\bm{x}dt=\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\left(ak((\bm{x},t),(\bm{x}_{1},t_{1}))+\bm{b}\bm{k}\left((\bm{x},t),\tilde{\mathcal{D}}\right)\right)^{2}d\bm{x}dt$
(39)

with $a=\bm{E}y_{1}+\bm{G}\tilde{\bm{y}}$ and
$\bm{b}^{\top}=\bm{H}^{\top}y_{1}+\tilde{\bm{y}}^{\top}\left(\bm{F}-\bm{\Delta}_{\tilde{\mathcal{D}}}^{-1}\right)$.

Expanding the square, we get three different integrals, denoted $A,B$ and $C$
with

$\displaystyle A$
$\displaystyle=\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}a^{2}k^{2}((\bm{x},t),(\bm{x}_{1},t_{1}))d\bm{x}dt,$
$\displaystyle B$
$\displaystyle=\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}2a\bm{b}k((\bm{x},t),(\bm{x}_{1},t_{1}))\bm{k}((\bm{x},t),\tilde{\mathcal{D}})d\bm{x}dt,$
$\displaystyle C$
$\displaystyle=\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\bm{b}\bm{k}\left((\bm{x},t),\tilde{\mathcal{D}}\right)\bm{k}^{\top}\left((\bm{x},t),\tilde{\mathcal{D}}\right)\bm{b}^{\top}d\bm{x}dt.$

Using Lemma B.1, computing an upper bound for $A$, $B$ and $C$ is immediate.
In fact,

$\displaystyle A$
$\displaystyle\leq\lambda^{2}a^{2}C((\bm{x}_{1},t_{1}),(\bm{x}_{1},t_{1}))$
(40) $\displaystyle B$ $\displaystyle\leq
2\lambda^{2}a\bm{b}C((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}})$ (41)
$\displaystyle C$
$\displaystyle\leq\lambda^{2}\bm{1}^{\top}\left(\bm{b}\bm{b}^{\top}\odot
C(\tilde{\mathcal{D}},\tilde{\mathcal{D}})\right)\bm{1}$ (42)

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
$\mathcal{D}=\left\\{((\bm{x}_{i},t_{i}),y_{i})\right\\}_{i\in\llbracket
1,n\rrbracket}$ be a dataset of observations before $t_{0}$. Let
$\tilde{\mathcal{D}}=\left\\{((\bm{x}_{i},t_{i}),y_{i})\right\\}_{i\in\llbracket
2,n\rrbracket}$ and $\mathcal{F}_{t_{0}}=\mathcal{S}\times[t_{0},+\infty)$ the
domain of future predictions. Then, we have

$\begin{split}\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}(\sigma_{\mathcal{D}}(\bm{x},t)-\sigma_{\tilde{\mathcal{D}}}(\bm{x},t))^{2}d\bm{x}dt\leq\lambda^{2}\bm{E}C((\bm{x}_{1},t_{1}),(\bm{x}_{1},&t_{1}))+\lambda^{2}\bm{c}C((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}})\\\
&+\lambda^{2}\bm{1}^{\top}\left(\bm{M}\odot
C(\tilde{\mathcal{D}},\tilde{\mathcal{D}})\right)\bm{1},\end{split}$ (43)

with $\bm{c}=\bm{G}+\bm{H}^{\top}$,
$\bm{M}=\bm{F}-\bm{\Delta}^{-1}_{\tilde{\mathcal{D}}}$,
$C(\mathcal{X},\mathcal{Y})=\left((k_{S}*k_{S})(\bm{x}_{j}-\bm{x}_{i})(k_{T}*k_{T})_{t_{0}-t_{i}}^{+\infty}(t_{j}-t_{i})\right)_{\begin{subarray}{c}(\bm{x}_{i},t_{i})\in\mathcal{X}\\\
(\bm{x}_{j},t_{j})\in\mathcal{Y}\end{subarray}}$, $\bm{1}$ the conformable
vector of ones, $\odot$ the Hadamard product, $(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]$.

###### Proof.

This proof is conceptually identical to the previous one. By Proposition A.3,

$\begin{split}\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}(\sigma_{\mathcal{D}}(\bm{x},t)-\sigma_{\tilde{\mathcal{D}}}(\bm{x},t))^{2}d\bm{x}dt\leq\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}(\bm{E}k((\bm{x},t),&(\bm{x}_{1},t_{1}))^{2}+k((\bm{x},t),(\bm{x}_{1},t_{1}))\bm{c}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})\\\
&+\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\bm{M}\bm{k}((\bm{x},t),\tilde{\mathcal{D}}))d\bm{x}dt\end{split}$
(44)

where $\bm{c}=\bm{G}+\bm{H}^{\top}$ and
$\bm{M}=\bm{F}-\bm{\Delta}^{-1}_{\tilde{\mathcal{D}}}$.

By linearity of the integral, the RHS of (44) can be split into three
different integrals $A$, $B$ and $C$ where

$\displaystyle A$
$\displaystyle=\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\bm{E}k((\bm{x},t),(\bm{x}_{1},t_{1}))^{2}d\bm{x}dt$
$\displaystyle B$
$\displaystyle=\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}k((\bm{x},t),(\bm{x}_{1},t_{1}))\bm{c}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})d\bm{x}dt$
$\displaystyle C$
$\displaystyle=\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\bm{k}^{\top}((\bm{x},t),\tilde{\mathcal{D}})\bm{M}\bm{k}((\bm{x},t),\tilde{\mathcal{D}})d\bm{x}dt$

Once again, Lemma B.1 can be used to compute an upper bound for $A$, $B$ and
$C$. In fact,

$\displaystyle A$
$\displaystyle\leq\lambda^{2}\bm{E}C((\bm{x}_{1},t_{1}),(\bm{x}_{1},t_{1}))$
(45) $\displaystyle B$
$\displaystyle\leq\lambda^{2}\bm{c}C((\bm{x}_{1},t_{1}),\tilde{\mathcal{D}})$
(46) $\displaystyle C$
$\displaystyle\leq\lambda^{2}\bm{1}^{\top}\left(\bm{M}\odot
C(\tilde{\mathcal{D}},\tilde{\mathcal{D}})\right)\bm{1}$ (47)

Adding (45), (46) and (47) together concludes the proof. ∎

Together, Lemmas B.1, B.2 and B.3 yield the first part of Theorem 4.1. For the
second part of Theorem 4.1, we prove the following lemma.

######  Lemma B.4.

Let $t_{0}$ be the present time and
$\mathcal{D}=\left\\{((\bm{x}_{i},t_{i}),y_{i})\right\\}_{i\in\llbracket
1,n\rrbracket}$ be a dataset of observations before $t_{0}$. Let
$\mathcal{F}_{t_{0}}=\mathcal{S}\times[t_{0},+\infty)$ the domain of future
predictions. Then, we have

$\begin{split}W^{2}_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset})\leq\lambda^{2}\left(\bm{a}^{\top}C\left(\mathcal{D},\mathcal{D}\right)\bm{a}+\bm{1}^{\top}\left(\bm{\Delta}^{-1}\odot
C\left(\mathcal{D},\mathcal{D}\right)\right)\bm{1}\right)\end{split}$ (48)

with $\bm{a}=\bm{\Delta}^{-1}\bm{y}$ and
$C(\mathcal{X},\mathcal{Y})=\left((k_{S}*k_{S})(\bm{x}_{j}-\bm{x}_{i})(k_{T}*k_{T})_{t_{0}-t_{i}}^{+\infty}(t_{j}-t_{i})\right)_{\begin{subarray}{c}(\bm{x}_{i},t_{i})\in\mathcal{X}\\\
(\bm{x}_{j},t_{j})\in\mathcal{Y}\end{subarray}}$, $\bm{1}$ the conformable
vector of ones, $\odot$ the Hadamard product, $(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]$.

###### Proof.

Recall that, according to $\mathcal{GP}_{\emptyset}$,
$f(\bm{x},t)\sim\mathcal{N}\left(0,\lambda\right)$ for any point
$(\bm{x},t)\in\mathcal{F}_{t_{0}}$. Consequently,

$\displaystyle W_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset})$
$\displaystyle=\left(\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\mu^{2}_{\mathcal{D}}(\bm{x},t)d\bm{x}dt+\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\left(\sqrt{\lambda}-\sigma_{\mathcal{D}}(\bm{x},t)\right)^{2}d\bm{x}dt\right)^{\frac{1}{2}}.$
(49)

These two integrals in (49) can be computed with the same techniques as above.
For the mean integral, we have

$\displaystyle\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\mu^{2}_{\mathcal{D}}(\bm{x},t)d\bm{x}dt$
$\displaystyle=\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\bm{y}^{\top}\bm{\Delta}^{-\top}\bm{k}^{\top}((\bm{x},t),\mathcal{D})\bm{k}((\bm{x},t),\mathcal{D})\bm{\Delta}^{-1}\bm{y}d\bm{x}dt$
$\displaystyle=\bm{y}^{\top}\bm{\Delta}^{-\top}\left(\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\bm{k}^{\top}((\bm{x},t),\mathcal{D})\bm{k}((\bm{x},t),\mathcal{D})d\bm{x}dt\right)\bm{\Delta}^{-1}\bm{y}$
$\displaystyle\leq\lambda^{2}\bm{y}^{\top}\bm{\Delta}^{-\top}C(\mathcal{D},\mathcal{D})\bm{\Delta}^{-1}\bm{y}.$
(50)

Regarding the variance integral in (49), we have

$\displaystyle\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\left(\sqrt{\lambda}-\sigma_{\mathcal{D}}(\bm{x},t)\right)^{2}d\bm{x}dt$
$\displaystyle=\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\left(\lambda-2\sqrt{\lambda}\sigma_{\mathcal{D}}(\bm{x},t)+\sigma^{2}_{\mathcal{D}}(\bm{x},t)\right)d\bm{x}dt$
$\displaystyle\leq\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\left(\lambda-\sigma^{2}_{\mathcal{D}}(\bm{x},t)\right)d\bm{x}dt$
(51)
$\displaystyle=\oint_{\mathcal{S}}\int_{t_{0}}^{\infty}\bm{k}^{\top}((\bm{x},t),\mathcal{D})\bm{\Delta}^{-1}\bm{k}((\bm{x},t),\mathcal{D})d\bm{x}dt$
(52) $\displaystyle\leq\lambda^{2}\bm{1}^{\top}\left(\bm{\Delta}^{-1}\odot
C\left(\mathcal{D},\mathcal{D}\right)\right)\bm{1}$ (53)

where (51) holds because
$\sqrt{\lambda}\sigma_{\mathcal{D}}(\bm{x},t)\geq\sigma^{2}_{\mathcal{D}}(\bm{x},t)$
and (52) holds because of (4).

Together, (50) and (53) conclude the proof. ∎

By combining Lemmas B.2 and B.3, the proof of Theorem 4.1 is immediate.

##  Appendix C Approximation Error

In Appendix B, we provided a computationally tractable upper bound of
$W_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\tilde{\mathcal{D}}})$ and
$W_{2}(\mathcal{GP}_{\mathcal{D}},\mathcal{GP}_{\emptyset})$. 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 $\mathcal{S}=[0,1]^{d}$. In
Appendix B, we have to approximate the Wasserstein distance because the
integration over $\mathcal{S}$ in (34) is difficult to compute. The upper
bound proposed in (35) is

$\oint_{[0,1]^{d}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x}\leq\oint_{\mathbb{R}^{d}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x}.$
(54)

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
$\mathcal{S}=[0,1]^{d}$ to $\mathbb{R}^{d}$ has a bounded impact on the
result. Clearly, a first absolute approximation error for (54) is

$\oint_{\left(\mathbb{R}\setminus[0,1]\right)^{d}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x}.$

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:

$\begin{split}A(\bm{x}_{i},\bm{x}_{j};l_{S})=\min&\left\\{\oint_{\left(\mathbb{R}\setminus[0,1]\right)^{d}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x},\right.\\\
&\left.1-\oint_{[0,1]^{d}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x}\right\\}.\end{split}$
(55)

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 $\left(\bm{x}_{i},\bm{x}_{j}\right)\in\mathcal{S}^{2}$, with
$\mathcal{S}=[0,1]^{d}$. Let $k_{S}$ be a SE kernel with lengthscale $l_{S}$.
Then,

$\operatorname*{arg\,max}_{l_{S}\in\mathbb{R}^{+}}A(\bm{x}_{i},\bm{x}_{j};l_{S})=\frac{1}{\sqrt{\pi}}e^{\frac{1}{2}W_{0}\left(\frac{\pi||\bm{x}_{i}-\bm{x}_{j}||_{2}^{2}}{2d}\right)},$
(56)

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

$\displaystyle\oint_{\left(\mathbb{R}\setminus[0,1]\right)^{d}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x}$
$\displaystyle=1-\oint_{[0,1]^{d}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x},$
$\displaystyle\oint_{\mathbb{R}^{d}}k_{S}(||\bm{x}-\bm{x}_{i}||_{2})k_{S}(||\bm{x}-\bm{x}_{j}||_{2})d\bm{x}$
$\displaystyle=1,$
$\displaystyle\pi^{\frac{d}{2}}l_{S}^{d}e^{\frac{-||\bm{x}_{i}-\bm{x}_{j}||_{2}^{2}}{4l_{S}^{2}}}$
$\displaystyle=1,$ (57)

where (57) uses a result derived in Appendix D and reported in Table 3.
Solving for $l_{S}$ concludes the proof. ∎

Figure 4: (Top row) Absolute approximation error (55) with respect to the
spatial lengthscale $l_{S}$ for a $1$, $3$ and $5$-dimensional spatial domain.
Both error terms in (55) are shown in orange and green dashed lines,
respectively. Finally, the critical lengthscale (56) is shown as a red
vertical line. In this example, $k_{S}$ is a SE correlation function. (Bottom
row) Relative approximation error with respect to the spatial lengthscale
$l_{S}$. The color codes are the same.

Let 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 $\mathcal{S}$ is computed with a Monte-Carlo numerical
integration technique. The results are shown in Figure 4. 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 4 shows that the same observations
apply to the relative errors.

Figure 4 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 4, which clearly
increases in magnitude with the dimensionality of the spatial domain.

Figure 5: Relative error between the criterion (9) and its approximation (12),
with respect to the spatial lengthscale $l_{S}$ for a $1$, $3$ and
$5$-dimensional spatial domain. The relative error computed with both terms in
(55) are shown as orange and green dashed lines, respectively. In this
example, $k_{S}$ is a SE correlation function.

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

##  Appendix D Convolutions of Usual Covariance Functions

Table 3: Analytic forms for the convolution of usual spatial covariance functions. $\Gamma$ is the Gamma function, $J_{\alpha}$ is a Bessel function of the first kind of order $\alpha$, $K_{\alpha}$ is a modified Bessel function of the second kind of order $\alpha$. Covariance Function $k_{S}$ | $(k_{S}*k_{S})(\bm{x})$
---|---
Squared-Exponential ($l_{S}$) | $\pi^{\frac{d}{2}}l_{S}^{d}e^{\frac{-||\bm{x}||^{2}_{2}}{4l_{S}^{2}}}$
Matérn ($\nu,l_{S}$) | $\frac{2^{\frac{d}{2}-2\nu+1}\pi^{\frac{d}{2}}\Gamma(\nu+\frac{d}{2})^{2}}{\Gamma(\nu)^{2}\Gamma(2\nu+d)}\left(\frac{\sqrt{2\nu}}{l_{S}}\right)^{2\nu-\frac{d}{2}}||\bm{x}||_{2}^{2\nu+\frac{d}{2}}K_{2\nu+\frac{d}{2}}\left(\frac{||\bm{x}||_{2}\sqrt{2\nu}}{l_{S}}\right)$
Table 4: Analytic forms for the convolution of usual temporal covariance functions on the interval $[t_{0}-t_{i},+\infty)$. Note that erf is the error function. Also, for the sake of brevity, the terms $C_{k_{1}k_{2}}$, $P_{k_{1}k_{2}}$ and $Q_{k_{1}k_{2}}$ are defined in Appendix D.2.2. Covariance Function $k_{T}$ | $(k_{T}*k_{T})_{t_{0}-t_{i}}^{+\infty}(t_{j}-t_{i})$
---|---
Squared-Exponential ($l_{T}$) | $\frac{\sqrt{\pi}l_{T}}{2}e^{\frac{-(t_{i}-t_{j})^{2}}{2l_{T}^{2}}}\left(1-\text{erf}\left(\frac{2t_{0}-t_{i}-t_{j}}{2l_{T}}\right)\right)$
Matérn ($\nu=p+\frac{1}{2},l_{T}$) | $\sum_{k_{1}=0}^{p}\sum_{k_{2}=0}^{p}C_{k_{1}k_{2}}e^{\frac{-\sqrt{2p+1}(2t_{0}-t_{i}-t_{j})}{l_{T}}}P_{k_{1}k_{2}}(t_{0},t_{i},t_{j})$

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

###  D.1 Spatial Covariance Functions

In this subsection, we compute specifically the analytic forms for the
convolution of usual spatial covariance functions. We rely on a direct
consequence of the convolution theorem, that is
$(k*k)(\bm{x})=\mathcal{F}^{-1}\left(\mathcal{F}^{2}\left(k\right)\right)(\bm{x})$,
with $\mathcal{F}(f)$ denoting the Fourier transform of $f$ and
$\mathcal{F}^{-1}(f)$ the inverse Fourier transform of $f$.

The Fourier transform $\mathcal{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 [25] for more details).
Furthermore, it is known that if $k$ is isotropic (i.e. it can be written as a
function of $r=||\bm{x}||_{2}$), then its spectral density $S(\bm{s})$ can be
written as a function of $s=||\bm{s}||_{2}$. In that case, the two functions
are linked by the pair of transforms (see [25])

$\displaystyle k(r)$
$\displaystyle=\frac{2\pi}{r^{\frac{d}{2}-1}}\int_{0}^{\infty}S(s)J_{\frac{d}{2}-1}(2\pi
rs)s^{\frac{d}{2}}ds$ (58) $\displaystyle S(s)$
$\displaystyle=\frac{2\pi}{s^{\frac{d}{2}-1}}\int_{0}^{\infty}k(r)J_{\frac{d}{2}-1}(2\pi
rs)r^{\frac{d}{2}}dr$ (59)

where $J_{\frac{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=||\bm{x}||_{2}$ and $s=||\bm{s}||_{2}$. Then,

$(k*k)(r)=\frac{2\pi}{r^{\frac{d}{2}-1}}\int_{0}^{\infty}S^{2}(s)J_{\frac{d}{2}-1}(2\pi
rs)s^{\frac{d}{2}}ds$ (60)

where $J_{\frac{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}$.

####  D.1.1 Squared-Exponential Covariance Function

######  Lemma D.2.

Let $k_{S}$ be a Squared-Exponential covariance function (see Table 1), with
lengthscale $l_{S}>0$. Then,

$(k_{S}*k_{S})(\bm{x})=\pi^{\frac{d}{2}}l_{S}^{d}e^{\frac{-||\bm{x}||^{2}_{2}}{4l_{S}^{2}}}.$
(61)

###### Proof.

The spectral density of a Squared-Exponential covariance function $k_{S}$ is
(see [25])

$S(s)=(2\pi l_{S}^{2})^{\frac{d}{2}}e^{-2\pi^{2}l_{S}^{2}s^{2}}.$ (62)

According to Corollary D.1, we have

$\displaystyle(k_{S}*k_{S})(\bm{x})$
$\displaystyle=\frac{2\pi}{r^{\frac{d}{2}-1}}\int_{0}^{\infty}S^{2}(s)J_{\frac{d}{2}-1}(2\pi
rs)s^{\frac{d}{2}}ds$
$\displaystyle=\frac{2\pi}{r^{\frac{d}{2}-1}}\int_{0}^{\infty}(2\pi
l_{S}^{2})^{d}e^{-4\pi^{2}l_{S}^{2}s^{2}}J_{\frac{d}{2}-1}(2\pi
rs)s^{\frac{d}{2}}ds$
$\displaystyle=\frac{\left(2\pi\right)^{d+1}l_{S}^{2d}}{r^{\frac{d}{2}-1}}\int_{0}^{\infty}e^{-4\pi^{2}l_{S}^{2}s^{2}}J_{\frac{d}{2}-1}(2\pi
rs)s^{\frac{d}{2}}ds$ (63)

where $r=||\bm{x}||_{2}$ and $J_{\frac{d}{2}-1}$ is a Bessel function of the
first kind of order $\frac{d}{2}-1$.

It is known (see [29]) that

$\int_{0}^{\infty}e^{-\alpha x^{2}}x^{\nu+1}J_{\nu}(\beta
x)dx=\frac{\beta^{\nu}}{(2\alpha)^{\nu+1}}e^{\frac{-\beta^{2}}{4\alpha}}.$

Therefore,

$\displaystyle(k_{S}*k_{S})(\bm{x})$
$\displaystyle=\frac{\left(2\pi\right)^{d+1}l_{S}^{2d}}{r^{\frac{d}{2}-1}}\frac{(2\pi
r)^{\frac{d}{2}-1}}{(8\pi^{2}l_{S}^{2})^{\frac{d}{2}}}e^{\frac{-4\pi^{2}r^{2}}{16\pi^{2}l_{S}^{2}}}$
$\displaystyle=\pi^{\frac{d}{2}}l_{S}^{d}e^{\frac{-r^{2}}{4l_{S}^{2}}}.$ (64)

Replacing $r$ by $||\bm{x}||_{2}$ in (64) concludes the proof. ∎

####  D.1.2 Matérn Covariance Function

######  Lemma D.3.

Let $k_{S}$ be a Matérn covariance function (see Table 1), with smoothness
parameter $\nu>0$ and lengthscale $l_{S}>0$. Then,

$(k_{S}*k_{S})(\bm{x})=\frac{2^{\frac{d}{2}-2\nu+1}\pi^{\frac{d}{2}}\Gamma(\nu+\frac{d}{2})^{2}}{\Gamma(\nu)^{2}\Gamma(2\nu+d)}\left(\frac{\sqrt{2\nu}}{l_{S}}\right)^{2\nu-\frac{d}{2}}||\bm{x}||_{2}^{2\nu+\frac{d}{2}}K_{2\nu+\frac{d}{2}}\left(\frac{||\bm{x}||_{2}\sqrt{2\nu}}{l_{S}}\right),$
(65)

where $\Gamma$ is the Gamma function and $K_{\alpha}$ is a modified Bessel
function of the second kind of order $\alpha$.

###### Proof.

The spectral density of a Matérn covariance function $k_{S}$ is (see [25])

$S(s)=\frac{2^{d}\pi^{\frac{d}{2}}\Gamma(\nu+\frac{d}{2})(2\nu)^{\nu}}{\Gamma(\nu)l_{S}^{2\nu}}\left(\frac{2\nu}{l_{S}^{2}}+4\pi^{2}s^{2}\right)^{-\nu-\frac{d}{2}}$
(66)

where $\Gamma$ is the Gamma function.

According to Corollary D.1, we have

$\displaystyle(k_{S}*k_{S})(\bm{x})$
$\displaystyle=\frac{2\pi}{r^{\frac{d}{2}-1}}\int_{0}^{\infty}S^{2}(s)J_{\frac{d}{2}-1}(2\pi
rs)s^{\frac{d}{2}}ds$
$\displaystyle=\frac{2\pi}{r^{\frac{d}{2}-1}}\int_{0}^{\infty}\frac{2^{2d}\pi^{d}\Gamma^{2}(\nu+\frac{d}{2})(2\nu)^{2\nu}}{\Gamma^{2}(\nu)l_{S}^{4\nu}}\left(\frac{2\nu}{l_{S}^{2}}+4\pi^{2}s^{2}\right)^{-2\nu-d}J_{\frac{d}{2}-1}(2\pi
rs)s^{\frac{d}{2}}ds$
$\displaystyle=\frac{2^{2d+1}\pi^{d+1}\Gamma^{2}(\nu+\frac{d}{2})(2\nu)^{2\nu}}{\Gamma^{2}(\nu)l_{S}^{4\nu}r^{\frac{d}{2}-1}}\int_{0}^{\infty}\left(\frac{2\nu}{l_{S}^{2}}+4\pi^{2}s^{2}\right)^{-2\nu-d}J_{\frac{d}{2}-1}(2\pi
rs)s^{\frac{d}{2}}ds$
$\displaystyle=\frac{2^{\frac{3d}{2}}\pi^{\frac{d}{2}}\Gamma^{2}(\nu+\frac{d}{2})(2\nu)^{2\nu}}{\Gamma^{2}(\nu)l_{S}^{4\nu}r^{\frac{d}{2}-1}}\int_{0}^{\infty}\left(\frac{2\nu}{l_{S}^{2}}+u^{2}\right)^{-2\nu-d}J_{\frac{d}{2}-1}(ru)u^{\frac{d}{2}}du$
(67)

where $J_{\frac{d}{2}-1}$ is a Bessel function of the first kind of order
$\frac{d}{2}-1$, and (67) comes from the change of variable $u=2\pi s$.

It is known (see [29]) that

$\int_{0}^{\infty}(a^{2}+x^{2})^{-(\mu+1)}J_{\alpha}(bx)x^{\alpha+1}dx=\frac{a^{\alpha-\mu}b^{\mu}}{2^{\mu}\Gamma(\mu+1)}K_{\mu-\alpha}(ab)$

where $K_{\mu-\alpha}$ is a modified Bessel function of the second kind of
order $\mu-\alpha$.

Therefore,

$\displaystyle(k_{S}*k_{S})(\bm{x})$
$\displaystyle=\frac{2^{\frac{3d}{2}}\pi^{\frac{d}{2}}\Gamma^{2}(\nu+\frac{d}{2})(2\nu)^{2\nu}}{\Gamma^{2}(\nu)l_{S}^{4\nu}r^{\frac{d}{2}-1}}\left(\frac{\sqrt{2\nu}}{l_{S}}\right)^{-2\nu-\frac{d}{2}}\frac{r^{2\nu+d-1}}{2^{2\nu+d-1}\Gamma(2\nu+d)}K_{2\nu+\frac{d}{2}}\left(\frac{r\sqrt{2\nu}}{l_{S}}\right)$
$\displaystyle=\frac{2^{\frac{d}{2}-2\nu+1}\pi^{\frac{d}{2}}\Gamma(\nu+\frac{d}{2})^{2}}{\Gamma(\nu)^{2}\Gamma(2\nu+d)}\left(\frac{\sqrt{2\nu}}{l_{S}}\right)^{2\nu-\frac{d}{2}}r^{2\nu+\frac{d}{2}}K_{2\nu+\frac{d}{2}}\left(\frac{r\sqrt{2\nu}}{l_{S}}\right).$
(68)

Replacing $r$ by $||\bm{x}||_{2}$ in (68) concludes the proof. ∎

###  D.2 Temporal Covariance Functions

In this section, we derive the analytic expression of the convolutions of the
most popular temporal covariance functions $k_{T}$ restricted to the interval
$[t_{0}-t_{i},+\infty)$. Therefore, we compute many integrals of the form

$\int_{t_{0}-t_{i}}^{+\infty}k_{T}(|t|)k_{T}(|t_{j}-t_{i}-t|)dt,$

which can be rewritten as

$\int_{t_{0}-t_{i}}^{+\infty}k_{T}(t)k_{T}(t+t_{i}-t_{j})dt.$ (69)

since $t\geq 0$ and $t_{j}-t_{i}-t\leq 0$ for all $t\in[t_{0}-t_{i},+\infty)$.
The form (69) will be used in every proof of this section.

####  D.2.1 Squared-Exponential Covariance Function

######  Lemma D.4.

Let $k_{T}$ be a Squared-Exponential covariance function (see Table 1), with
lengthscale $l_{S}>0$. Then,

$(k_{T}*k_{T})^{+\infty}_{t_{0}-t_{i}}(t_{i}-t_{j})=\frac{\sqrt{\pi}l_{T}}{2}e^{\frac{-(t_{i}-t_{j})^{2}}{2l_{T}^{2}}}\left(1-\text{erf}\left(\frac{2t_{0}-t_{i}-t_{j}}{2l_{T}}\right)\right)$
(70)

where erf is the error function.

###### Proof.

Since $k_{T}$ is a Squared-Exponential function, (69) becomes

$\displaystyle\int_{t_{0}-t_{i}}^{+\infty}e^{\frac{-t^{2}}{2l_{T}^{2}}}e^{\frac{-(t-t_{j}+t_{i})^{2}}{2l_{T}^{2}}}dt$
$\displaystyle=$
$\displaystyle\int_{t_{0}-t_{i}}^{+\infty}e^{\frac{-(2t^{2}-2(t_{j}-t_{i})t+t_{j}^{2}+t_{i}^{2}-2t_{i}t_{j})}{2l_{T}^{2}}}dt.$
(71)

It is known (see [29]) that

$\int
e^{-(ax^{2}+2bx+c)}dx=\frac{1}{2}\sqrt{\frac{\pi}{a}}e^{\frac{b^{2}-ac}{a}}\text{erf}\left(\sqrt{a}x+\frac{b}{\sqrt{a}}\right)$

where erf is the error function.

Therefore, (71) becomes

$\displaystyle\frac{\sqrt{\pi}l_{T}}{2}e^{\frac{-(t_{i}-t_{j})^{2}}{2l_{T}^{2}}}\left(1-\text{erf}\left(\frac{2t_{0}-t_{i}-t_{j}}{2l_{T}}\right)\right).$

∎

####  D.2.2 Matérn Covariance Function

######  Lemma D.5.

Let $k_{T}$ be a Matérn covariance function (see Table 1), with smoothness
parameter $\nu=p+\frac{1}{2},p\in\mathbb{N}$ and lengthscale $l_{T}>0$. Then,

$(k_{T}*k_{T})^{l_{T}+t_{0}-t_{i}}_{t_{0}-t_{i}}(t_{i}-t_{j})=\sum_{k_{1}=0}^{p}\sum_{k_{2}=0}^{p}C_{k_{1}k_{2}}e^{\frac{-\sqrt{2p+1}(2t_{0}-t_{i}-t_{j})}{l_{T}}}P_{k_{1}k_{2}}(t_{0},t_{i},t_{j})$
(72)

where

$\displaystyle C_{k_{1}k_{2}}$
$\displaystyle=\left(\frac{p!}{(2p)!}\right)^{2}\frac{(p+k_{1})!(p+k_{2})!}{k_{1}!k_{2}!(p-k_{1})!(p-k_{2})!}\left(\frac{2\sqrt{2p+1}}{l_{T}}\right)^{2p-k_{1}-k_{2}-1},$
(73) $\displaystyle P_{k_{1}k_{2}}(t_{0},t_{i},t_{j})$
$\displaystyle=\sum_{k_{3}=0}^{2p-k_{1}-k_{2}}\left(\frac{l_{T}}{2\sqrt{2p+1}}\right)^{k_{3}}P^{(k_{3})}(t_{0}-t_{i}),$
(74) $\displaystyle P(t)$ $\displaystyle=t^{p-k_{1}}(t-t_{j}+t_{i})^{p-k_{2}}$
(75)

and $P^{(k)}$ the $k$th derivative of $P(t)$ with respect to $t$.

###### Proof.

The Matérn covariance function has a simpler form when its smoothness
parameter $\nu$ is a half-integer, that is $\nu=p+\frac{1}{2},p\in\mathbb{N}$
(see [28]). In that case,

$k_{T}(t)=e^{\frac{-\sqrt{2p+1}t}{l_{T}}}\frac{p!}{(2p)!}\sum_{k_{1}=0}^{p}\frac{(p+k_{1})!}{k_{1}!(p-k_{1})!}\left(\frac{2\sqrt{2p+1}t}{l_{T}}\right)^{p-k_{1}}.$

Therefore,

$\displaystyle
k_{T}(t)k_{T}(t+t_{i}-t_{j})=\frac{2\sqrt{2p+1}}{l_{T}}\sum_{k_{1}=0}^{p}\sum_{k_{2}=0}^{p}C_{k_{1}k_{2}}e^{\frac{-\sqrt{2p+1}(2t-t_{j}+t_{i})}{l_{T}}}P(t)$
(76)

with $C_{k_{1}k_{2}}$ defined in (73) and $P(t)$ defined in (75).

Integrating (76), we get

$\displaystyle\frac{2\sqrt{2p+1}}{l_{T}}\sum_{k_{1}=0}^{p}\sum_{k_{2}=0}^{p}C_{k_{1}k_{2}}e^{\frac{-\sqrt{2p+1}(t_{i}-t_{j})}{l_{T}}}\int_{t_{0}-t_{i}}^{+\infty}e^{\frac{-2\sqrt{2p+1}t}{l_{T}}}P(t)dt$
(77)

thanks to the linearity of the integral.

It is known (see [29]) that

$\int
P(x)e^{ax}dx=\frac{e^{ax}}{a}\sum_{k=0}^{m}(-1)^{k}\frac{P_{m}^{(k)}(x)}{a^{k}}$

where $P_{m}$ is a polynomial of degree $m$ and $P_{m}^{(k)}$ is the $k$th
derivative of $P_{m}$.

Therefore,

$\displaystyle\int
e^{\frac{-2\sqrt{2p+1}t}{l_{T}}}P(t)dt=-\frac{l_{T}}{2\sqrt{2p+1}}e^{-\frac{2\sqrt{2p+1}t}{l_{T}}}\sum_{k_{3}=0}^{2p-k_{1}-k_{2}}\frac{P^{(k_{3})}(t)l_{T}^{k_{3}}}{\left(2\sqrt{2p+1}\right)^{k_{3}}}$
(78)

Combining (77) and (78) we get

$(k_{T}*k_{T})^{+\infty}_{t_{0}-t_{i}}(t_{i}-t_{j})=\sum_{k_{1}=0}^{p}\sum_{k_{2}=0}^{p}C_{k_{1}k_{2}}e^{\frac{-\sqrt{2p+1}(2t_{0}-t_{i}-t_{j})}{l_{T}}}P_{k_{1}k_{2}}(t_{0},t_{i},t_{j})$

with $C_{k_{1}k_{2}}$ defined in (73) and $P_{k_{1}k_{2}}$ defined in (74).

This concludes the proof. ∎

##  Appendix E Relative Quantification of Relevancy

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((\bm{x},t),(\bm{x}^{\prime},t^{\prime}))$ with the formulas provided in
Tables 3 and 4 for the SE covariance functions. We get

$\displaystyle C((\bm{x},t),(\bm{x}^{\prime},t^{\prime}))$
$\displaystyle=\pi^{\frac{d}{2}}l_{S}^{d}e^{\frac{-||\bm{x}-\bm{x}^{\prime}||_{2}^{2}}{4l_{S}^{2}}}\frac{\sqrt{\pi}}{2}l_{T}e^{\frac{-(t-t^{\prime})^{2}}{2l_{T}^{2}}}\left(1-\text{erf}\left(\frac{2t_{0}-t-t^{\prime}}{2l_{T}}\right)\right)$
$\displaystyle=\frac{1}{2}\pi^{\frac{d+1}{2}}l_{S}^{d}l_{T}e^{\frac{-||\bm{x}-\bm{x}^{\prime}||_{2}^{2}}{4l_{S}^{2}}-\frac{-(t-t^{\prime})^{2}}{2l_{T}^{2}}}\left(1-\text{erf}\left(\frac{2t_{0}-t-t^{\prime}}{2l_{T}}\right)\right)$
$\displaystyle=\frac{1}{2}\pi^{\frac{d+1}{2}}l_{S}^{d}l_{T}C^{*}((\bm{x},t),(\bm{x}^{\prime},t^{\prime})).$
(79)

The dependency on the covariance function hyperparameters
$\lambda,l_{S},l_{T}$ appears clearly in (79). Both $l_{S}$ and $l_{T}$ are
used, not only to scale the spatial distance $||\bm{x}-\bm{x}^{\prime}||_{2}$
and the temporal distance $|t-t^{\prime}|$ in $C^{*}$, but also as a scaling
constant of the magnitude of the output of the product of convolutions itself.
Because $C((\bm{x},t),(\bm{x}^{\prime},t^{\prime}))$ is involved in every term
of (10) and (11) in Theorem 4.1,
$\frac{1}{2}\pi^{\frac{d+1}{2}}l_{S}^{d}l_{T}$ can be factored out of (10) and
(11). Overall, both equations have in common the factor
$\frac{1}{2}\pi^{\frac{d+1}{2}}\lambda^{2}l_{S}^{d}l_{T}$. Clearly, this shows
how the covariance function hyperparameters
$\bm{\theta}=\left(\lambda,l_{S},l_{T}\right)$ 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
$\frac{1}{2}\pi^{\frac{d+1}{2}}\lambda^{2}l_{S}^{d}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 (79) no
longer exists.

##  Appendix F Removal Budget

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 form

$\begin{cases}b_{0}&=1,\\\ b_{t+\Delta t}&=b_{t}(1+\alpha)^{\Delta
t/l_{T}}\end{cases}$ (80)

with $l_{T}$ the temporal lengthscale (see Assumption 3.2) and $\alpha$ 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 (80).

More interestingly, let us discuss the exponent $\Delta t/l_{T}$. Arguably, an
alternative, more intuitive form of the removal budget would be

$\begin{cases}b_{0}&=1,\\\ b_{t+\Delta t}&=b_{t}(1+\alpha)^{\Delta
t}\end{cases}.$ (81)

Table 5: Comparison of removal budgets (80) and (81) when doing experiments of different durations on the Hartmann3d synthetic function. All experiments use the same time domain $[0,1]$. Duration $D$ |  $1$ Second | Lengthscale $l_{T}$ | Lengthscale $l_{T}$ | Budget (80) | Budget (81)
---|---|---|---|---|---
(seconds) | (axis unit) | (axis unit) | (seconds) | $(1+\alpha)^{D/l_{T}}$ | $(1+\alpha)^{D}$
300 | $1/300$ | $3/5$ | 180 | $(1+\alpha)^{5/3}$ | $(1+\alpha)^{300}$
600 | $1/600$ | $3/5$ | 480 | $(1+\alpha)^{5/3}$ | $(1+\alpha)^{600}$
1800 | $1/1800$ | $3/5$ | 1080 | $(1+\alpha)^{5/3}$ | $(1+\alpha)^{1800}$

Although easier to understand, the budget (81) 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 (81), which
depends on the number of elapsed seconds only, also varies. Conversely, the
budget (80) 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 (81) becomes really troublesome when it comes to
making a recommendation for the hyperparameter $\alpha$. If the analysis in
Section 5.1 had used the budget (81), its recommendation $\alpha^{*}$ 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 (80), is
a single number that is valid regardless of experiment duration. This is a
much more general insight.

##  Appendix G Empirical Results

###  G.1 Experimental Settings

In each experiment, the $d$-dimensional spatial domain is scaled in
$\mathcal{S}^{\prime}=[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 $\mathcal{S}^{\prime}\times\left[0,\frac{1}{40}\right]$.
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 time taken to perform the tasks (i) and
(ii) constitutes the response time of the DBO algorithm, denoted by $\Delta
t$. Clearly, $\Delta 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+\Delta t$ (except for ABO if it
decides to sample $f$ later).

The solutions have been implemented using BOTorch [30] (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 [31] (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).

###  G.2 Benchmarks and Figures

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.

Figure 6: (Left) Average response time and average regrets of the DBO
solutions during the optimization of the Rastrigin synthetic function. (Right)
Dataset sizes of the DBO solutions during the optimization of the Rastrigin
synthetic function.

##### Rastrigin.

The Rastrigin function is $d$-dimensional, and has the form

$f(\bm{x})=ad+\sum_{i=1}^{d}x_{i}^{2}-a\cos\left(2\pi x_{i}\right).$

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

Figure 7: (Left) Average response time and average regrets of the DBO
solutions during the optimization of the Schwefel synthetic function. (Right)
Dataset sizes of the DBO solutions during the optimization of the Schwefel
synthetic function.

##### Schwefel.

The Schwefel function is $d$-dimensional, and has the form

$f(\bm{x})=418.9829d-\sum_{i=1}^{d}x_{i}\sin\left(\sqrt{|x_{i}|}\right).$

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

Figure 8: (Left) Average response time and average regrets of the DBO
solutions during the optimization of the Styblinski-Tang synthetic function.
(Right) Dataset sizes of the DBO solutions during the optimization of the
Styblinski-Tang synthetic function.

##### Styblinski-Tang.

The Syblinski-Tang function is $d$-dimensional, and has the form

$f(\bm{x})=\frac{1}{2}\sum_{i=1}^{d}x_{i}^{4}-16x_{i}^{2}+5x_{i}.$

For the numerical evaluation, we set $d=4$ and we optimized the function on
the domain $[-5,5]^{d}$. The results are provided in Figure 8.

Figure 9: (Left) Average response time and average regrets of the DBO
solutions during the optimization of the Eggholder synthetic function. (Right)
Dataset sizes of the DBO solutions during the optimization of the Eggholder
synthetic function.

##### Eggholder.

The Eggholder function is $2$-dimensional, and has the form

$f(\bm{x})=-(x_{2}+47)\sin\left(\sqrt{|x_{2}+\frac{x_{1}}{2}+47|}\right)-x_{1}\sin\left(\sqrt{|x_{1}-x_{2}-47|}\right).$

For the numerical evaluation, we optimized the function on the domain
$[-512,512]^{2}$. The results are provided in Figure 9.

Figure 10: (Left) Average response time and average regrets of the DBO
solutions during the optimization of the Ackley synthetic function. (Right)
Dataset sizes of the DBO solutions during the optimization of the Ackley
synthetic function.

##### Ackley.

The Ackley function is $d$-dimensional, and has the form

$f(\bm{x})=-a\exp\left(-b\sqrt{\frac{1}{d}\sum_{i=1}^{d}x_{i}^{2}}\right)-\exp\left(\frac{1}{d}\sum_{i=1}^{d}\cos(cx_{i})\right)+a+\exp(1).$

For the numerical evaluation, we set $a=20$, $b=0.2$, $c=2\pi$, $d=4$ and we
optimized the function on the domain $[-32,32]^{d}$. The results are provided
in Figure 10.

Figure 11: (Left) Average response time and average regrets of the DBO
solutions during the optimization of the Rosenbrock synthetic function.
(Right) Dataset sizes of the DBO solutions during the optimization of the
Rosenbrock synthetic function.

##### Rosenbrock.

The Rosenbrock function is $d$-dimensional, and has the form

$f(\bm{x})=\sum_{i=1}^{d-1}100(x_{i+1}-x_{i}^{2})^{2}+(x_{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 11.

Figure 12: (Left) Average response time and average regrets of the DBO
solutions during the optimization of the Shekel synthetic function. (Right)
Dataset sizes of the DBO solutions during the optimization of the Shekel
synthetic function.

##### Shekel.

The Shekel function is $4$-dimensional, and has the form

$f(\bm{x})=-\sum_{i=1}^{m}\left(\sum_{j=1}^{4}(x_{j}-C_{ji})^{2}+\beta_{i}\right)^{-1}.$

For the numerical evaluation, we set $m=10$,
$\bm{\beta}=\frac{1}{10}\left(1,2,2,4,4,6,3,7,5,5\right)$,

$\bm{C}=\begin{pmatrix}4&1&8&6&3&2&5&8&6&7\\\ 4&1&8&6&7&9&3&1&2&3.6\\\
4&1&8&6&3&2&5&8&6&7\\\ 4&1&8&6&7&9&3&1&2&3.6\end{pmatrix},$

and we optimized the function on the domain $[0,10]^{4}$. The results are
provided in Figure 12.

Figure 13: (Left) Average response time and average regrets of the DBO
solutions during the optimization of the Hartmann-3 synthetic function.
(Right) Dataset sizes of the DBO solutions during the optimization of the
Hartmann-3 synthetic function.

##### Hartmann-3.

The Hartmann-3 function is $3$-dimensional, and has the form

$f(\bm{x})=-\sum_{i=1}^{4}\alpha_{i}\exp\left(-\sum_{j=1}^{3}A_{ij}(x_{j}-P_{ij})^{2}\right).$

For the numerical evaluation, we set
$\bm{\alpha}=\left(1.0,1.2,3.0,3.2\right)$,

$\bm{A}=\begin{pmatrix}3&10&30\\\ 0.1&10&35\\\ 3&10&30\\\ 0.1&10&35\\\
\end{pmatrix},\bm{P}=10^{-4}\begin{pmatrix}3689&1170&2673\\\ 4699&4387&7470\\\
1091&8732&5547\\\ 381&5743&8828\end{pmatrix},$

and we optimized the function on the domain $[0,1]^{3}$. The results are
provided in Figure 13.

Figure 14: (Left) Average response time and average regrets of the DBO
solutions during the optimization of the Hartmann-6 synthetic function.
(Right) Dataset sizes of the DBO solutions during the optimization of the
Hartmann-6 synthetic function.

##### Hartmann-6.

The Hartmann-6 function is $6$-dimensional, and has the form

$f(\bm{x})=-\sum_{i=1}^{4}\alpha_{i}\exp\left(-\sum_{j=1}^{6}A_{ij}(x_{j}-P_{ij})^{2}\right).$

For the numerical evaluation, we set
$\bm{\alpha}=\left(1.0,1.2,3.0,3.2\right)$,

$\bm{A}=\begin{pmatrix}10&3&17&3.50&1.7&8\\\ 0.05&10&17&0.1&8&14\\\
3&3.5&1.7&10&17&8\\\ 17&8&0.05&10&0.1&14\\\
\end{pmatrix},\bm{P}=10^{-4}\begin{pmatrix}1312&1696&5569&124&8283&5886\\\
2329&4135&8307&3736&1004&9991\\\ 2348&1451&3522&2883&3047&6650\\\
4047&8828&8732&5743&1091&381\end{pmatrix},$

and we optimized the function on the domain $[0,1]^{6}$. The results are
provided in Figure 14.

Figure 15: (Left) Average response time and average regrets of the DBO
solutions during the optimization of the Powell synthetic function. (Right)
Dataset sizes of the DBO solutions during the optimization of the Powell
synthetic function.

##### Powell.

The Powell function is $d$-dimensional, and has the form

$f(\bm{x})=\sum_{i=1}^{d/4}(x_{4i-3}+10x_{4i-2})^{2}+5(x_{4i-1}-x_{4i})^{2}+(x_{4i-2}-2x_{4i-1})^{4}+10(x_{4i-3}-x_{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 15.

Figure 16: (Left) Average response time and average regrets of the DBO
solutions during the Temperature real-world experiment. (Right) Dataset sizes
of the DBO solutions during the Temperature real-world experiment.

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

Figure 17: (Left) Average response time and average regrets of the DBO
solutions during the WLAN real-world experiment. (Right) Dataset sizes of the
DBO solutions during the WLAN real-world experiment.

##### 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 [32] 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}(\bm{x},t)=W\log_{2}(1+S_{ij}(\bm{x},t)).$

Then, the objective function is

$f(\bm{x},t)=\sum_{i=1}^{4}\sum_{j\in\mathcal{N}_{i}}C_{ij}(\bm{x},t),$

with $\mathcal{N}_{i}$ the end-users associated with node $i$.

For the numerical evaluation, we optimized the power levels $\bm{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 17.

###  G.3 Discussion on Empirical Performance

In this section, we discuss the performance achieved by all the DBO solutions
on the benchmarks introduced in the previous section.

##### 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., $\lambda,l_{S},\sigma^{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 3).

##### 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 11) and
Powell (see Figure 15). Conversely, ABO obtains competitive performance on
functions with smaller temporal lengthscales, e.g. Eggholder (see Figure 9) or
Shekel (see Figure 12). 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 $\sigma^{2}$) are inferred on the fly, the MLE explains the variance in
the observations with an increasingly large observational noise level
$\sigma^{2}$ as time goes by. However, by construction of ET-GP-UCB, the
greater $\sigma^{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 14),
Shekel (see Figure 12) or Ackley (see Figure 10). 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 15).

##### 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
$\epsilon$, and the reset is triggered each time the dataset size exceeds
$N(\epsilon)$ (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 9),
Hartmann3d (see Figure 13) or Temperature (see Figure 16).

##### TV-GP-UCB.

TV-GP-UCB directly accounts for temporal correlations by computing a specific
covariance function controlled by an hyperparameter $\epsilon$. 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 6), Schwefel (see Figure 7) or Shekel (see Figure 12).
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 15) or Rosenbrock (see Figure 11), 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 10), Shekel (see Figure 12) or Temperature (see Figure 16). 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 3), but also on most of the benchmarks themselves, e.g., Schwefel
(see Figure 7), Ackley (see Figure 10), Shekel (see Figure 12), Hartmann-6
(see Figure 14) or Powell (see Figure 15).

###  G.4 Animated Visualizations

In this section, we describe and discuss the two animated visualizations
provided as supplementary material for the paper. These videos show W-DBO
optimizing two 2-dimensional synthetic functions. They depict W-DBO’s
predictions, collected observations and deleted observations into the spatio-
temporal domains of the functions.

One of the videos depict the optimization of the Six-Hump Camel function on
the domain $[-2,2]^{2}$. The SHC function is

$\text{SHC}(x,t)=\left(4-2.1x^{2}+\frac{x^{4}}{3}\right)x^{2}+xt+\left(-4+4t^{2}\right)t^{2}.$

To study how W-DBO reacts to sudden changes in the objective function, the
other video depict the optimization of the piecewise function

$f(x,t)=\begin{cases}\text{SHC}(x,t)&\text{if }t<-\frac{1}{2},\\\
\text{SHC}(t,x)&\text{otherwise.}\end{cases}$

Figure 18: Snapshot from one of the videos showing the optimization conducted
by W-DBO. The normalized temporal dimension is shown on the x-axis and the
normalized spatial dimension is shown on the y-axis. The observations that are
in the dataset are depicted as red dots, while the deleted observations are
depicted as black crosses. The maximal arguments
$\left\\{\operatorname*{arg\,max}_{x\in\mathcal{S}}f(x,t),t\in\mathcal{T}\right\\}$
are depicted with a cyan curve. The predictions of W-DBO are shown with a
contour plot. Finally, the present time is depicted as a black vertical line
labelled $t_{0}$.

A snapshot from the latter can be found in Figure 18. It illustrates that the
benefits brought by W-DBO are substantial, since the algorithm is able to
track $\max_{x\in\mathcal{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).

##  Appendix H Limitations

For the sake of completeness, we explicitly discuss the limitations of W-DBO
in this appendix. Three 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.

  * •

We also introduce a structure for spatio-temporal correlations with Assumption
3.2. Although less restrictive than the one enforced by [21, 22] and
equivalent to the one in [20], 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.