# Incremental Ensemble Gaussian Processes

Qin Lu, Georgios V. Karanikolas, and Georgios B. Giannakis The authors are
with Dept. of Electrical and Computer Engineering and Digital Technology
Center, University of Minnesota Minneapolis, MN 55455. E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>The first two authors are equally
contributed.

###### Abstract

Belonging to the family of Bayesian nonparametrics, Gaussian process (GP)
based approaches have well-documented merits not only in learning over a rich
class of nonlinear functions, but also in quantifying the associated
uncertainty. However, most GP methods rely on a single preselected kernel
function, which may fall short in characterizing data samples that arrive
sequentially in time-critical applications. To enable online kernel
adaptation, the present work advocates an incremental ensemble (IE-) GP
framework, where an EGP meta-learner employs an ensemble of GP learners, each
having a unique kernel belonging to a prescribed kernel dictionary. With each
GP expert leveraging the random feature-based approximation to perform online
prediction and model update with scalability, the EGP meta-learner capitalizes
on data-adaptive weights to synthesize the per-expert predictions. Further,
the novel IE-GP is generalized to accommodate time-varying functions by
modeling structured dynamics at the EGP meta-learner and within each GP
learner. To benchmark the performance of IE-GP and its dynamic variant in the
adversarial setting where the modeling assumptions are violated, rigorous
performance analysis has been conducted via the notion of regret, as the norm
in online convex optimization. Last but not the least, online unsupervised
learning for dimensionality reduction is explored under the novel IE-GP
framework. Synthetic and real data tests demonstrate the effectiveness of the
proposed schemes.

###### Index Terms:

Gaussian processes, ensemble learning, online prediction, random features,
regret analysis

##  1 Introduction

Gaussian processes (GPs) cross-fertilize merits of kernel methods and Bayesian
models to benefit several learning tasks, including regression,
classification, ranking, and dimensionality reduction [1]. In GP-based
approaches, a Gaussian _prior_ is assumed over a learning function $f(\cdot)$
with covariance (kernel) capturing similarities among
$\\{f(\mathbf{x}_{t})\\}$ dependent on inputs $\\{\mathbf{x}_{t}\\}$. Given
observed outputs $\\{y_{t}\\}$ linked to the latent function $f(\cdot)$ via
the conditionally independent per-datum likelihood
$p(y_{t}|f(\mathbf{x}_{t}))$, Bayes rule produces the _posterior_ distribution
of $f(\cdot)$, based on which task-specific inference can be effected on the
unseen data. Besides learning functions with rich expressiveness, the Bayesian
framework of GP-based approaches further quantifies uncertainty of the
function estimate, which is of utmost importance in safety-critical
applications. For instance in medical diagnosis [2], human intervention would
be called for when machine operated decisions are accompanied by high
uncertainty.

In spite of the intriguing performance, applicability of plain-vanilla GPs in
the big data regime is discouraged by the cubic computational complexity in
the number of training samples [1]. To relieve the scalability issue, various
attempts have been made, including efficient numerical operation [3, 4], and
structured approximants of the kernel matrix [5, 6]. Of special interest to us
is the random feature (RF) based approach, which, leveraging the spectral
properties of stationary kernels, converts the nonparametric GP paradigm to a
parametric one [7, 8]. Such a parametric approach readily accommodates online
processing of data samples [9], which is necessitated in time-critical
applications. For instance, the detection of spam emails is performed on an
email-by-email basis in real time. Albeit accommodating online operation, the
performance of existing RF-based GP methods hinges on the single preselected
kernel, that may fall short in characterizing upcoming data. Henceforth,
online kernel adaptation is essential to real-time decision making.

On the theoretical horizon, to benchmark performance of online approaches,
analysis is usually conducted via the notion of regret, the norm in online
convex optimization [10] and online learning with experts [11], to combat with
the adversarial setting where the generative assumptions are violated.
Although several scalable GP approaches have been developed for the online
operation [12, 13], regret analysis has not been touch upon except for the
plain-vanilla GP [14].

In accordance with the aforementioned desiderata, the goal of the current work
is to pursue algorithmic developments of scalable GPs that could enable kernel
adaptation to cope with function dynamics in the online scenario, as well as
benchmark performance of the resultant approaches via regret-based analysis.

###  1.1 Related works

To contextualize the current contribution, the following existing works will
be outlined.

(Online) Scalable GPs. Approaches to effect scalability in GPs rely on
advanced numerical methods [3, 4], special kernel functions [15, 16], or low
rank approximants of the kernel matrix [5, 17, 6, 7]. A well-known low-rank
scheme summarizes the $T$ training samples via $q(\ll T)$ pseudo data with
inducing inputs that are employed for inference in the testing phase [5, 17,
6]. This global summary amounts to approximating the original GP prior with a
kernel matrix having low rank $q$, thus reducing the complexity of batch
computations to $\mathcal{O}(Tq^{2})$. Rather than the spatial sampling,
another less explored low-rank approach leverages spectral components of
shift-invariant kernels to yield the random feature (RF) based kernel
approximation [8]. Converting the nonparametric GP prior to a parametric one,
the resultant RF-based GP approaches can afford complexity comparable to the
inducing points-based approximants [7, 18]. To accommodate time-critical
applications, online scalable GP approaches have been developed relying on
stochastic optimization or online variational inference; see, e.g., [12, 13,
19]. Albeit ensuring scalability, these approaches rely on a single GP kernel,
which may limit expressiveness of the sought function. Also, theoretical
analyses that quantify the robustness of the online solvers to the adversarial
setting are largely unexplored.

Expert-based GPs. An ensemble of (local or distributed) GP experts, each
relying on a unique kernel to summarize a subset of the training samples, has
been leveraged to lower computational complexity, or/and account for
nonstationarity of the learning function. Depending on how data samples are
distributed and how predictions over experts are aggregated, well-known
examples include the naive-local-experts [20], product-of-experts [21, 22],
mixture-of-experts [23, 24], and most recently sum-product networks [25] based
approaches. In spite of the advantages enjoyed in different lines of work,
existing expert-based GP approaches operate in _batch_ mode, thus falling
short in dealing with time-critical applications that welcome online decision-
making.

(Online) Multi-kernel learning. Parallel to the probabilistic GP paradigm,
kernel-based learning has been pursued also in the _deterministic_ reproducing
kernel Hilbert space (RKHS). Faced with inscalability arising from abudance of
training data, kernel-based approaches also resort to low-rank approximants of
the kernel matrix, including the RF-based approximation [8]. Bypassing kernel
selection via cross validation, data-driven multi-kernel learning enjoys well-
documented performances; see, e.g., [26, 27]. To further accommodating online
operation, scalable kernel-based learning has been investigated for a single
[28] as well as for an ensemble of learners [29, 30]. Most recently, online
RF-based approaches based on an ensemble of RKHS learners have been reported
along with their regret-based performance for static and dynamic settings
[30].

GP latent variable model (LVM). Leveraging GPs to model the mapping from the
hidden low-dimensional input space to high-dimensional observations, GPLVMs
are established probabilistic approaches to nonlinear dimensionality reduction
[31]. Scalable GPLVMs have been devised by relying on inducing points-based
approximations in the batch setting [32, 33], as well as the variational and
online variants [33, 34, 35, 36]. However, to the best of our knowledge, RF-
based counterparts, in spite of the application in kernel principle component
analysis (PCA) (see, e.g. [37]), have not been touched upon in the realm of
GPLVMs. Regarding ensemble learning, a GPLVM scheme which can be broadly
categorized in this area is [35], where different from the proposed approach
the goal is to track the latent state of a dynamical system. Ensemble methods
are, nonetheless, more commonplace in the context of probabilistic PCA for
_linear_ dimensionality reduction; see [38] for the seminal work and [39] for
an online variant.

###  1.2  Contributions

Relative to the aforementioned past works, the present paper aims at bringing
together the fields of scalable GPs and online learning with expert advice
[11]. The pursuit lies in algorithmic development, as well as performance
analysis via the measure of regret to account for the violations of the
generative models. The detailed contributions are highlighted as follows.

  * c1)

Towards online kernel adaptation, the present work advocates an incremental
(I) approach based on a weighted ensemble (E) of GP learners with scalable RF-
based kernel approximations. The novel IE-GP learns the unknown function and
jointly adapts to the appropriate EGP kernel on-the-fly.

  * c2)

To cope with learning nonstationary functions, dynamic IE-GP variants have
been devised to capture structured dynamics at the EGP meta-learner and
individual GP learners via a hidden Markov model and the state-space models,
respectively.

  * c3)

To account for data being adversarially chosen in the online setting, the
performances of IE-GP and its dynamic variant are compared with some benchmark
functions with data in hindsight via static and switching regret analyses. In
both cases, the cumulative regrets over $T$ slots are of order
$\mathcal{O}(\log T)$, implying no regret on average.

  * c4)

Complementary to the supervised function learning task, online unsupervised
learning for dimensionality reduction (a.k.a. latent variable model) is
investigated under the proposed IE-GP paradigm.

  * c5)

Extensive experimental results are provided to validate the merits of the
proposed methods in regression, classification and dimensionality reduction
tasks.

Relative to the conference precursor [40], the novelty lies in the following
four aspects: 1) A switching (S) IE-GP approach is devised by modeling
dynamics at the EGP meta-learner via a first-order Markov chain; 2) The
performance of the proposed SIE-GP is analysed via the notion of switching
regret; 3) EGP-based online unsupervised learning with RFs for scalability is
further explored; 4) Experimental section has been significantly expanded via
the inclusion of classification and dimensionality reduction tests.

Notation. Scalars are denoted by lowercase, column vectors by bold lowercase,
and matrices by bold uppercase fonts. Superscripts $\leavevmode\nobreak\
{}^{\top}$ and $\leavevmode\nobreak\ {}^{-1}$ denote transpose, and matrix
inverse, respectively; while $\mathbf{0}_{N}$ stands for the $N\times 1$ all-
zero vector; and $\mathcal{N}(\mathbf{x};\bm{\mu},\mathbf{K})$ for the
probability density function (pdf) of a Gaussian random vector $\mathbf{x}$
with mean $\bm{\mu}$, and covariance matrix $\bf K$. Subscript
“${t+1|\mathbf{t}}$” signifies that prediction for slot $t+1$ relies on the
_batch_ of samples up to and including $t$, while “${t+1|t}$” stands for a
single-step predictor. $I(x)$ represents the indicator function, which is $1$
if $x$ is true, and $0$ otherwise.

##  2 Preliminaries and background

As a prelude to our online EGP approach that will also introduce context and
notation, this section deals with batch and scalable learning based on a
single GP.

###  2.1 Non-scalable batch GP-based learning

Given data $\\{\mathbf{x}_{\tau},y_{\tau}\\}$, the goal is to learn a function
$f(\cdot)$ that links the $d\times 1$ input $\mathbf{x}_{\tau}$ with the
scalar output $y_{\tau}$ as $\mathbf{x}_{\tau}\rightarrow
f(\mathbf{x}_{\tau})\rightarrow y_{\tau}$. Postulating $f$ with a GP prior as
$f\sim\mathcal{GP}(0,\kappa(\mathbf{x},\mathbf{x}^{\prime}))$, where
$\kappa(\cdot,\cdot)$ is a kernel function measuring pairwise similarity of
any two inputs, the joint prior pdf of function evaluations
$\mathbf{f}_{t}:=[f(\mathbf{x}_{1}),\ldots,f(\mathbf{x}_{t})]^{\top}$ at any
inputs
$\mathbf{X}_{t}:=\left[\mathbf{x}_{1},\ldots,\mathbf{x}_{t}\right]^{\top}$ is
Gaussian distributed as [1]

$p(\mathbf{f}_{t}|\mathbf{X}_{t})=\mathcal{N}(\mathbf{f}_{t};{\bf 0}_{t},{\bf
K}_{t})\ \ \ \forall t$ (1)

where ${\bf K}_{t}$ is a $t\times t$ covariance matrix with
$(\tau,\tau^{\prime})$th entry $[{\bf K}_{t}]_{\tau,\tau^{\prime}}={\rm
cov}(f(\mathbf{x}_{\tau}),f(\mathbf{x}_{\tau^{\prime}})):=\kappa(\mathbf{x}_{\tau},\mathbf{x}_{\tau^{\prime}})$.

To estimate $f$, we rely on the observed outputs
$\mathbf{y}_{t}:=[y_{1},\ldots,y_{t}]^{\top}$ that are linked with
$\mathbf{f}_{t}$ via the conditional likelihood
$p(\mathbf{y}_{t}|\mathbf{f}_{t},\mathbf{X}_{t})=\prod_{\tau=1}^{t}p(y_{\tau}|f(\mathbf{x}_{\tau}))$
that is assumed known. Through Bayes’ rule, the latter will yield the
posterior $p(\mathbf{f}_{t}|\mathbf{y}_{t},\mathbf{X}_{t})\propto
p(\mathbf{f}_{t}|\mathbf{X}_{t})p(\mathbf{y}_{t}|\mathbf{f}_{t},\mathbf{X}_{t})$.
For Gaussian process regression (GPR) the conditional likelihood is assumed
normal with mean $\mathbf{f}_{t}$ and covariance matrix $\sigma_{n}^{2}{\bf
I}_{t}$, that is,
$p(\mathbf{y}_{t}|\mathbf{f}_{t},\mathbf{X}_{t})=\mathcal{N}(\mathbf{y}_{t};\mathbf{f}_{t},\sigma_{n}^{2}{\bf
I}_{t})$, which along with the GP prior in (1) yields the Gaussian posterior
$p(\mathbf{f}_{t}|\mathbf{y}_{t},\mathbf{X}_{t})$. For non-Gaussian
likelihoods, sampling or approximate inference techniques will be called for
to carry out the analytically intractable posterior
$p(\mathbf{f}_{t}|\mathbf{y}_{t},\mathbf{X}_{t})$ [1].

Prediction with a single GP. Given training data
$\\{\mathbf{X}_{t},\mathbf{y}_{t}\\}$ and a new test input $\mathbf{x}_{t+1}$,
we have from (1) that $p(f(\mathbf{x}_{t+1})|\mathbf{f}_{t},\mathbf{X}_{t})$
is Gaussian with known mean and covariance. Together with the known posterior
$p(\mathbf{f}_{t}|\mathbf{y}_{t},\mathbf{X}_{t})$, the so-termed predictive
pdf of $f(\mathbf{x}_{t+1})$ can be obtained as [1]

$\displaystyle
p(f(\mathbf{x}_{t+1})|\mathbf{y}_{t},\mathbf{X}_{t})\\!=\\!\\!\int\\!\\!p(f(\mathbf{x}_{t+1})|\mathbf{f}_{t},\mathbf{X}_{t})p(\mathbf{f}_{t}|\mathbf{y}_{t},\mathbf{X}_{t})d\mathbf{f}_{t}$
(2)

which is generally non-Gaussian if
$p(\mathbf{f}_{t}|\mathbf{y}_{t},\mathbf{X}_{t})$ is non-Gaussian, and thus
necessitates Monte Carlo (MC) sampling to estimate it. Alternatively,
$p(\mathbf{f}_{t}|\mathbf{y}_{t},\mathbf{X}_{t})$ can be approximated by a
Gaussian, yielding a Gaussian approximation for (2) as well. Of course,
$p(f(\mathbf{x}_{t+1})|\mathbf{y}_{t},\mathbf{X}_{t})$ is Gaussian for GPR,
with its mean and covariance matrix available in closed form.

Using the pdf in (2) and the known $p(y_{t+1}|f(\mathbf{x}_{t+1}))$, it is
also possible to find the predictive pdf of $y_{t+1}$ as

$\displaystyle p(y_{t+1}|\mathbf{y}_{t},\mathbf{X}_{t+1})=$ (3)
$\displaystyle\qquad\quad\int\\!\\!p(y_{t+1}|f(\mathbf{x}_{t+1}))p(f(\mathbf{x}_{t+1})|\mathbf{y}_{t},\mathbf{X}_{t})df(\mathbf{x}_{t+1})$

which generally requires MC sampling or $p(f_{t+1}|\mathbf{y}_{t};$
$\mathbf{X}_{t+1})$ to be (at least approximately) Gaussian. Either way, (3)
yields the data predictive pdf that fully quantifies the uncertainty of
$y_{t+1}$.

Specifically for GPR, we have

$\displaystyle
p(y_{t+1}|\mathbf{y}_{t},\mathbf{X}_{t+1})=\mathcal{N}(y_{t+1};\hat{y}_{t+1|\mathbf{t}},\sigma_{t+1|\mathbf{t}}^{2})$
(4)

where the mean and the variance of the predictor are given by [1]

$\displaystyle\hat{y}_{t+1|\mathbf{t}}$
$\displaystyle=\mathbf{k}_{t+1}^{\top}(\mathbf{K}_{t}+\sigma_{n}^{2}\mathbf{I}_{t})^{-1}\mathbf{y}_{t}$
(5a) $\displaystyle\sigma_{t+1|\mathbf{t}}^{2}$
$\displaystyle=\\!\kappa(\mathbf{x}_{t+1},\mathbf{x}_{t+1}\\!)\\!-\\!\mathbf{k}_{t+1}^{\top}\\!(\mathbf{K}_{t}\\!+\\!\sigma_{n}^{2}\mathbf{I}_{t})^{-1}\mathbf{k}_{t+1}\\!+\\!\sigma_{n}^{2}$
(5b)

with
$\mathbf{k}_{t+1}:=[\kappa(\mathbf{x}_{1},\mathbf{x}_{t+1}),\ldots,\kappa(\mathbf{x}_{t},\mathbf{x}_{t+1})]^{\top}$.
Clearly, this GP predictor is not scalable, since the complexity
$\mathcal{O}(t^{3})$ for inverting the $t\times t$ matrix in (5) will become
prohibitively high as $t$ grows.

###  2.2 Scalable RF learning with a single GP

Various attempts have been made to effect scalability in GP-based learning;
see, e.g., [17, 6, 7]. Most existing approaches amount to summarizing the
training data via a much smaller number of pseudo data with inducing inputs,
thereby obtaining a low-rank approximant of $\mathbf{K}_{t}$ [17]. However,
finding the locations of these inducing inputs entails involved optimization
procedure. Targeting a low-rank approximant that bypasses such intricate
training practice, we rely here on a standardized shift-invariant
$\bar{\kappa}(\mathbf{x},\mathbf{x}^{\prime})=\bar{\kappa}(\mathbf{x}-\mathbf{x}^{\prime})$,
whose inverse Fourier transform is

$\displaystyle\bar{\kappa}(\mathbf{x}-\mathbf{x}^{\prime})=\int\pi_{\bar{\kappa}}(\mathbf{v})e^{j\mathbf{v}^{\top}(\mathbf{x}-\mathbf{x}^{\prime})}d\mathbf{v}:=\mathbb{E}_{\pi_{\bar{\kappa}}}\left[e^{j\mathbf{v}^{\top}(\mathbf{x}-\mathbf{x}^{\prime})}\right]$
(6)

where $\pi_{\bar{\kappa}}$ is the power spectral density (PSD), and the last
equality follows after normalizing so that $\pi_{\bar{\kappa}}(\mathbf{v})$
integrates to $1$, thus allowing one to view it as a pdf.

Since $\bar{\kappa}$ is real, the expectation in (6) is given by
$\mathbb{E}_{\pi_{\bar{\kappa}}}\left[\cos(\mathbf{v}^{\top}(\mathbf{x}-\mathbf{x}^{\prime}))\right]$,
which, upon drawing a sufficient number, say $n_{\text{RF}}$, of independent
and identically distributed (i.i.d.) samples
$\\{\mathbf{v}_{j}\\}_{j=1}^{n_{\text{RF}}}$ from
$\pi_{\bar{\kappa}}(\mathbf{v})$, can be approximated by 111Quantities with
$\check{}$ involve RF approximations.

$\displaystyle\check{\bar{\kappa}}(\mathbf{x},\mathbf{x}^{\prime}):=\frac{1}{n_{\text{RF}}}\sum_{j=1}^{n_{\text{RF}}}\cos\left(\mathbf{v}_{j}^{\top}(\mathbf{x}-\mathbf{x}^{\prime})\right)\;.$
(7)

Define the $2n_{\text{RF}}\times 1$ random feature (RF) vector as [7]

$\displaystyle{\mbox{\boldmath$\phi$}}_{\mathbf{v}}(\mathbf{x})$ (8)
$\displaystyle:=\\!\frac{1}{\sqrt{n_{\text{RF}}}}\\!\left[\sin(\mathbf{v}_{1}^{\top}\mathbf{x}),\cos(\mathbf{v}_{1}^{\top}\mathbf{x}),\ldots,\sin(\mathbf{v}_{n_{\text{RF}}}^{\top}\mathbf{x}),\cos(\mathbf{v}_{n_{\text{RF}}}^{\top}\mathbf{x})\right]^{\top}$

which allows us to rewrite $\check{\bar{\kappa}}$ in (7) with
$\check{\bar{\kappa}}(\mathbf{x},\mathbf{x}^{\prime})={\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{\top}(\mathbf{x}){\mbox{\boldmath$\phi$}}_{\mathbf{v}}(\mathbf{x}^{\prime})$;
and thus, the parametric approximant

$\displaystyle{\check{f}}(\mathbf{x})={\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{\top}(\mathbf{x}){\mbox{\boldmath$\theta$}},\quad{\mbox{\boldmath$\theta$}}\sim\mathcal{N}({\mbox{\boldmath$\theta$}};\mathbf{0}_{2n_{\text{RF}}},\sigma_{\theta}^{2}\mathbf{I}_{2n_{\text{RF}}})$
(9)

can be viewed as coming from a realization of the Gaussian $\theta$ combined
with ${\mbox{\boldmath$\phi$}}_{\mathbf{v}}$ to yield the GP prior in (1) with
$\kappa=\sigma_{\theta}^{2}\bar{\kappa}$, where $\sigma_{\theta}^{2}$ is the
magnitude of $\kappa$. Clearly, for any $\mathbf{X}_{t}$, the prior pdf of
$\check{\mathbf{f}}_{t}$ is then

$\displaystyle
p(\check{\mathbf{f}}_{t}|\mathbf{X}_{t})=\mathcal{N}(\check{\mathbf{f}}_{t};\mathbf{0}_{t},\check{\mathbf{K}}_{t}),\quad\check{\mathbf{K}}_{t}=\sigma_{\theta}^{2}{\mbox{\boldmath$\Phi$}}_{t}{\mbox{\boldmath$\Phi$}}_{t}^{\top}$
(10)

where
${\mbox{\boldmath$\Phi$}}_{t}:=\left[{\mbox{\boldmath$\phi$}}_{\mathbf{v}}(\mathbf{x}_{1}),\ldots,{\mbox{\boldmath$\phi$}}_{\mathbf{v}}(\mathbf{x}_{t})\right]^{\top}$,
and $\check{\mathbf{K}}_{t}$ is then a low rank ($2n_{\text{RF}}$) approximant
of $\mathbf{K}_{t}$ in (1) for $t>2n_{\text{RF}}$.

With the parametric form of ${\check{f}}(\mathbf{x})$ in (9), the likelihood
$p(\mathbf{y}_{t}|\check{\mathbf{f}}_{t},\mathbf{X}_{t})$ is also parametrized
by $\theta$. This together with the Gaussian prior of $\theta$ (cf. (9)),
yields the posterior
$p({\mbox{\boldmath$\theta$}}|\mathbf{y}_{t},\mathbf{X}_{t})$, based on which
we can predict $f$ and $y$ at new test input $\mathbf{x}$. Specifically, upon
replacing $p(f(\mathbf{x}_{t+1})|\mathbf{f}_{t},\mathbf{X}_{t})$ and
$p(\mathbf{f}_{t}|\mathbf{y}_{t},\mathbf{X}_{t})$ in (2) by
$p(\check{f}(\mathbf{x}_{t+1})|{\mbox{\boldmath$\theta$}})=\delta(\check{f}(\mathbf{x}_{t+1})-{\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{\top}(\mathbf{x}_{t+1}){\mbox{\boldmath$\theta$}})$
and $p({\mbox{\boldmath$\theta$}}|\mathbf{y}_{t},\mathbf{X}_{t})$,
respectively, we obtain the predictive pdf of the RF-based
$\check{f}(\mathbf{x}_{t+1})$, which further leads to the predictive pdf of
$y_{t+1}$ in (3) after replacing $f(\mathbf{x}_{t+1})$ by
$\check{f}(\mathbf{x}_{t+1})$. For GPR, the predictive pdf of $y_{t+1}$ is

$\displaystyle{p}(y_{t+1}|\mathbf{y}_{t},\mathbf{X}_{t+1})=\mathcal{N}(y_{t+1};\hat{\check{y}}_{t+1|\mathbf{t}},\check{\sigma}^{2}_{t+1|\mathbf{t}})$
(11)

where

$\displaystyle\hat{\check{y}}_{t+1|\mathbf{t}}$
$\displaystyle={\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{\top}(\mathbf{x}_{t+1})\left({\mbox{\boldmath$\Phi$}}_{t}^{\top}{\mbox{\boldmath$\Phi$}}_{t}+\frac{\sigma_{n}^{2}}{\sigma_{\theta}^{2}}\mathbf{I}_{2n_{\text{RF}}}\right)^{-1}\\!\\!\\!\\!{\mbox{\boldmath$\Phi$}}_{t}^{\top}\mathbf{y}_{t}$
(12a) $\displaystyle\check{\sigma}^{2}_{t+1|\mathbf{t}}$
$\displaystyle={\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{\top}(\mathbf{x}_{t+1})\\!\\!\left(\\!\\!\frac{{\mbox{\boldmath$\Phi$}}_{t}^{\top}{\mbox{\boldmath$\Phi$}}_{t}}{\sigma_{n}^{2}}\\!+\\!\frac{\mathbf{I}_{2n_{\text{RF}}}}{\sigma_{\theta}^{2}}\\!\\!\right)^{-1}\\!\\!\\!\\!\\!\\!\\!\\!{\mbox{\boldmath$\phi$}}_{\mathbf{v}}(\mathbf{x}_{t+1})\\!+\\!\sigma_{n}^{2}.$
(12b)

This _batch_ predictor incurs complexity
$\mathcal{O}(t(2n_{\text{RF}})^{2}+(2n_{\text{RF}})^{3})$, which is dominated
by $\mathcal{O}(t(2n_{\text{RF}})^{2})$ for $t\gg 2n_{\text{RF}}$. This linear
(in $t$) complexity is apparently much more affordable than the plain-vanilla
GP predictor (5).

The RF-based function approximant $\check{f}$ easily accommodates _online_
operation [9], which is called for in many time-critical applications,
including time series prediction [41], and robot localization [42]. While the
RF-based online approach for GPR is offered in [9], its performance hinges on
a _preselected_ kernel for the GP prior, which may fall short in
characterizing upcoming data samples. Next, we will broaden the scope of a
single GP prior by an ensemble (E) of GPs to enable real-time kernel
adaptation. Besides serving the role of a non-Gaussian prior, EGP will turn
out to be scalable too, after adopting once again the RF approximation.

##  3 Online scalable ensemble GPs

Towards data-driven kernel selection in the online setting, an EGP meta-
learner employs an ensemble of $M$ GP experts (a.k.a. models or learners),
each of which places a unique GP prior on $f$ as
$f|m\sim\mathcal{GP}(0,\kappa^{m}(\mathbf{x},\mathbf{x}^{\prime}))$, where
$m\in\mathcal{M}:=\\{1,\ldots,M\\}$ is the expert index and $\kappa^{m}$ is a
shift-invariant kernel selected from a known kernel dictionary
$\mathcal{K}:=\\{\kappa^{1},\ldots,\kappa^{M}\\}$. Here, $\mathcal{K}$ should
be constructed as large as computational constraints allow, depending on
resources and the learning task. Per expert $m$, the prior pdf of function
values at $\mathbf{X}_{t}$ is

$\displaystyle
p(\mathbf{f}_{t}|i\\!=\\!m,\mathbf{X}_{t})\\!=\\!\mathcal{N}(\mathbf{f}_{t};{\bf
0}_{t},{\bf K}^{m}_{t}),\ [{\bf
K}_{t}^{m}]_{\tau,\tau^{\prime}}\\!:=\\!\kappa^{m}(\mathbf{x}_{\tau},\mathbf{x}_{\tau^{\prime}}\\!)$

where the hidden random variable $i$ is introduced to denote the expert index.
The ensemble prior pdf of $\mathbf{f}_{t}$ that accounts for all GP experts is
given by the Gaussian mixture (GM)

$\displaystyle
p(\mathbf{f}_{t}|\mathbf{X}_{t})=\sum_{m=1}^{M}w^{m}\mathcal{N}(\mathbf{f}_{t};{\bf
0}_{t},\mathbf{K}^{m}_{t})\;,\;\;\;\;\;\;\sum_{m=1}^{M}w^{m}=1$ (13)

where the unknown weights $\\{w^{m}\\}_{m=1}^{M}$, viewed as probabilities of
the GP experts to be present in the EGPs, are to be learned from data that
arrive sequentially.

Seeking a scalable predictor, each expert $m$ relies on the RF-based function
approximant (9) with the per-expert parameter vector
${\mbox{\boldmath$\theta$}}^{m}$ and RF vector
${\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{m}(\mathbf{x})$ constructed as in (8)
using $\\{\mathbf{v}_{j}^{m}\\}_{j=1}^{n_{\text{RF}}}$. Vectors
$\\{\mathbf{v}_{j}^{m}\\}_{j=1}^{n_{\text{RF}}}$ here are drawn i.i.d. from
$\pi_{\bar{\kappa}}^{m}(\mathbf{v})$, which is the PSD of the standardized
kernel $\bar{\kappa}^{m}$, relating to $\kappa^{m}$ through the magnitude
$\sigma_{\theta^{m}}^{2}$ as
$\kappa^{m}=\sigma_{\theta^{m}}^{2}\bar{\kappa}^{m}$. The per expert $m$
generative model for output $y$ is then

$\displaystyle{p}({\mbox{\boldmath$\theta$}}^{m})$
$\displaystyle=\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\mathbf{0}_{2n_{\text{RF}}},\sigma_{\theta^{m}}^{2}\mathbf{I}_{2n_{\text{RF}}})$
(14a) $\displaystyle p(y|{\mbox{\boldmath$\theta$}}^{m},\mathbf{x})$
$\displaystyle=p(y|{\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{m\top}(\mathbf{x}){\mbox{\boldmath$\theta$}}^{m})\;.$
(14b)

Focusing on the incremental (I) setting, each expert $m$ interleaves
prediction of $y_{t+1}$ based on
$p({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})$, and update
of the parameter posterior upon the arrival of $y_{t+1}$ per slot. To assess
the per-expert contribution, the EGP meta-learner relies on the posterior
probability $w_{t}^{m}:={\rm Pr}(i=m|\mathbf{y}_{t};\mathbf{X}_{t})$. As we
shall see next, the resultant IE-GP proceeds in two steps, namely prediction
and correction, by propagating per-expert weights and posterior pdfs
$\\{w_{t}^{m},{p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})\\}_{m=1}^{M}$
from slot to slot.

Prediction. Upon receiving $\mathbf{x}_{t+1}$, each expert $m$ constructs the
RF vector using ${\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{m}(\mathbf{x}_{t+1})$
as in (8). With
${p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})$ available
from slot $t$, the per-expert predictive pdf of $y_{t+1}$ can be obtained by
invoking the sum-product probability rule

$\displaystyle{p}(y_{t+1}|\mathbf{y}_{t},i=m,\mathbf{X}_{t+1})\\!\\!=\\!\\!\int\\!\\!p(y_{t+1}|{\mbox{\boldmath$\theta$}}^{m}\\!,\mathbf{x}_{t})p({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})d{\mbox{\boldmath$\theta$}}^{m}.$
(15)

Leveraging again the sum-product rule, the EGP meta-learner seeks the ensemble
predictive pdf as

$\displaystyle{p}(y_{t+1}|\mathbf{y}_{t},\mathbf{X}_{t+1})\\!$
$\displaystyle=\\!\\!\sum_{m=1}^{M}\\!\\!{\rm
Pr}(i\\!=\\!m|\mathbf{y}_{t},\\!\mathbf{X}_{t}){p}(y_{t+1}\\!|\mathbf{y}_{t},\\!i\\!=\\!m,\\!\mathbf{X}_{t+1}\\!)$
$\displaystyle=\\!\\!\sum_{m=1}^{M}w_{t}^{m}{p}(y_{t+1}|\mathbf{y}_{t},\\!i\\!=\\!m,\mathbf{X}_{t+1})$
(16)

which takes an intuitive form as a weighted combination of predictions from
the individual GP experts. Having available the predictive pdf, we are ready
to update the posterior pdf of the RF model parameter vector.

Correction. With the arrival of $y_{t+1}$, each expert $m$ updates the
posterior pdf of ${\mbox{\boldmath$\theta$}}^{m}$ via Bayes’ rule as

${p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t+1},\mathbf{X}_{t+1})=\frac{{p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t}\\!){p}(y_{t+1}|{\mbox{\boldmath$\theta$}}^{m},\mathbf{x}_{t+1})}{{p}(y_{t+1}|\mathbf{y}_{t},i=m,\mathbf{X}_{t+1})}$
(17)

where ${p}(y_{t+1}|{\mbox{\boldmath$\theta$}}^{m},\mathbf{x}_{t+1})$ is the
known likelihood (cf. (14b)), and
${p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})$ is
available from slot $t$. For later use, the so-termed Bayesian loss incurred
by expert $m$ at slot $t+1$ is (cf. [43])

$\displaystyle
l_{t+1|t}^{m}:=-\log{p}(y_{t+1}|\mathbf{y}_{t},i=m,\mathbf{X}_{t+1})$ (18)

whose ensemble version is given by

$\displaystyle\ell_{t+1|t}$
$\displaystyle:=-\log{p}(y_{t+1}|\mathbf{y}_{t},\mathbf{X}_{t+1})$
$\displaystyle\
=-\log\\!\\!\sum_{m=1}^{M}\\!\\!w_{t}^{m}\exp\\!\left(-l_{t+1|t}^{m}\\!\right).$
(19)

Simultaneously, the EGP meta-leaner obtains the updated weight
$w_{t+1}^{m}:={\rm Pr}(i=m|\mathbf{y}_{t+1},\mathbf{X}_{t+1})$ as

$\displaystyle w_{t+1}^{m}$ $\displaystyle=\frac{{\rm
Pr}(i=m|\mathbf{y}_{t},\mathbf{X}_{t}){p}(y_{t+1}|\mathbf{y}_{t},i=m,\mathbf{X}_{t+1})}{{p}(y_{t+1}|\mathbf{y}_{t},\mathbf{X}_{t+1})}$
$\displaystyle=w_{t}^{m}\exp(\ell_{t+1|t}-l_{t+1|t}^{m})$ (20)

where $w_{t}^{m}$ is available from slot $t$. Intuitively, large
$l_{t+1|t}^{m}$ implies small $\ell_{t+1|t}-l_{t+1|t}^{m}$, and thus
$w_{t+1}^{m}$ relative to the rest will be smaller than that at slot $t$.

Summarizing, our scalable IE-GP algorithm for general likelihoods (and thus
posteriors) relies on (15)-(17) to transition from slot $t$ to slot $t+1$.
Next, we specialize our novel IE-GP to GPR that enjoys closed-form pdf and
weight updates, as well as non-Gaussian likelihoods that entail Laplace
approximation [44] to evaluate the (possibly high-dimensional) integrals in
the prediction and correction steps.

Algorithm 1 IE-GP for GPR

1:Input: $\kappa^{m}$, $m=1,\ldots,M$, and number of RFs $n_{\text{RF}}$.

2:Initialization:

3:for $m=1,2,\ldots,M$ do

4: Draw $n_{\text{RF}}$ random vectors
$\\{\mathbf{v}_{i}^{m}\\}_{i=1}^{n_{\text{RF}}}$;

5: $w_{0}^{m}=1/M$;
$\hat{{\mbox{\boldmath$\theta$}}}_{0}^{m}=\mathbf{0}_{2D}$;
${\mbox{\boldmath$\Sigma$}}_{0}^{m}=\sigma_{\theta^{m}}^{2}\mathbf{I}_{2D}$;

6:end for

7:for $t=1,2,\ldots,T$ do

8: Receive input datum $\mathbf{x}_{t}$ ;

9: for $m=1,2,\ldots,M$ do

10: Construct RF ${\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{m}(\mathbf{x}_{t})$
via (8);

11: Obtain per-expert pdf of $y_{t}$ via (22);

12: Update $w_{t}^{m}$ via (25);

13: Update per-expert pdf of ${\mbox{\boldmath$\theta$}}^{m}$ via (26);

14: end for

15:end for

###  3.1 Closed-form updates for GPR

For GPR, the likelihood per expert is given by
${p}(y_{t}|{\mbox{\boldmath$\theta$}}^{m},\mathbf{x}_{t})=\mathcal{N}(y_{t};{\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{m\top}(\mathbf{x}_{t}){\mbox{\boldmath$\theta$}}^{m},\sigma_{n}^{2})$,
which together with the per-expert Gaussian prior
$p({\mbox{\boldmath$\theta$}}^{m})$ (cf. (14a)), yields the Gaussian posterior
at the end of slot $t$ expressed as

$\displaystyle{p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})=\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m},{\mbox{\boldmath$\Sigma$}}^{m}_{t})$
(21)

with mean $\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m}$ and covariance matrix
${\mbox{\boldmath$\Sigma$}}^{m}_{t}$ per expert $m$.

Building on (21) and (15), the predictive pdf of $y_{t+1}$ from expert $m$ is
also Gaussian

$\displaystyle{p}(y_{t+1}|\mathbf{y}_{t},i\\!=\\!m,\mathbf{X}_{t+1})\\!=\\!\mathcal{N}\left(y_{t+1};\hat{y}_{t+1|t}^{m},(\sigma_{t+1|t}^{m})^{2}\right)$
(22)

where the predicted mean and variance are

$\displaystyle\hat{y}_{t+1|t}^{m}$
$\displaystyle={\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{t+1})\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m}$
(23a) $\displaystyle(\sigma_{t+1|t}^{m})^{2}$
$\displaystyle={\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{t+1}){\mbox{\boldmath$\Sigma$}}^{m}_{t}{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{t+1})+\sigma_{n}^{2}\;.$
(23b)

Thus, the ensemble predictive pdf of $y_{t+1}$ in (16) specialized to GPR is a
GM, based on which the EGP meta-learner obtains the minimum mean-square error
(MMSE) predictor of $y_{t+1}$ together with the associated variance as

$\displaystyle\hat{y}_{t+1|t}$
$\displaystyle=\sum_{m=1}^{M}w_{t}^{m}\hat{y}_{t+1|t}^{m}$ (24a)
$\displaystyle\sigma_{t+1|t}^{2}$
$\displaystyle=\sum_{m=1}^{M}w_{t}^{m}[(\sigma_{t+1|t}^{m})^{2}+(\hat{y}_{t+1|t}\\!-\\!\hat{y}_{t+1|t}^{m})^{2}].$
(24b)

When $y_{t+1}$ becomes available, the EGP meta-learner updates the per-expert
weight as (cf. (20) and (22))

$\displaystyle
w_{t+1}^{m}=\frac{w_{t}^{m}\mathcal{N}\left(y_{t+1};\hat{y}_{t+1|t}^{m},(\sigma_{t+1|t}^{m})^{2}\right)}{\sum_{m^{\prime}=1}^{M}w_{t}^{m^{\prime}}\mathcal{N}\left(y_{t+1};\\!\hat{y}_{t+1|t}^{m^{\prime}},(\sigma_{t+1|t}^{m^{\prime}})^{2}\right)}\;.$
(25)

With the per-expert Gaussian likelihood, the arrival of $y_{t+1}$ also
propagates Gaussianity of the posterior pdf of
${\mbox{\boldmath$\theta$}}^{m}$ from slot $t$ to $t+1$, expressed as

$\displaystyle{p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t+1},\mathbf{X}_{t+1})$
$\displaystyle=\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\hat{{\mbox{\boldmath$\theta$}}}_{t+1}^{m},{\mbox{\boldmath$\Sigma$}}^{m}_{t+1})$
(26)

where the per-expert mean $\hat{{\mbox{\boldmath$\theta$}}}_{t+1}^{m}$ and
covariance matrix ${\mbox{\boldmath$\Sigma$}}^{m}_{t+1}$ are

$\displaystyle\hat{{\mbox{\boldmath$\theta$}}}_{t+1}^{m}$
$\displaystyle=\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m}\\!+\\!(\sigma_{t+1|t}^{m})^{-2}{\mbox{\boldmath$\Sigma$}}^{m}_{t}{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{t+1})(y_{t+1}\\!-\\!\hat{y}_{t+1|t}^{m})$
(27a) $\displaystyle{\mbox{\boldmath$\Sigma$}}_{t+1}^{m}$
$\displaystyle={\mbox{\boldmath$\Sigma$}}_{t}^{m}\\!\\!-\\!(\sigma_{t+1|t}^{m})^{-2}{\mbox{\boldmath$\Sigma$}}^{m}_{t}{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}\\!(\mathbf{x}_{t+1}){\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}\\!\\!(\mathbf{x}_{t+1}){\mbox{\boldmath$\Sigma$}}^{m}_{t}.$
(27b)

Accounting for all $M$ expert updates, our scalable IE-GP approach to GPR (see
Algorithm 1) has per-iteration complexity of
$\mathcal{O}(M(2n_{\text{RF}})^{2})$; hence, scalability is not compromised by
the ensemble approach that also offers a richer model for the learning
function. The deterministic RKHS online approach (termed “Raker” in [30])
relies on first-order gradient descent to update $\theta$ at per-iteration
complexity of $\mathcal{O}(Mn_{\text{RF}}d)$, which is lower than our second-
order update in (27). Our probabilistic IE-GP approach offers numerically
improved performance that is also analytically quantifiable through the
predictor variance (24b) in (24a).

###  3.2 Laplace approximation for non-Gaussian likelihood

When the likelihood is non-Gaussian as in classification, Poisson regression
or ordinal regression, the per-expert parameter posterior is no longer
available in closed form. Fortunately, the so-termed Laplace approximation
[44] can be leveraged to carry out the prediction and correction steps with
tractability in IE-GP. Specifically, each expert maintains a Gaussian
approximant for the parameter vector as
$p({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})\approx\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\hat{{\mbox{\boldmath$\theta$}}}^{m}_{t},{\mbox{\boldmath$\Sigma$}}_{t}^{m})$,
which readily allows the predicted pdf of $y_{t+1}$ (15) to be calculated via
Monte-Carlo sampling, or through probit approximation for binary
classification with the logistic likelihood; see, e.g., [45, Chapter 8.4.4.2].

Given $y_{t+1}$, the EGP meta-learner updates the per-expert weight as in (20)
with the per-expert loss at hand from the prediction step. Also, each expert
relies on Laplace approximation to seek an updated Gaussian approximant as
${p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t+1},\mathbf{X}_{t+1})\approx\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\hat{{\mbox{\boldmath$\theta$}}}_{t+1}^{m},{\mbox{\boldmath$\Sigma$}}^{m}_{t+1})$,
where $\hat{{\mbox{\boldmath$\theta$}}}^{m}_{t+1}$ and
${\mbox{\boldmath$\Sigma$}}^{m}_{t+1}$ are respectively the mode of the log
posterior and the corresponding information matrix inverse, that can be
obtained by solving the following optimization problem using Newton’s
iteration [45]

$\displaystyle\hat{{\mbox{\boldmath$\theta$}}}_{t+1}^{m}=\underset{{\mbox{\boldmath$\theta$}}^{m}}{\arg\max}\log{p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})+\log
p(y_{t+1}|{\mbox{\boldmath$\theta$}}^{m},\mathbf{x}_{t+1})$
$\displaystyle({\mbox{\boldmath$\Sigma$}}^{m}_{t+1})^{-1}=({\mbox{\boldmath$\Sigma$}}^{m}_{t})^{-1}\\!\\!-\\!\nabla^{2}_{\theta^{m}}\log
p(y_{t+1}|{\mbox{\boldmath$\theta$}}^{m}\\!\\!,\mathbf{x}_{t+1})\Big{\rvert}_{{\mbox{\boldmath$\theta$}}^{m}=\hat{{\mbox{\boldmath$\theta$}}}_{t+1}^{m}}.$

##  4 Regret analysis

The pdf ${p}(y_{t+1}|\mathbf{y}_{t},\mathbf{X}_{t+1})$ in (16) provide an
online performance metric for ${\hat{y}}_{t+1|t}$, from which its mean and
variance can be also obtained (even in closed form, cf. (24b)). These metrics
however, rely on the assumption of knowing the prior pdf of $f$, and the
conditional data likelihood. To guard against having _imperfect knowledge_ of
these pdfs (the norm in adversarial settings), regret analysis is well
motivated along the lines of online convex optimization [10] and online
learning with expert advice [11]. This is the subject of this section that
aims to benchmark performance of our IE-GP predictor relative to the best
function estimator with data in hindsight when the generative assumptions are
violated.

To this end, let $\mathcal{L}(f(\mathbf{x}_{\tau});y_{\tau}):=-\log
p(y_{\tau}|f(\mathbf{x}_{\tau}))$ be the per-slot negative log-likelihood
(NLL). For any fixed function estimator $\hat{f}^{*}(\cdot)$, the incurred
loss over $T$ slots is
$\sum_{\tau=1}^{T}\mathcal{L}({\hat{f}}^{*}(\mathbf{x}_{\tau});y_{\tau})$.
With the EGP prior in (13), the best function estimate (benchmark) with data
$\\{\mathbf{X}_{T},\mathbf{y}_{T}\\}$ available in hindsight, are obtained
with the optimal weights $\\{w^{m}\\}$ in the EGP prior by maximizing the
batch function posterior,
$p(\mathbf{f}_{T}|\mathbf{y}_{T},\mathbf{X}_{T})\propto
p(\mathbf{f}_{T}|\mathbf{X}_{T})p(\mathbf{y}_{T}|\mathbf{f}_{T},\mathbf{X}_{T})$,
as

$\displaystyle(\hat{\mathbf{f}}_{T},\\{\hat{w}^{m}\\})=\underset{\begin{subarray}{c}\mathbf{f}_{T},\\{w^{m}\\}\\\
\sum_{m}w^{m}=1\end{subarray}}{\arg\max}p(\mathbf{y}_{T}|\mathbf{f}_{T},\mathbf{X}_{T})\\!\\!\sum_{m=1}^{M}\\!w^{m}p(\mathbf{f}_{T}|i\\!=\\!m,\mathbf{X}_{T})$

whose solution is $\hat{w}^{m^{*}}=1$ and $\hat{w}^{m}=0$ for $m\neq m^{*}$.
This implies that only one GP expert $m^{*}$ is active in the benchmark
function estimate for $\tau=1,\ldots,T$. The optimal estimate by expert
$m^{*}$ are then given by

$\displaystyle\hat{\mathbf{f}}_{T}=\underset{\mathbf{f}_{T}}{\arg\max}\ \
p(\mathbf{f}_{T}|i=m^{*},\mathbf{X}_{T})p(\mathbf{y}_{T}|\mathbf{f}_{T},\mathbf{X}_{T})\;.$
(29)

As every positive semidefinite kernel $\kappa^{m}$ is associated with a unique
RHKS $\mathcal{H}^{m}$[14], the optimal function estimator
$\hat{f}^{m^{*}}(\cdot)$ is extracted from (29) as

$\displaystyle m^{*}\in$
$\displaystyle\>\underset{m\in\mathcal{M}}{\arg\min}\sum_{\tau=1}^{T}{\mathcal{L}}({\hat{f}}^{m}(\mathbf{x}_{\tau});y_{\tau})+\frac{1}{2}\|{\hat{f}}^{m}\|_{\mathcal{H}^{m}}^{2}$
(30)

where the optimal function estimator per expert $\hat{f}^{m}(\cdot)$,
$m=1,\ldots,M$, is obtained as

$\displaystyle{\hat{f}}^{m}(\cdot)\in$
$\displaystyle\>\underset{f^{m}\in\mathcal{H}^{m}}{\arg\min}\sum_{\tau=1}^{T}{\mathcal{L}}(f^{m}(\mathbf{x}_{\tau});y_{\tau})+\frac{1}{2}\|f^{m}\|_{\mathcal{H}^{m}}^{2}\>.$

With the best fixed function estimator $\hat{f}^{m^{*}}(\cdot)$ at hand, the
static regret over $T$ slots is then defined as [43]

$\displaystyle\mathcal{R}^{\rm
ST}(T):=\sum_{\tau=1}^{T}\ell_{\tau|\tau-1}-\sum_{\tau=1}^{T}\mathcal{L}({\hat{f}}^{m^{*}}(\mathbf{x}_{\tau});y_{\tau})$
(31)

where $\ell_{\tau|\tau-1}$, defined in (19), captures the ensemble online
Bayesian loss incurred by IE-GP.

(a)

(b)

(c)

Figure 1: Schematic diagrams for (a) IE-GP, (b) SIE-GP, and (c) DIE-GP. The
EGP meta-learner (Owl) assesses the contributions of GP experts (Einsteins)
via the posterior probabilities of time-invariant $i$ or time-varying $i_{t}$.
The GP experts summarize past data in the posterior of
${\mbox{\boldmath$\theta$}}^{m}$ or ${\mbox{\boldmath$\theta$}}^{m}_{t}$.

Although the cumulative online loss in the first sum of (31) has different
form than that of the benchmark, they are comparable by the data likelihood,
where the function is nonrandom. In other words, the online Bayesian loss is
obtained by taking the expectation of the likelihood wrt the online predictive
pdf of the function, thus eliminating the randomness of the function in the
likelihood.

To proceed, we will need the following assumptions.

  * (as1)

The NLL $\mathcal{L}(z_{\tau};y_{\tau})$ is continuously twice differentiable
with $|\frac{d^{2}}{dz_{\tau}^{2}}\mathcal{L}(z_{\tau};y_{\tau})|\leq
c,\forall z_{\tau}$;

  * (as2)

The NLL $\mathcal{L}(z_{t};y_{t})$ is convex and has bounded derivative wrt
$z_{t}$; that is, $|\frac{d}{dz_{t}}\mathcal{L}(z_{t};y_{t})|\leq L$;

  * (as3)

Kernels $\\{\bar{\kappa}^{m}\\}_{s=1}^{M}$ are shift-invariant, standardized
and bounded, that is
$\bar{\kappa}^{m}(\mathbf{x}_{t},\mathbf{x}_{t^{\prime}})\leq
1,\forall\mathbf{x}_{t},\mathbf{x}_{t^{\prime}}$.

Differentiability and convexity of the NLL in (as1)-(as2) are satisfied by
most forms of likelihood in GP-based learning, including the Gaussian
likelihood in GPR, and the logistic one for classification. Conditions in
(as3) hold for a wide class of kernels including Gaussian, Laplace and Cauchy
ones [8]. As the derivations rely on the general form of IE-GP (cf. (15)-(17))
that corresponds to general likelihoods, the regret bound established here
applies to general learning tasks.

To establish the static regret bound of IE-GP, we will need the following
intermediate lemma.

Lemma 1. Under (as1), and with prior of ${\mbox{\boldmath$\theta$}}^{m}$ given
by (14a), the following bound holds concerning the cumulative online Bayesian
loss incurred by the IE-GP and the counterpart from a single RF-based GP
expert with fixed ${\mbox{\boldmath$\theta$}}^{m}_{*}$

$\displaystyle\sum_{\tau=1}^{T}\ell_{\tau|\tau-1}-\sum_{\tau=1}^{T}\mathcal{L}({\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}^{m}_{*};y_{\tau})$
$\displaystyle\leq\frac{\|{\mbox{\boldmath$\theta$}}^{m}_{*}\|^{2}}{2\sigma_{\theta^{m}}^{2}}+n_{\text{RF}}\log\left(1+\frac{T\sigma_{\theta^{m}}^{2}}{2n_{\text{RF}}}\right)+\log
M\;.$ (32)

Proof:  See Sec. 9.1.

Lemma 1 bounds the cumulative online Bayesian loss of IE-GP relative to any
single RF-based GP learner with a fixed strategy. Next, we will work towards
the ultimate static regret by further bounding the loss of RF-based function
estimator relative to the best function estimator in the original RKHS for
each expert.

Theorem 1. Under as(1)-as(3) and with $\hat{f}^{m^{*}}$ belonging to the RHKS
$\mathcal{H}^{m^{*}}$ induced by $\kappa^{m^{*}}$, for a fixed $\epsilon>0$,
the following bound holds with probability at least
$1-2^{8}(\frac{\sigma_{m^{*}}}{\epsilon})^{2}\exp\left(\frac{-n_{\text{RF}}\epsilon^{2}}{4d+8}\right)$

$\displaystyle\sum_{\tau=1}^{T}\ell_{\tau|\tau-1}-\sum_{\tau=1}^{T}\mathcal{L}(\hat{f}^{m^{*}}(\mathbf{x}_{\tau});y_{\tau})\leq$
(33)
$\displaystyle\frac{(1+\epsilon)C^{2}}{2\sigma_{\theta^{m^{*}}}^{2}}+n_{\text{RF}}\log\left(1+\frac{Tc\sigma_{\theta^{m^{*}}}^{2}}{2n_{\text{RF}}}\right)+\log
M+\epsilon LTC$

where $C$ is a constant, and
$\sigma_{m^{*}}^{2}:=\mathbb{E}_{\pi^{m^{*}}_{\bar{\kappa}}}[\|\mathbf{v}^{m^{*}}\|^{2}]$
is the second-order moment of $\mathbf{v}^{m^{*}}$. Setting
$\epsilon=\mathcal{O}(\log T/T)$, the static regret in (31) boils down to

$\displaystyle\mathcal{R}^{\rm ST}(T)=\mathcal{O}(\log T)\;.$ (34)

Proof. See Sec. 9.2.

Theorem 1 asserts that IE-GP incurs no regret on average with cumulative
static regret $\mathcal{O}(\log T)$ over $T$ slots, thereby demonstrating its
robustness to the adversarial setting. It is also worth highlighting that this
regret bound is tighter than that of the deterministic RKHS-based online
multi-kernel counterpart [30] with regret $\mathcal{O}(\sqrt{T})$ in the
static setting.

##  5 EGP for dynamic learning

In the proposed IE-GP, the EGP meta-learner relies on the posterior
probability of a static random variable $i$ to assess contributions of the GP
experts, each of which models the learning function via the time-invariant
parameter vector ${\mbox{\boldmath$\theta$}}^{m}$. Such a stationary setting
implies that IE-GP handles no dynamics in the unknown function. This is also
manifested in Sec. 4 that with batch data in hindsight the optimal function
estimate is associated with one of the GP experts (cf. (29)). To further
enable learning for dynamic functions, the rest of this section will explore
extensions to accommodate time-varying $i_{t}$ and
${\mbox{\boldmath$\theta$}}^{m}_{t}$ for the EGP meta-learner and individual
GP learners, respectively.

###  5.1 Dynamics at EGP meta-learner

Capitalizing on time-dependent $i_{t}\in\mathcal{M}$ to denote the index of
the contributing expert, the EGP meta-learner models the evolution of $i_{t}$
via a Markov chain with prior transition probability $q_{mm^{\prime}}:={\rm
Pr}(i_{t+1}=m|i_{t}=m^{\prime})$ for $m,m^{\prime}\in\mathcal{M}$. Such a
dynamic model allows the learning function to jump among the candidate spaces
associated with the GP experts, yielding the so-termed switching (S) IE-GP
hereafter. The values of $\\{q_{mm^{\prime}}\\}_{m,m^{\prime}}\in[0,1]$ are
user-defined parameters. It is worth pointing out that IE-GP can be regarded
as a special case of SIE-GP with $q_{mm}=1$.

The novel SIE-GP differs from IE-GP in the weight update. To illustrate this,
the per-expert posterior weight in SIE-GP is first adapted as
$w_{t|t}^{m}:={\rm Pr}(i_{t}=m|\mathbf{y}_{t},\mathbf{X}_{t})$ given time-
varying $i_{t}$. Before propagating to $w_{t+1|t+1}^{m}$, the EGP meta-learner
leverages the aforementioned Markov transition model to predict the weight for
GP model $m$ at slot $t+1$ via $w_{t+1|t}^{m}:={\rm
Pr}(i_{t+1}\\!=\\!m|\mathbf{y}_{t},\mathbf{X}_{t})$, which, with
$\\{w_{t|t}^{m}\\}$ available, can be obtained as

$\displaystyle w_{t+1|t}^{m}$
$\displaystyle=\\!\\!\sum_{m^{\prime}=1}^{M}\\!{\rm
Pr}(i_{t+1}=m,i_{t}=m^{\prime}|\mathbf{y}_{t},\mathbf{X}_{t})$
$\displaystyle=\\!\\!\sum_{m^{\prime}=1}^{M}{\rm
Pr}(i_{t+1}\\!=\\!m|i_{t}\\!=\\!m^{\prime},\mathbf{y}_{t},\mathbf{X}_{t}){\rm
Pr}(i_{t}\\!=\\!m^{\prime}|\mathbf{y}_{t},\mathbf{X}_{t})$
$\displaystyle=\\!\\!\sum_{m^{\prime}=1}^{M}q_{m,m^{\prime}}w_{t|t}^{m^{\prime}}$
(35)

where the last equality holds since the evolution of $i_{t}$ is independent
from $\mathbf{y}_{t}$ and $\mathbf{X}_{t}$. With non-zero
$\\{q_{mm^{\prime}}\\}$, the prediction rule in (35) allows for activation of
the previously inactive expert $m$ ($w_{t|t}^{m}\approx 0$), thereby
accommodating switching among candidate GP models or function spaces.

With the predicted weights (35) and the per-expert predictive pdf (15)
available, the EGP meta-learner leverages the sum-product probability rule to
predict the pdf of $y_{t+1}$

$\displaystyle{p}(y_{t+1}|\mathbf{y}_{t},\mathbf{X}_{t+1})=\sum_{m=1}^{M}{p}(y_{t+1},i_{t+1}=m|\mathbf{y}_{t},\\!\mathbf{X}_{t+1})$
$\displaystyle=\sum_{m=1}^{M}{\rm
Pr}(i_{t+1}=m|\mathbf{y}_{t},\mathbf{X}_{t}){p}(y_{t+1}|\mathbf{y}_{t},i_{t+1}\\!=\\!m,\\!\mathbf{X}_{t+1})$
$\displaystyle=\sum_{m=1}^{M}w_{t+1|t}^{m}\
{p}(y_{t+1}|\mathbf{y}_{t},i_{t+1}\\!=\\!m,\mathbf{X}_{t+1})$ (36)

where $w_{t+1|t}^{m}$ replaces $w_{t}^{m}$ in (16) for IE-GP. To facilitate
the upcoming regret analysis, the aggregated online loss for SIE-GP that
accounts for the per-expert loss (18) is defined as

$\displaystyle\ell_{t+1|t}^{\rm
SW}:=-\log\sum_{m=1}^{M}w_{t+1|t}^{m}\exp\left(-l_{t+1|t}^{m}\right)$ (37)

where the superscript “SW”, denoting the switching scenario, is used to
distinguish from the loss from IE-GP in (19).

Upon acquiring $y_{t+1}$, the EGP meta-leaner then updates the per-expert
weight as

$\displaystyle w_{t+1|t+1}^{m}$
$\displaystyle:=\frac{w_{t+1|t}^{m}p(y_{t+1}|\mathbf{y}_{t},i_{t+1}\\!=\\!m,\mathbf{X}_{t+1})}{\sum_{m^{\prime}=1}^{M}w_{t+1|t}^{m^{\prime}}p(y_{t+1}|\mathbf{y}_{t},i_{t+1}\\!=\\!m^{\prime},\mathbf{X}_{t+1})}$
$\displaystyle\ =w_{t+1|t}^{m}\exp(\ell_{t+1|t}^{\rm SW}-l_{t+1|t}^{m})\;.$
(38)

To sum up, relative to IE-GP, the meta-leaner in SIE-GP performs an additional
weight prediction step (35), which contributes to the ensemble predictive pdf
(36) and the weight update (38). Next, the regret analysis of SIE-GP will be
conducted in line with Sec. 4 to account for the adversarial setting.

####  5.1.1 Switching regret analysis

With the underlying assumption that the active model per slot changes over
time, the notion of switching or shifting regret [11, Chapter 5.2] is
leveraged to analyse the performance of SIE-GP in the adversarial setting
where the generative assumptions are violated. Specifically, SIE-GP is
compared with an arbitrary sequence of benchmark functions
$\\{\hat{f}^{i_{\tau}}\in\mathcal{H}^{i_{\tau}}\\}_{\tau=1}^{T}$ with data in
hindsight, yielding the switching regret defined as

$\displaystyle\mathcal{R}^{\rm
SW}(T):=\sum_{\tau=1}^{T}\ell_{\tau|\tau-1}^{\rm
SW}-\underset{i_{1},\ldots,i_{T}}{\min}\sum_{\tau=1}^{T}\mathcal{L}(\hat{f}^{i_{\tau}}(\mathbf{x}_{\tau});y_{\tau})$
(39)

where, as in static regret (31), the loss incurred by the comparator in the
switching case is measured via the NLL.

Aiming at establishing an upper bound for (39), the following two additional
assumptions will be entailed.

  * (as4)

$q_{mm}=q_{0}$, $q_{mm^{\prime}}=\frac{q_{1}}{M-1}$ for
$m,m^{\prime}\in\mathcal{M}$, $q_{0}+q_{1}=1$, and $0\leq
q_{1}<\frac{1}{2}<q_{0}\leq 1$;

  * (as5)

The number of switches for sequence $\\{i_{1},\ldots,i_{T}\\}$ is upper
bounded by $S$, i.e., $\sum_{\tau=1}^{T}I(i_{\tau}\neq i_{\tau+1})\leq S$, and
$S\ll T$.

The transition probabilities dictated by (as4) yield a weight evolution
strategy that is similar to the one given by the fixed-share forecaster in
online expert-based learning [11, Chapter 5.2], where the conditions in (as5)
are also leveraged in the regret analysis to bound the variation of the
benchmark functions. Before obtaining the upper bound for (39) in Theorem 2,
the following two intermediate lemmas will be established first.

Lemma 2. Under as(4)-as(5), the following bound holds true concerning the
cumulative ensemble switching loss and the single expert-based online
counterpart for any sequence $\\{i_{1},\ldots,i_{T}\\}$

$\displaystyle\sum_{\tau=1}^{T}\\!\ell^{\rm
SW}_{\tau|\tau-1}\\!-\\!\\!\sum_{\tau=1}^{T}l^{i_{\tau}}_{\tau|\tau-1}\leq\log
M\\!+\\!S\log T\\!-\\!S\log S\\!+\\!S\;.$ (40)

Proof: See Sec. 9.3.

Lemma 3. Under (as1) and for any sequence $\\{i_{1},\ldots,i_{T}\\}$, the
following bound holds regarding the difference of the cumulative single
expert-based online loss and the counterpart incurred by RF-based benchmark
functions with fixed parameters
$\\{{\mbox{\boldmath$\theta$}}^{m}_{*}\\}_{m=1}^{M}$

$\displaystyle\sum_{\tau=1}^{T}l^{i_{\tau}}_{\tau|\tau-1}-\sum_{\tau=1}^{T}\mathcal{L}({\mbox{\boldmath$\phi$}}^{i_{\tau}\\!\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}^{i_{\tau}}_{*};y_{\tau})$
$\displaystyle\leq\sum_{m=1}^{M}\frac{\|{\mbox{\boldmath$\theta$}}^{m}_{*}\|^{2}}{2\sigma_{\theta^{m}}^{2}}+n_{\text{RF}}M\log\left(1+\frac{Tc\sigma_{\theta^{*}}^{2}}{2n_{\text{RF}}M}\right)$
(41)

where
$\sigma_{\theta^{*}}^{2}:=\max_{m\in\mathcal{M}}\sigma_{\theta^{m}}^{2}$.

Proof: See Sec. 9.4.

Theorem 2. Under as(1)-as(5) and with $\hat{f}^{m}$ belonging to the RHKS
$\mathcal{H}^{m}$ induced by $\kappa^{m}$, for a fixed $\epsilon>0$, the
following bound holds with probability at least
$1-2^{8}(\frac{\sigma_{*}}{\epsilon})^{2}\exp\left(\frac{-n_{\text{RF}}\epsilon^{2}}{4d+8}\right)$

$\displaystyle\sum_{\tau=1}^{T}\ell^{\rm
SW}_{\tau|\tau-1}\\!\\!-\\!\\!\sum_{\tau=1}^{T}\\!\mathcal{L}(\hat{f}^{i_{\tau}}(\mathbf{x}_{\tau});y_{\tau}\\!)\\!\leq\log\\!M\\!+\\!S\log
T\\!\\!-\\!S\log\\!S\\!+\\!S$ $\displaystyle+\epsilon
LTC^{\prime}\\!+\\!\\!\sum_{m=1}^{M}\\!\\!\frac{(1+\epsilon){C^{\prime}}^{2}}{2\sigma_{\theta^{m^{*}}}^{2}}\\!+\\!n_{\text{RF}}M\log\left(1+\frac{Tc\sigma_{\theta^{*}}^{2}}{2n_{\text{RF}}M}\right)$
(42)

where $C^{\prime}$ is some constant, and
$\sigma_{*}^{2}:=\max_{m\in\mathcal{M}}\sigma_{m}^{2}=\max_{m\in\mathcal{M}}\
\mathbb{E}_{\pi^{m^{*}}_{\bar{\kappa}}}[\|\mathbf{v}^{m^{*}}\|^{2}]$. Setting
$\epsilon=\mathcal{O}(\log T/T)$, the switching regret in (39) boils down to

$\displaystyle\mathcal{R}^{\rm SW}(T)=\mathcal{O}(\log T)\;.$ (43)

Proof: See Sec. 9.5.

Even in the presence of model switching, the advocated SIE-GP suffers from
diminishing average regret by explicitly accounting for such switching
dynamics.

###  5.2 Dynamics within each GP expert

The aforementioned SIE-GP accounts for dynamics of expert switching at the EGP
meta-leaner. To further handle a dynamic learning function within each expert
$m$, a time-varying parameter vector ${\mbox{\boldmath$\theta$}}^{m}_{t}$ will
be considered instead of time-invariant ${\mbox{\boldmath$\theta$}}^{m}$ in
IE-GP, yielding the dynamic (D) IE-GP approach. Specifically, DIE-GP captures
dynamics in ${\mbox{\boldmath$\theta$}}^{m}_{t}$ via the random walk model

$\displaystyle{\mbox{\boldmath$\theta$}}^{m}_{t+1}={\mbox{\boldmath$\theta$}}^{m}_{t}+{\mbox{\boldmath$\epsilon$}}_{t+1}^{m}\vspace{-0.5cm}$
(44)

where the noise ${\mbox{\boldmath$\epsilon$}}_{t+1}^{m}$ is white and Gaussian
distributed with mean zero and covariance matrix
$\sigma_{\epsilon^{m}}^{2}\mathbf{I}_{2D}$.

Rather than updating
$p({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})$ as in IE-GP,
expert $m$ in DIE-GP propagates
$p({\mbox{\boldmath$\theta$}}^{m}_{t}|\mathbf{y}_{t},\mathbf{X}_{t})$ across
slots. Taking into account (44), expert $m$ first predicts the pdf of
${\mbox{\boldmath$\theta$}}_{t+1}^{m}$ at the beginning of slot $t+1$ as

$\displaystyle
p(\\!{\mbox{\boldmath$\theta$}}_{t+1}^{m}|\mathbf{y}_{t},\mathbf{X}_{t+1})=\int
p({\mbox{\boldmath$\theta$}}_{t+1}^{m}|{\mbox{\boldmath$\theta$}}_{t}^{m})p({\mbox{\boldmath$\theta$}}_{t}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})d{\mbox{\boldmath$\theta$}}_{t}^{m}\\!\\!$
(45)

which replaces
$p({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}_{t},\mathbf{X}_{t})$ in (15) and
(17) to obtain the predictive pdf
$p(y_{t+1}|\mathbf{y}_{t},\mathbf{X}_{t+1})$, and the posterior
$p({\mbox{\boldmath$\theta$}}^{m}_{t+1}|\mathbf{y}_{t+1},\mathbf{X}_{t+1})$ in
the dynamic setting. Specifically for GPR with per-expert Gaussian posterior
$p({\mbox{\boldmath$\theta$}}^{m}_{t}|\mathbf{y}_{t},\mathbf{X}_{t})=\mathcal{N}({\mbox{\boldmath$\theta$}}_{t}^{m};\hat{{\mbox{\boldmath$\theta$}}}^{m}_{t},\\!{\mbox{\boldmath$\Sigma$}}^{m}_{t})$
, the predictive pdf in (45) is
$p({\mbox{\boldmath$\theta$}}_{t+1}^{m}|\mathbf{y}_{t},\mathbf{X}_{t+1})=\mathcal{N}({\mbox{\boldmath$\theta$}}_{t+1}^{m};\\!\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m},{\mbox{\boldmath$\Sigma$}}_{t}^{m}\\!+\\!\sigma_{\epsilon^{m}}^{2}\mathbf{I}_{2D})$.

Remark 1. The dynamics in EGP meta-learner and GP learners can be readily
combined to yield the DSIE-GP generalization. See Fig. 1 for the schematic
diagrams of (D,S)IE-GP.

Algorithm 2 IE-GPLVM

1:Initialization:

2:for $m=1,\ldots,M$ do

3: Draw vectors
$\\{\mathbf{v}_{j}^{m}\\}_{j=1}^{n_{\text{RF}}}\sim\pi_{\bar{\kappa}^{m}}(\mathbf{v})$;

4: Embed $\mathbf{Y}_{t_{0}}\\!\\!\rightarrow\hat{\mathbf{X}}_{t_{0}}^{m}$ and
obtain hyperparameters (cf. (55));

5:
$\mathbf{B}^{m}_{t_{0}}={\hat{{\mbox{\boldmath$\Phi$}}}_{t_{0}}}^{m\top}\mathbf{Y}_{t_{0}}$
with
${\hat{{\mbox{\boldmath$\Phi$}}}_{t_{0}}}^{m}=[{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\hat{\mathbf{x}}_{1}^{m})\ldots{\mbox{\boldmath$\phi$}}_{\mathbf{v}}(\hat{\mathbf{x}}_{t_{0}}^{m})]^{\top}$;

6: $\mathbf{R}_{t_{0}}^{m}=$
`CholeskyFactor`$({\hat{{\mbox{\boldmath$\Phi$}}}_{t_{0}}}^{m\top}{\hat{{\mbox{\boldmath$\Phi$}}}_{t_{0}}}^{m}+\sigma_{n}^{2}\mathbf{I}_{2n_{\text{RF}}})$;

7:end for

8:for $t=t_{0}+1,t_{0}+2,\ldots$ do

9: Receive datum $\mathbf{y}_{t:}$;

10: for $m=1,\ldots,M$ do

11: Obtain embedding $\hat{\mathbf{x}}_{t}^{m}$ based on (52)

12:
$\mathbf{B}_{t}^{m}=\mathbf{B}_{t-1}^{m}+{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}({\hat{\mathbf{x}}_{t}^{m}})\mathbf{y}^{\top}_{t:}$

13:
$\mathbf{R}_{t}^{m}=\verb|CholeskyUpdate|(\mathbf{R}_{t-1}^{m},{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}({\hat{\mathbf{x}}_{t}^{m}}))$

14: end for

15: Obtain $\hat{\mathbf{x}}_{t}=\hat{\mathbf{x}}_{t}^{m*}$ based on (53);

16: Update $w_{t+1}^{m}$ based on (54);

17:end for

##  6 EGPs for online unsupervised learning

Rather than supervised learning, this section deals with EGP-based latent
variable model (LVM) for unsupervised dimensionality reduction. Consider first
the GPLVM context, where the $D\times 1$ observation
$\mathbf{y}_{\tau:}:=[y_{\tau,1}\ldots y_{\tau,D}]^{\top}$222Notice the
introduction of the colon in the subscript, indicating a multivariate
observation at time $\tau$; not to be confused with $\bm{\mathbf{y}}_{\tau}$
in previous sections.  is linked with the unobserved low-dimensional input
$\mathbf{x}_{\tau}\in\mathbb{R}^{d}$ ($d<D$) via [31]

$\displaystyle y_{\tau j}=f_{j}(\bm{\mathbf{x}}_{\tau})+n_{\tau j},\quad
j=1,\ldots,D$ (46)

where $f_{j}(\cdot)\sim\mathcal{GP}(0,\kappa)$, and $\\{n_{tj}\\}$ are assumed
to be drawn i.i.d. from ${\mathcal{N}}(0,\sigma_{n}^{2})$.

Given $t$ observations
$\bm{\mathbf{Y}}_{t}:=[\bm{\mathbf{y}}_{1:}\ldots\bm{\mathbf{y}}_{t:}]^{\top}\equiv[\bm{\mathbf{y}}_{t}^{(1)}\ldots\bm{\mathbf{y}}_{t}^{(D)}]$,
where $\bm{\mathbf{y}}_{t}^{(j)}:=[y_{1,j}\ldots y_{t,j}]^{\top}$, the
estimate of the low-dimensional embedding
$\bm{\mathbf{X}}_{t}:=[\bm{\mathbf{x}}_{1}\ldots\bm{\mathbf{x}}_{t}]^{\top}$
is sought together with the kernel hyperparameters $\bm{\alpha}$ by solving
[31]

$\displaystyle(\hat{\bm{\mathbf{X}}}_{t},\hat{\bm{\alpha}})=\underset{\bm{\mathbf{X}}_{t},\bm{\alpha}}{\arg\max}\
\log p(\bm{\mathbf{Y}}_{t}|\bm{\mathbf{X}}_{t};\bm{\alpha})+\log
p(\mathbf{X}_{t})$ (47)

where
$p(\bm{\mathbf{Y}}_{t}|\bm{\mathbf{X}}_{t};\bm{\alpha})=\prod_{j=1}^{D}{\mathcal{N}}(\bm{\mathbf{y}}_{t}^{(j)};\bm{\mathbf{0}},\bm{\mathbf{K}}_{t}+\sigma_{n}^{2}\bm{\mathbf{I}}_{t})$
is the so-termed marginal likelihood (ML), and a standard choice for the prior
of $\bm{\mathbf{X}}_{t}$ is
$p(\bm{\mathbf{X}}_{t})=\prod_{\tau=1}^{t}{\mathcal{N}}(\bm{\mathbf{x}}_{\tau};\bm{\mathbf{0}},\sigma_{x}^{2}\bm{\mathbf{I}}_{d})$.

The routine to solve (47) entails inverting the $t\times t$ kernel matrix,
thus incurring unaffordable complexity as in the supervised setting [31]. To
effect scalablility via the RF approximation and accommodate online kernel
adaptation, an IE-GP based LVM is well motivated in accordance with the
preceding discussion.

###  6.1 IE-GPLVM

In the novel IE-GPLVM, the EGP meta-learner employs an ensmeble of GP experts
to independently seek low-dimensional embeddings of the observations. As
before, each expert $m$ will leverage the $\kappa^{m}$-induced RF mapping
${\mbox{\boldmath$\phi$}}_{\bm{\mathbf{v}}}^{m}(\cdot)$ (8) to yield the per-
datum conditional likelihood

$\displaystyle
p(\mathbf{y}_{\tau:}|{\mbox{\boldmath$\Theta$}}^{m},\bm{\mathbf{x}}_{\tau})$
$\displaystyle=\prod_{j=1}^{D}{\mathcal{N}}(y_{\tau
j};\bm{\phi}_{\mathbf{v}}^{m\top}(\bm{\mathbf{x}}_{\tau}){\mbox{\boldmath$\theta$}}_{j}^{m},\sigma_{n}^{2})$
(48)

where the RF-based parameter vectors over $D$ output channels are collected in
${\mbox{\boldmath$\Theta$}}^{m}:=[{\mbox{\boldmath$\theta$}}_{1}^{m},\ldots,{\mbox{\boldmath$\theta$}}_{D}^{m}]^{\top}$,
whose prior pdf is given by

$\displaystyle p({\mbox{\boldmath$\Theta$}}^{m})$
$\displaystyle=\prod_{j=1}^{D}{\mathcal{N}}({\mbox{\boldmath$\theta$}}_{j}^{m};\bm{\mathbf{0}},\sigma_{\theta^{m}}^{2}\bm{\mathbf{I}}_{2n_{\text{RF}}})\;.$
(49)

It is worth mentioning that such an RF-based GPLVM can be regarded as the
nonlinear dual form of probabilistic PCA [31].

To incrementally project subsequent observations to the low-dimensional
embeddings, each expert $m$ in IE-GPLVM relies on the generative model
(48)–(49) to summarize past outputs $\mathbf{Y}_{t}$ and the estimated inputs
$\hat{\mathbf{X}}_{t}^{m}$ in the posterior pdf

$\displaystyle
p({\mbox{\boldmath$\Theta$}}^{m}|\mathbf{Y}_{t},\hat{\mathbf{X}}_{t}^{m})=\prod_{j=1}^{D}\mathcal{N}({\mbox{\boldmath$\theta$}}^{m}_{j};\hat{{\mbox{\boldmath$\theta$}}}^{m}_{t,j},{\mbox{\boldmath$\Sigma$}}_{t}^{m})$
(50)

where the parameter vectors associated with different output channels share
the same covariance matrix ${\mbox{\boldmath$\Sigma$}}_{t}^{m}$. Further, the
per-expert weight assessed by the EGP meta-learner in IE-GPLVM is adapted with
estimated inputs as $w_{t}^{m}:={\rm
Pr}(i=m|\bm{\mathbf{Y}}_{t},\\{\hat{\bm{\mathbf{X}}}_{t}^{\nu}\\}_{\nu=1}^{M}))$.
In the same spirit as IE-GP, each iteration in IE-GPLVM alternates between
estimation of $\bm{\mathbf{x}}_{\tau}$ from $\bm{\mathbf{y}}_{\tau:}$, and
correction of the per-expert weight $w_{t}^{m}$ and posterior pdf (50). Note
that instead of the moments of (50), matrices
$\mathbf{A}_{t}^{m}:=({\mbox{\boldmath$\Sigma$}}_{t}^{m})^{-1}$ and
$\mathbf{B}_{t}^{m}:=({\mbox{\boldmath$\Sigma$}}_{t}^{m})^{-1}\hat{{\mbox{\boldmath$\Theta$}}}_{t}^{m}$
will be equivalently updated with Chelosky decomposition performed via
$\mathbf{A}_{t}^{m}=\mathbf{R}_{t}^{m\top}\mathbf{R}_{t}^{m}$ for numerical
stability (cf. Alg. 2). However, for consistency with previous sections, the
following discussion still uses the moments
$\hat{{\mbox{\boldmath$\Theta$}}}_{t}^{m}$ and
${\mbox{\boldmath$\Sigma$}}_{t}^{m}$.

Estimation. To estimate $\bm{\mathbf{x}}_{t+1}$ based on
$\bm{\mathbf{y}}_{t+1:}$, expert $m$ capitalizes on the parameter posterior
(50) to obtain the conditional likelihood (similar to (22)) as
$p(\bm{\mathbf{y}}_{t+1:}|\bm{\mathbf{Y}}_{t},i\\!=\\!m,\hat{\mathbf{X}}_{t}^{m},\bm{\mathbf{x}}_{t+1})={\mathcal{N}}(\bm{\mathbf{y}}_{t+1:};\hat{\mathbf{y}}_{t+1:}^{m}(\bm{\mathbf{x}}_{t+1}),(\sigma_{t+1}^{m}(\bm{\mathbf{x}}_{t+1}))^{2}\bm{\mathbf{I}}_{D})$,
where the first two moments are functions of $\mathbf{x}_{t+1}$ as

$\displaystyle\hat{\mathbf{y}}_{t+1:}^{m}(\bm{\mathbf{x}}_{t+1})$
$\displaystyle=\hat{{\mbox{\boldmath$\Theta$}}}_{t}^{m}{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{t+1})$
(51a) $\displaystyle(\sigma_{t+1}^{m}(\bm{\mathbf{x}}_{t+1}))^{2}$
$\displaystyle={\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{t+1}){\mbox{\boldmath$\Sigma$}}_{t}^{m}{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{t+1})+\sigma_{n}^{2}\;.$
(51b)

Further imposing a prior on $\bm{\mathbf{x}}_{t+1}$, the maximum-a-posteriori
(MAP) estimate of $\bm{\mathbf{x}}_{t+1}$ is given by

$\displaystyle\hat{\bm{\mathbf{x}}}_{t+1}^{m}\\!=\\!\underset{\bm{\mathbf{x}}_{t+1}}{\arg\max}\
\log
p(\bm{\mathbf{y}}_{t+1:}|\bm{\mathbf{Y}}_{t},\\!i\\!=\\!m,\hat{\mathbf{X}}_{t}^{m}\\!\\!,\bm{\mathbf{x}}_{t+1})\\!+\\!\log
p(\bm{\mathbf{x}}_{t+1}\\!).$ (52)

With the per-expert embeddings
$\\{\hat{\bm{\mathbf{x}}}_{t+1}^{m}\\}_{m=1}^{M}$ at hand, the EGP meta-
learner seeks the final estimate as
$\hat{\bm{\mathbf{x}}}_{t+1}=\hat{\bm{\mathbf{x}}}_{t+1}^{m^{*}}$, where

$\displaystyle m^{*}\\!=\\!\underset{m\in{\cal M}}{\arg\max}\
w_{t}^{m}{p}(\bm{\mathbf{y}}_{t+1:}|\mathbf{Y}_{t},\\!i\\!=\\!m,\hat{\bm{\mathbf{X}}}_{t}^{m}\\!\\!,\hat{\bm{\mathbf{x}}}_{t+1}^{m})\;p(\hat{\bm{\mathbf{x}}}_{t+1}^{m}\\!)\;.$
(53)

It can be readily verified that $(m^{*},\hat{\bm{\mathbf{x}}}_{t+1}^{m^{*}})$
corresponds to the MAP solution of
$p(\mathbf{x}_{t+1},i=m|\mathbf{y}_{t+1:},\mathbf{Y}_{t},\\{\hat{\bm{\mathbf{X}}}_{t}^{\nu}\\}_{\nu=1}^{M})$.333
Observe that
$\displaystyle\max_{\mathbf{x}_{t+1},m}\;p(\mathbf{x}_{t+1},i=m|\mathbf{y}_{t+1:},\mathbf{Y}_{t},\\{\hat{\bm{\mathbf{X}}}_{t}^{\nu}\\}_{\nu=1}^{M})$
$\displaystyle\equiv\max_{\mathbf{x}_{t+1},m}p(\mathbf{y}_{t+1:},\mathbf{x}_{t+1},i\\!=\\!m|\mathbf{Y}_{t},\\{\hat{\bm{\mathbf{X}}}_{t}^{\nu}\\}_{\nu=1}^{M})$
$\displaystyle\equiv\max_{\mathbf{x}_{t+1},m}p(i\\!=\\!m|\mathbf{Y}_{t},\\{\hat{\bm{\mathbf{X}}}_{t}^{\nu}\\}_{\nu=1}^{M}){p}(\mathbf{y}_{t+1:}|\mathbf{Y}_{t},i\\!=\\!m,\hat{\bm{\mathbf{X}}}_{t}^{m}\\!\\!,\mathbf{x}_{t+1})p(\mathbf{x}_{t+1})$

Correction. Upon obtaining the estimate $\hat{\bm{\mathbf{x}}}_{t+1}^{m}$, the
EGP meta-learner updates the per-expert weight

$\displaystyle w_{t+1}^{m}$ $\displaystyle:={\rm
Pr}(i=m|\bm{\mathbf{Y}}_{t+1},\\{\hat{\mathbf{X}}_{t+1}^{\nu}\\}_{\nu=1}^{M})$
$\displaystyle\ \propto
w_{t}^{m}\;{p}(\bm{\mathbf{y}}_{t+1:}|\mathbf{Y}_{t},i=m,\hat{\mathbf{X}}_{t}^{m},\hat{\mathbf{x}}_{t+1}^{m})$
(54)

which is quite intuitive in that it favors experts which assign higher
likelihood to the observation $\bm{\mathbf{y}}_{t+1:}$. In the long run, we
expect the probability mass to concentrate at the expert(s) whose embeddings
best describe the observed data. IE-GPLVM is thus effectively performing
online kernel selection, adapting to the data as they become available.

Meanwhile, the pair
$\\{\hat{\bm{\mathbf{x}}}_{t+1}^{m},\bm{\mathbf{y}}_{t+1:}\\}$ allows expert
$m$ to update the posterior pdf of ${\mbox{\boldmath$\Theta$}}^{m}$ as
$p({\mbox{\boldmath$\Theta$}}^{m}|\mathbf{Y}_{t+1},\hat{\mathbf{X}}_{t+1}^{m})=\prod_{j=1}^{D}\mathcal{N}({\mbox{\boldmath$\theta$}}^{m}_{j};\hat{{\mbox{\boldmath$\theta$}}}^{m}_{t+1,j},{\mbox{\boldmath$\Sigma$}}_{t+1}^{m})$.
As previously mentioned, the update of the moments is implemented by
propagation of $\mathbf{B}_{t}^{m}$ and $\mathbf{R}_{t}^{m}$ as shown in Alg.
2, where CholeskyFactor computes the Cholesky factor (CF) of its argument and
CholeskyUpdate performs a rank-one CF update.

A few remarks are in order.

Remark 2. To obtain the kernel hyperparameters, expert $m$ will leverage the
first $t_{0}$ samples $\mathbf{Y}_{t_{0}}$ to solve the optimisation problem

$\displaystyle(\hat{\bm{\mathbf{X}}}_{t_{0}}^{m},\hat{\bm{\alpha}}^{m})=\underset{\bm{\mathbf{X}}_{t},\bm{\alpha}}{\arg\max}\
\log p(\bm{\mathbf{Y}}_{t_{0}}|\bm{\mathbf{X}}_{t_{0}},i=m;\bm{\alpha})+\log
p(\mathbf{X}_{t_{0}})$ (55)

where the RF-based likelihood for expert $m$ is
$p(\bm{\mathbf{Y}}_{t_{0}}|i\\!=\\!m,\bm{\mathbf{X}}_{t_{0}};\bm{\alpha})=\prod_{j=1}^{D}{\mathcal{N}}(\bm{\mathbf{y}}_{t_{0}}^{(j)};\bm{\mathbf{0}},\sigma_{\theta^{m}}^{2}\bm{\mathbf{\Phi}}_{t_{0}}^{m}\bm{\mathbf{\Phi}}^{m\top}_{t_{0}}+\sigma_{n}^{2}\bm{\mathbf{I}}_{t_{0}})$
with
$\bm{\mathbf{\Phi}}_{t_{0}}^{m}:=[{\mbox{\boldmath$\phi$}}_{\bm{\mathbf{v}}}^{m}(\bm{\mathbf{x}}_{1})\ldots{\mbox{\boldmath$\phi$}}_{\bm{\mathbf{v}}}^{m}(\bm{\mathbf{x}}_{t_{0}})]^{\top}$.
With $\hat{\bm{\mathbf{X}}}_{t_{0}}^{m}$ at hand, expert $m$ further relies on
the generative model (48)–(49) with inputs replaced by the estimates, to
obtain
$p({\mbox{\boldmath$\Theta$}}^{m}|\mathbf{Y}_{t_{0}},\hat{\mathbf{X}}_{t_{0}}^{m})=\prod_{j=1}^{D}\mathcal{N}({\mbox{\boldmath$\theta$}}^{m}_{j};\hat{{\mbox{\boldmath$\theta$}}}^{m}_{{t_{0}},j},{\mbox{\boldmath$\Sigma$}}_{{t_{0}}}^{m})$,
based on which the upcoming data will be processed incrementally (cf. Alg. 2).

Remark 3. In practice, when solving (52), $\bm{\mathbf{x}}_{t+1}$ is
initialized at the embedding corresponding to the point in
$\bm{\mathbf{Y}}_{t}$ that is the nearest neighbor (NN) of
$\bm{\mathbf{y}}_{t+1:}$ [34]. Note, however, that since the per incoming
observation complexity for obtaining the NN, at slot $t$ scales linearly with
$t$, it can significantly surpass the (constant wrt $t$) likelihood evaluation
complexity of ${\mathcal{O}}(n_{\text{RF}}^{2})$. To obtain a scalable
algorithm, we will rely on an approximate kNN search scheme that relies on a
hierarchical graph construct [46]. This approach can be shown to
(approximately, at the limit) achieve ${\mathcal{O}}(\log t)$ complexity; see
[46] for a detailed description.

Remark 4. Dynamic variants of IE-GPLVM can also been pursued along the lines
of Sec. 5.

##  7 Numerical tests

To assess performance, real-data tests are presented here for regression,
classification as well as dimensionality reduction tasks. The supplementary
file also contains synthetic tests that validate the regret bounds in Theorems
1 and 2.

(a)

(b)

(c)

Figure 2: Log scale nMSE plots on (a) Tom’s hardware; (b) SARCOS; and, (c) Air quality datasets. Table I: Statistics of the datasets Task | Datasets | $T$ | $d$ | $D$
---|---|---|---|---
Regression | Tom’s hardware | 9725 | 96 | 1
Regression | SARCOS [1] | 44484 | 21 | 1
Regression | Air Quality [47] | 7322 | 12 | 1
Classification | Banana [48] | 5300 | 2 | 1
Classification | Musk [49] | 6598 | 166 | 1
Classification | Ionosphere [50] | 351 | 34 | 1
GPLVMs | USPS [51] | 5474 | NA | 256
GPLVMs | Oil flow [52] | 1000 | NA | 12
GPLVMs | MNIST | 70000 | NA | 784

(a)

(b)

(c)

Figure 3: Predictive negative log-likelihood on (a) Tom’s hardware; (b)
SARCOS; and, (c) Air quality datasets.

###  7.1 Regression

Regression tests were performed on the SARCOS dataset [1], widely used for
evaluating GP-based approaches, as well as on the Air quality [47], Tom’s
hardware and Twitter datasets [53] from the UCI repository [54]. The
statistics of the datasets are summarized in Table I. We compared the proposed
(D)IE-GP approaches with AdaRaker [30], Incremental Sparse Spectrum Gaussian
Process Regression (I-SSGPR) [9], and the Streaming Sparse Gaussian Process
(SSGP) approach [13], in terms of normalized mean-square error (nMSE) and
running time. SIE-GP is not included for comparison since the datasets exhibit
no switching behavior among the candidate GP models. With $s_{y}^{2}$ denoting
the sample variance of $\mathbf{y}_{T}$, the nMSE is defined as
$\text{nMSE}_{t}:=t^{-1}\sum_{t^{\prime}=1}^{t}(y_{t^{\prime}}-\hat{y}_{t^{\prime}|t^{\prime}-1})^{2}/{s_{y}^{2}}$.

For all RF-based approaches (namely (D)IE-GP, AdaRaker and I-SSGPR) we used
$2n_{\text{RF}}=100$ and the reported results correspond to the run which
resulted in the median nMSE among $101$ runs for the corresponding method.
Finally, all reported runtimes include hyperparameter learning/model
initialization computations performed on the first $1,000$ samples. If for
some expert $m$ and time instance $t$ we have that $w_{t}^{m}=0$, it follows
that $w_{t^{\prime}}^{m}=0$ for all $t^{\prime}>t$ (cf. (20)). Experts with
$w_{t}^{m}<10^{-16}$ were deemed inactive for $t^{\prime}>t$; thus, we set
$w_{t^{\prime}}^{m}=0$ for $t^{\prime}>t$, and avoided unnecessary
prediction/correction steps.

(a)

(b)

(c)

Figure 4: Normalized running times for regression on (a) Tom’s hardware; (b)
SARCOS; and, (c) Air quality datasets.

The kernel dictionary for (D)IE-GP and AdaRaker comprised radial basis
functions (RBFs) with variances from the set $\\{10^{k}\\}_{k=-4}^{6}$. The
automatic relevance determination (ARD) kernel was used for I-SSGPR, as in
[9]. The per kernel noise and prior variances (as well as ARD length scales
for I-SSGPR), were estimated by maximizing the marginal likelihood of the
first $1,000$ samples using the `minimize` function from the GPML toolbox
[55]. The aforementioned samples were not used in the deployment phase. In
DIE-GP, $\sigma_{\epsilon^{m}}^{2}=0.001$ was used for all $m$ and in all
experiments. Regarding SSGP, the ARD kernel was used, the batch size was set
to $300$, the number of inducing points was $100$ and the first $1,000$
samples were used for obtaining an initial model, all as per the original work
[13].

(a)

(b)

(c)

Figure 5: Cumulative classification errors on (a) Banana, (b) Musk, and (c)
Ionosphere datasets.

The nMSE performance of the tested approaches on the Tom’s hardware dataset is
plotted in Fig. 2(a). The proposed (D)IE-GP approaches outperform the
competing alternatives in terms of nMSE while also featuring the lowest
running time, which corresponds to less than $0.3\%$ of that of the most
closely competing (in terms of nMSE) alternative (cf. Fig. 4(a)). The results
on the SARCOS dataset are depicted in Fig. 2(b). Our IE-GP remains competitive
whereas the proposed dynamic variant (DIE-GP) features the lowest nMSE, while
also achieving both faster convergence as well as a runtime that is an order
of magnitude lower than that of the second best (in terms of nMSE) approach
(cf. Fig. 4(b)). These results further highlight the computational efficiency
of the proposed approaches. Similar observations can be made on the Air
quality (cf. Figs. 2,4(c)).

To further demonstrate (D)IE-GP’s uncertainty quantification performance,
tests were conducted among GP-based approaches regarding the predictive
negative log-likelihood (pnLL) as ${\rm pnLL}_{t}:=-\log
p(y_{t}|\mathbf{y}_{t-1},\mathbf{X}_{t})$, which is computable from (16). As
illustrated in Fig. 3, (D)IE-GP always outperform I-SSGPR; while they
outperform SSGP in SARCOS and air-quality datasets; they are comparable in the
Twitter dataset; and perform inferior to SSGP on the Tom’s hardware dataset,
even though SSGP is two orders of magnitude slower than (D)IE-GP.

###  7.2 Classification

Coupled with the logistic likelihood, our IE-GP and the switching variant were
tested for binary classification using Laplace approximation (cf. Sec. 3.2).
The performance of (S)IE-GP were also compared with AdaRaker [30] and SSGP
[13] in terms of classification error and running time on the Banana, Musk,
and Ionosphere datasets whose statistics are provided in Table I. DIE-GP are
not included for comparison since it achieves similar performance relative to
IE-GP. For (S)IE-GP and AdaRaker, the value of $n_{\text{RF}}$ was set to
$15$, and the kernel dictionary is the same as in the regression test.
Regarding SSGP, the number of inducing points was $30$, the batch size was
chosen to be $40$, and the first $20\%$ of the samples were used for model
initialization. For (S)IE-GP, the kernel magnitude $\sigma_{\theta^{m}}^{2}$,
the only hyperparameter per expert, was obtained by maximizing the marginal
likelihood using Laplace apprximation [1, Chapter 5.5.1].

The cumulative classification error and running time of the four competing
approaches are plotted in Figs. 5–6. Clearly, (S)IE-GP outperforms SSGP and
AdaRaker in both classification accuracy and computational efficiency on the
Banana and Musk datasets. Although achieving lower classification error than
(S)IE-GP on the Ionosphere dataset, SSGP runs more than two orders of
magnitude slower. Also, it is worth mentioning that the performance of SSGP
depends on the batch size. Decreasing its value from $40$ to $20$ yields SSGP
performing inferior to (S)IE-GP in classification accuracy. Regarding the two
proposed approaches, SIE-GP achieves comparable classification accuracy
relative to the static IE-GP since the tested datasets exhibit negligible
switching dynamics among the candidate models. In addition, SIE-GP’s higher
running time is explained by the fact that all the GP experts update at all
slots to detect possible model switching, whereas IE-GP implements expert
shutdown for computational efficiency.

(a)

(b)

(c)

Figure 6: Normalized running times for classification on (a) Banana, (b) Musk,
and (c) Ionosphere datasets.

(a)

(b)

(c)

Figure 7: Classification error versus runtime plots for dimensionality
reduction on the (a) MNIST, (b) oil, and (c) USPS datasets. Figure 8:
Visualization of the embedding attained by IE-GPLVM on the USPS dataset.
Colors represent different digits.

###  7.3 Dimensionality reduction

Tests for dimensionality reduction were performed on several benchmark
datasets, including MNIST ($D=784$) as well as the oil flow data ($D=12$)
[52], and the subset of the USPS handwritten digits set ($D=256$) comprising
digits $0-4$. . The latter two datasets were also used in the original GPLVM
paper [31]. Several competing alternatives were considered. GPLVM based
methods comprise the original GPLVM [31, 32], a variational inference based
scheme (varGPLVM) [33], as well as an online GPLVM variant (onGPLVM); see Alg.
2 in [34]. PCA based alternatives encompass online PCA [56, 57], and (batch)
kernel PCA [58]. The embedding dimensionality was set to $d=2$, and the
results presented correspond to the median across 11 trials. Regarding the
proposed IE-GPLVM scheme, $n_{\text{RF}}=50$ random features were used, each
expert relied on a RBF kernel with variance taken from the set
$\\{2^{k}\\}_{k=-3}^{3}$, $t_{0}$ was set to $5\%$ of the number of samples
for MNIST and to $10\%$ for the remaining (smaller) datasets, and (55) was
additionally optimized over $\sigma_{n}^{2}$. For the GPLVM based methods, the
RBF kernel was used, $100$ inducing points were utilized, and the maximum
number of iterations was set to $1,000$. Initializations were provided by
means of PPCA embeddings. Finally, for kernel PCA (kPCA), a grid search was
performed over RBF kernels with variances in $\\{2^{k}\\}_{k=-10}^{10}$ and
the lowest error rate achieved is reported. Note that for kPCA the reported
runtime does not include the computational time required for the grid search.

The error rate of the nearest neighbor classification rule was used as the
performance metric, when applied to the resultant embeddings; see e.g. [31].
The results for the three tested datasets are summarized in Fig. 7, in the
form of error rate versus runtime plots. In the `MNIST` dataset, the proposed
`IE-GPLVM` approach achieves both the lowest overall error rate, as well as
runtime among GPLVM schemes. The only method achieving a somewhat similar
error rate, namely onGPLVM, has a runtime that is more than $20$ times higher
relative to our approach. In the `oil` dataset, our IE-GPLVM achieves similar
error rate (less than $1\%$) to `GPLVM` and `varGPLVM`, while being more than
one and two orders of magnitude faster, respectively. In the `USPS` set, our
approach achieves the lowest overall error rate. The only schemes that achieve
error rates in the same order of magnitude, namely `varGPLVM` and `onGPLVM`,
have runtimes that are $14$ and $21$ times higher, respectively. Finally, in
both experiments, PCA based schemes, although computationally efficient, yield
high error rates. A visualization of the embedding attained is provided for
the proposed IE-GPLVM on the `USPS` dataset (Fig. 8). We can observe that good
separation between clusters of different digits is achieved, in line with the
low classification error rate in Fig. 7.

##  8 Conclusions

This paper put forth an incremental scheme that leverages an ensemble of
scalable RF-based parametric GP learners to jointly infer the unknown function
along with its performance, and a data-driven kernel combination. Dynamic
function learning was enabled through modeling structured dynamics for the EGP
meta-learner and individual GP learners. On the theoretical aspect, regret
analysis was conducted to benchmark even in adversarial settings the novel IE-
GP and its dynamic variant relative to benchmark strategies with data in
hindsight. Further, EGP-based latent variable model is devised for online
kernel-adaptive dimensionality reduction. Extensive experimental results are
provided to illustrate the superior performance of the novel IE-GP schemes.

##  9 Proofs

###  9.1 Proof of Lemma 1

To prove Lemma 1, we will first upper bound the cumulative online Bayesian
loss associated with IE-GP, $\sum_{\tau=1}^{T}\ell_{\tau|\tau-1}$, relative to
that incurred by any RF-based GP expert $m$, namely
$\sum_{\tau=1}^{T}l_{\tau|\tau-1}^{m}$. Reorganizing (20), we have
$\exp(-\ell_{\tau|\tau-1})/\exp(-l_{\tau|\tau-1}^{m})=w_{\tau-1}^{m}/w_{\tau}^{m}$,
multiplying which (20) from $\tau=1$ to $T$, it follows that
$\exp(-\sum_{\tau=1}^{T}\ell_{\tau|\tau-1}+\sum_{\tau=1}^{T}l_{\tau|\tau-1}^{m})=1/(Mw_{T}^{m})$,
yielding

$\displaystyle\sum_{\tau=1}^{T}\ell_{\tau|\tau-1}\\!-\\!\sum_{\tau=1}^{T}l_{\tau|\tau-1}^{m}\\!=\log
M\\!+\\!\log w_{T}^{m}\\!\overset{(a)}{\leq}\\!\log M$ (56)

where $(a)$ holds because $w_{T}^{m}\in[0,1]$.

Next, we will bound the difference between
$\sum_{\tau=1}^{T}l_{\tau|\tau-1}^{m}$ and the cumulative loss incurred by a
fixed strategy ${\mbox{\boldmath$\theta$}}^{m}_{*}$, for any expert
$m\in\mathcal{M}$. For notational brevity, we will drop expert index $m$ in
the remaining of the proof.

Upon defining the cumulative loss over $T$ slots with a time-invariant
$\theta$ as
$\mathcal{L}_{\theta}:=\sum_{\tau=1}^{T}\mathcal{L}({\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{\top}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}};y_{\tau})=-\log
p(\mathbf{y}_{T}|{\mbox{\boldmath$\theta$}},\mathbf{X}_{T})$, the expected
cumulative loss over any pdf $q({\mbox{\boldmath$\theta$}})$ is [43]

$\displaystyle\bar{\mathcal{L}}_{q_{\theta}}:=\mathbb{E}_{q}[\mathcal{L}_{\theta}]=\int_{{\mbox{\boldmath$\theta$}}}q({\mbox{\boldmath$\theta$}})\mathcal{L}_{\theta}d{\mbox{\boldmath$\theta$}}\;.$

On the other hand, the following equality holds for the cumulative online
Bayesian loss based on Bayes’ rule

$\displaystyle\sum_{\tau=1}^{T}\ell_{\tau|\tau-1}$
$\displaystyle=\sum_{\tau=1}^{T}-\log
p(y_{\tau}|\mathbf{y}_{\tau-1},\mathbf{X}_{\tau})=-\log
p(\mathbf{y}_{T}|\mathbf{X}_{T})\;.$

Let
$q({\mbox{\boldmath$\theta$}})=\mathcal{N}({\mbox{\boldmath$\theta$}};{\mbox{\boldmath$\theta$}}_{*},\xi^{2}\mathbf{I}_{2n_{\text{RF}}})$
with the variational parameter $\xi$ to be tuned later, and
$p({\mbox{\boldmath$\theta$}})=\mathcal{N}({\mbox{\boldmath$\theta$}};\mathbf{0},\sigma_{\theta}^{2}\mathbf{I}_{2n_{\text{RF}}})$.
It then follows that

$\displaystyle\sum_{\tau=1}^{T}l_{\tau|\tau-1}-\bar{\mathcal{L}}_{q_{\theta}}=\int
q({\mbox{\boldmath$\theta$}})\log\frac{p(\mathbf{y}_{T}|{\mbox{\boldmath$\theta$}},\mathbf{X}_{T})}{p(\mathbf{y}_{T}|\mathbf{X}_{T})}d{\mbox{\boldmath$\theta$}}$
$\displaystyle\overset{(a)}{=}\int
q({\mbox{\boldmath$\theta$}})\log\frac{p({\mbox{\boldmath$\theta$}}|\mathbf{y}_{T},\mathbf{X}_{T})}{p({\mbox{\boldmath$\theta$}})}d{\mbox{\boldmath$\theta$}}$
$\displaystyle=\int
q({\mbox{\boldmath$\theta$}})\log\frac{q({\mbox{\boldmath$\theta$}})}{p({\mbox{\boldmath$\theta$}})}d{\mbox{\boldmath$\theta$}}-\int
q({\mbox{\boldmath$\theta$}})\log\frac{q({\mbox{\boldmath$\theta$}})}{p({\mbox{\boldmath$\theta$}}|\mathbf{y}_{T},\mathbf{X}_{T})}d{\mbox{\boldmath$\theta$}}\
\ $ $\displaystyle=\\!{\rm
KL}(q(\\!{\mbox{\boldmath$\theta$}})\|p(\\!{\mbox{\boldmath$\theta$}})\\!)\\!-\\!{\rm
KL}(q(\\!{\mbox{\boldmath$\theta$}})\|p(\\!{\mbox{\boldmath$\theta$}}|\mathbf{y}_{T},\mathbf{X}_{T})\\!)\\!\overset{(b)}{\leq}\\!{\rm
KL}(q(\\!{\mbox{\boldmath$\theta$}})\|p(\\!{\mbox{\boldmath$\theta$}})\\!)$
$\displaystyle=2n_{\text{RF}}\log\sigma_{\theta}+\frac{\|{\mbox{\boldmath$\theta$}}_{*}\|^{2}+2n_{\text{RF}}\xi^{2}}{2\sigma_{\theta}^{2}}\\!-\\!n_{\text{RF}}\\!-\\!2n_{\text{RF}}\log\xi$
(57)

where $(a)$ holds since Bayes’ rule yields
$p(\mathbf{y}_{T}|{\mbox{\boldmath$\theta$}},\mathbf{X}_{T})p({\mbox{\boldmath$\theta$}})=p(\mathbf{y}_{T}|\mathbf{X}_{T})p({\mbox{\boldmath$\theta$}}|\mathbf{y}_{T},\mathbf{X}_{T})$;
and, $(b)$ comes from the fact that ${\rm
KL}(q(\\!{\mbox{\boldmath$\theta$}})\|p(\\!{\mbox{\boldmath$\theta$}}|\mathbf{y}_{T},\\!\mathbf{X}_{T})\\!)\geq
0$.

The last step towards bounding
$\sum_{\tau=1}^{T}l_{\tau|\tau-1}-\mathcal{L}_{\theta_{*}}$ is to establish an
upper bound for $\bar{\mathcal{L}}_{q_{\theta}}-\mathcal{L}_{\theta_{*}}$. To
this end, let
$z_{\tau}={\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{\top}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}$
and
$z_{\tau}^{*}={\mbox{\boldmath$\phi$}}^{\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}_{*}$.
Taking the Taylor’s expansion of $\mathcal{L}(z_{\tau};y_{\tau})$ around
$z_{\tau}^{*}$, yields

$\displaystyle\mathcal{L}(z_{\tau};y_{\tau})$
$\displaystyle=\mathcal{L}(z_{\tau}^{*};y_{\tau})+\frac{d\mathcal{L}(z_{\tau}^{*};y_{\tau})}{dz_{\tau}}(z_{\tau}-z_{\tau}^{*})$
$\displaystyle\ \ \ \
+\frac{d^{2}}{dz_{\tau}^{2}}\mathcal{L}(h(z_{\tau});y_{\tau})\frac{(z_{\tau}-z_{\tau}^{*})^{2}}{2}$
(58)

where $h(z_{\tau})$ is some function lying between $z_{\tau}$ and
$z_{\tau}^{*}$. Taking the expectation of (58) wrt
$q({\mbox{\boldmath$\theta$}})$, leads to

$\displaystyle\mathbb{E}_{q}[\mathcal{L}(z_{\tau};y_{\tau})]\\!-\\!\mathcal{L}(z_{\tau}^{*};y_{\tau})$
$\displaystyle=\mathbb{E}_{q}\left[\frac{d^{2}}{dz_{\tau}^{2}}\mathcal{L}(h(z_{\tau});y_{\tau})\frac{(z_{\tau}-z_{\tau}^{*})^{2}}{2}\right]$
$\displaystyle\overset{(a)}{\leq}c\mathbb{E}\left[\frac{(z_{\tau}-z_{\tau}^{*})^{2}}{2}\right]\overset{(b)}{\leq}\frac{c\xi^{2}}{2}$
(59)

where $(a)$ makes use of (as1) that
$\left|\frac{d^{2}}{dz^{2}}\mathcal{L}(z;y)\right|\leq c\ \forall z$, and
$(b)$ relies on the equality
$\|{\mbox{\boldmath$\phi$}}_{\mathbf{v}}(\mathbf{x}_{\tau})\|^{2}=1$.

Summing (59) from $\tau=1$ to $T$, we have

$\displaystyle\bar{\mathcal{L}}_{q_{\theta}}\leq\mathcal{L}_{\theta_{*}}+\frac{Tc\xi^{2}}{2}$
(60)

which, in conjunction with (57), yields the inequality

$\displaystyle\sum_{\tau=1}^{T}l_{\tau|\tau-1}-\mathcal{L}_{\theta_{*}}$ (61)
$\displaystyle\leq\frac{Tc\xi^{2}}{2}\\!+\\!2n_{\text{RF}}\log\sigma_{\theta}\\!+\\!\frac{\|{\mbox{\boldmath$\theta$}}_{*}\|^{2}+2n_{\text{RF}}\xi^{2}}{2\sigma_{\theta}^{2}}\\!-\\!n_{\text{RF}}\\!-\\!2n_{\text{RF}}\log\xi\;.$

Replacing the RHS of (61), a convex function of $\xi$, with the minimal value
taken at
$\xi^{2}=\frac{2n_{\text{RF}}\sigma_{\theta}^{2}}{2n_{\text{RF}}+Tc\sigma_{\theta}^{2}}$,
simplifies (61) for any expert $m\in\mathcal{M}$ to

$\displaystyle\sum_{\tau=1}^{T}l_{\tau}^{m}\\!-\\!\\!\sum_{\tau=1}^{T}\mathcal{L}\\!\left(\\!{\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}\\!(\mathbf{x}_{\tau}\\!){\mbox{\boldmath$\theta$}}^{m}_{*};y_{\tau}\\!\right)\\!\leq\\!\frac{\|{\mbox{\boldmath$\theta$}}^{m}_{*}\|^{2}}{2\sigma_{\theta^{m}}^{2}}\\!+\\!n_{\text{RF}}\log\\!\left(\\!1\\!+\\!\frac{Tc\sigma_{\theta^{m}}^{2}}{2n_{\text{RF}}}\\!\right)$

which, together with (56), readily prove Lemma 1.

###  9.2 Proof of Theorem 1

For a given shift-invariant standardized kernel $\bar{\kappa}^{m}$, the
maximum point-wise error of the RF kernel approximant is uniformly bounded
with probability at least
$1-2^{8}(\frac{\sigma_{m}}{\epsilon})^{2}\exp\left(\frac{-n_{\text{RF}}\epsilon^{2}}{4d+8}\right)$
[8]

$\displaystyle\sup_{\mathbf{x}_{i},\mathbf{x}_{j}\in\mathcal{X}}\left|{\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{i}){\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{j})-\bar{\kappa}^{m}(\mathbf{x}_{i},\mathbf{x}_{j})\right|<\epsilon$
(62)

where $\epsilon$ is a given constant, $n_{\text{RF}}$ is the number of
spectral feature vectors, $d$ is the dimension of $\mathbf{x}$, and
$\sigma_{m}^{2}:=\mathbb{E}_{\pi_{\bar{\kappa}}^{m}}[\|\mathbf{v}^{m}\|^{2}]$
is the second-order moment of the RF vector $\mathbf{v}^{m}$.

The optimal function estimator in $\mathcal{H}^{m}$ incurred by $\kappa^{m}$
is

$\displaystyle\hat{f}^{m}(\mathbf{x}):=\sum_{\tau=1}^{T}\hat{\alpha}^{m}_{\tau}\kappa^{m}(\mathbf{x},\mathbf{x}_{\tau})=\sigma_{\theta^{m}}^{2}\sum_{\tau=1}^{T}\hat{\alpha}^{m}_{\tau}\bar{\kappa}^{m}(\mathbf{x},\mathbf{x}_{\tau})$
(63)

and its RF-based approximant is
$\check{f}^{m}_{*}(\mathbf{x}):={\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}){\mbox{\boldmath$\theta$}}^{m}_{*}$
with
${\mbox{\boldmath$\theta$}}^{m}_{*}:=\sigma_{\theta^{m}}^{2}\sum_{\tau=1}^{T}\hat{\alpha}_{\tau}^{m}{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{\tau})$.
We then have that

$\displaystyle\left|\sum_{\tau=1}^{T}\mathcal{L}\left(\check{f}^{m}_{*}(\mathbf{x}_{\tau});y_{\tau}\right)-\sum_{\tau=1}^{T}\mathcal{L}\left(\hat{f}^{m}(\mathbf{x}_{\tau});y_{\tau}\right)\right|$
(64)
$\displaystyle\overset{(a)}{\leq}\sum_{\tau=1}^{T}\left|\mathcal{L}\left(\check{f}^{m}_{*}(\mathbf{x}_{\tau});y_{\tau}\right)-\mathcal{L}\left(\hat{f}^{m}(\mathbf{x}_{\tau});y_{\tau}\right)\right|$
$\displaystyle\overset{(b)}{\leq}\sum_{\tau=1}^{T}L\sigma_{\theta^{m}}^{2}\\!\left|\sum_{\tau^{\prime}=1}^{T}\\!\hat{\alpha}_{\tau^{\prime}}^{m}{\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{m\top}\\!(\mathbf{x}_{\tau}){\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{\tau^{\prime}})\\!-\\!\\!\\!\sum_{\tau^{\prime}=1}^{T}\\!\\!\hat{\alpha}_{\tau^{\prime}}^{m}\bar{\kappa}^{m}(\mathbf{x}_{\tau},\mathbf{x}_{\tau^{\prime}})\right|$
$\displaystyle\overset{(c)}{\leq}\sum_{\tau=1}^{T}L\sigma_{\theta^{m}}^{2}\sum_{\tau^{\prime}=1}^{T}\left|\hat{\alpha}_{\tau^{\prime}}^{m}\right|\left|{\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{m\top}(\mathbf{x}_{\tau}){\mbox{\boldmath$\phi$}}_{\mathbf{v}}^{m}(\mathbf{x}_{\tau^{\prime}})-\bar{\kappa}^{m}(\mathbf{x}_{\tau},\mathbf{x}_{\tau^{\prime}})\right|$

where $(a)$ follows from the triangle inequality; $(b)$ makes use of (as2),
which states the convexity and bounded derivative of $\mathcal{L}(z;y)$ wrt
$z$, and $(c)$ results from the Cauchy-Schwarz inequality. Leveraging (62) to
upper bound the RHS of (64), we find

$\displaystyle\left|\sum_{\tau=1}^{T}\mathcal{L}\left(\check{f}^{m}_{*}(\mathbf{x}_{\tau});y_{\tau}\right)-\sum_{\tau=1}^{T}\mathcal{L}\left(\hat{f}^{m}(\mathbf{x}_{\tau});y_{\tau}\right)\right|$
$\displaystyle\leq\sum_{\tau=1}^{T}L\sigma_{\theta^{m}}^{2}\epsilon\sum_{\tau=1}^{T}\left|\hat{\alpha}_{\tau}^{m}\right|\leq\epsilon
LTC,\ {\rm w.h.p.}$ (65)

where
$C:=\underset{m\in\mathcal{M}}{\max}\sum_{\tau=1}^{T}\sigma_{\theta^{m}}^{2}|\hat{\alpha}_{\tau}^{m}|$.
It thus holds w.h.p. that

$\displaystyle\sum_{\tau=1}^{T}\\!\mathcal{L}\\!\left(\\!{\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){{\mbox{\boldmath$\theta$}}^{m}_{*}};y_{\tau}\\!\right)\\!-\\!\sum_{\tau=1}^{T}\\!\mathcal{L}\\!\left(\\!\hat{f}^{m}(\mathbf{x}_{\tau});y_{\tau}\\!\right)\leq\epsilon
LTC.$ (66)

On the other hand, the uniform convergence bound in (62) and (as3) imply
w.h.p. that

$\displaystyle\underset{\mathbf{x}_{\tau},\mathbf{x}_{\tau^{\prime}}\in\mathcal{X}}{\sup}\\!\\!{\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}\\!\\!(\mathbf{x}_{\tau}){\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{\tau^{\prime}})\\!\leq\\!1\\!+\\!\epsilon$
(67)

based on which

$\displaystyle\|{\mbox{\boldmath$\theta$}}^{m}_{*}\|^{2}:=\left\|\sigma_{\theta^{m}}^{2}\sum_{\tau=1}^{T}\hat{\alpha}^{m}_{\tau}{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{\tau})\right\|^{2}$
(68)
$\displaystyle=\sigma_{\theta^{m}}^{4}\sum_{\tau=1}^{T}\sum_{\tau^{\prime}=1}^{T}\hat{\alpha}^{m}_{\tau}\hat{\alpha}^{m}_{\tau^{\prime}}{\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{\tau^{\prime}})\leq(1+\epsilon)C^{2}\;.$

Hence, in conjunction with (68), (66) and Lemma 1, it follows for any
$m\in\mathcal{M}$ that

$\displaystyle\sum_{\tau=1}^{T}\ell_{\tau|\tau-1}-\sum_{\tau=1}^{T}\mathcal{L}(\hat{f}^{m}(\mathbf{x}_{\tau});y_{\tau})$
(69) $\displaystyle\leq
n_{\text{RF}}\log\left(1+\frac{Tc\sigma_{\theta^{m}}^{2}}{2n_{\text{RF}}}\right)+\log
M+\frac{(1+\epsilon)C^{2}}{2\sigma_{\theta^{m}}^{2}}+\epsilon LTC$

thus completing the proof of Theorem 1 with $m=m^{*}$.

###  9.3 Proof of Lemma 2

Inspired by [11], the proof of Lemma 2 will be conducted relying on the notion
of compound experts, each associated with a sequence of contributing models
over $T$ slots, denoted as $\mathbf{i}_{T}=[i_{1},\ldots,i_{T}]^{\top}$. Let
$\bar{w}_{t}(\mathbf{i}_{T})$ represent the posterior weight of compound
expert $\mathbf{i}_{T}$ at slot $t$ as $\bar{w}_{t}(\mathbf{i}_{T}):={\rm
Pr}(\mathbf{i}_{T}|\mathbf{y}_{t},\mathbf{X}_{t})$, which is updated at slot
$t+1$ as

$\displaystyle\bar{w}_{t+1}(\mathbf{i}_{T})$
$\displaystyle=\frac{\bar{w}_{t}(\mathbf{i}_{T})p(y_{t+1}|\mathbf{y}_{t},\mathbf{i}_{T},\mathbf{X}_{t+1})}{p(y_{t+1}|\mathbf{y}_{t},\mathbf{X}_{t+1})}$
$\displaystyle\overset{(a)}{=}\bar{w}_{t}(\mathbf{i}_{T})\exp\left(\bar{\ell}_{t+1|t}-l_{t+1|t}^{i_{t+1}}\right)$
(70)

where equality $(a)$ holds since
$p(y_{t+1}|\mathbf{y}_{t},\mathbf{i}_{T},\mathbf{X}_{t+1})=p(y_{t+1}|\mathbf{y}_{t},i_{t+1},\mathbf{X}_{t+1})$,
and $\bar{\ell}_{t+1|t}$ signifies the ensemble online loss accounting for all
the compound experts as

$\displaystyle\bar{\ell}_{t+1|t}$
$\displaystyle:=-\log\sum_{i_{1},\ldots,i_{T}}\bar{w}_{t}(\mathbf{i}_{T})p(y_{t+1}|\mathbf{y}_{t},\mathbf{i}_{T},\mathbf{X}_{t+1})$
$\displaystyle\
=-\log\sum_{i_{t+1}}\bar{w}_{t+1}^{i_{t+1}}\exp(-l_{t+1|t}^{i_{t+1}})$ (71)

with marginal weight at slot $t+1$ given by

$\displaystyle\bar{w}_{t+1}^{i_{t+1}}:=\sum_{i_{1},\ldots,i_{t},i_{t+2},\ldots,i_{T}}\bar{w}_{t}(\mathbf{i}_{T})\;.$
(72)

Next, the equivalence of $\bar{\ell}_{t+1|t}$ and ${\ell}_{t+1|t}^{\rm SW}$
(37) will be established via proving the equality
$\bar{w}_{t+1}^{i_{t+1}}=w_{t+1|t}^{i_{t+1}}$ by induction. Towards this, we
will first leverage (70) to obtain

$\displaystyle\bar{w}_{t}(\mathbf{i}_{T})\propto\bar{w}_{0}(\mathbf{i}_{T})\exp\left(\sum_{\tau=1}^{t}-l_{\tau|\tau-1}^{i_{\tau}}\right)$

based on which, (72) can be rewritten as

$\displaystyle\bar{w}_{t+1}^{i_{t+1}}$
$\displaystyle\propto\sum_{i_{1},\ldots,i_{t}}\bar{w}_{0}(\mathbf{i}_{t+1})\exp\left(\sum_{\tau=1}^{t}-l_{\tau|\tau-1}^{i_{\tau}}\right)$
$\displaystyle=\sum_{i_{t}}\frac{\bar{w}_{0}(\mathbf{i}_{t+1})}{\bar{w}_{0}(\mathbf{i}_{t})}e^{-l_{t|t-1}^{i_{t}}}\sum_{i_{1},\ldots,i_{t-1}}\bar{w}_{0}(\mathbf{i}_{t})\exp\left(\sum_{\tau=1}^{t-1}-l_{\tau|\tau-1}^{i_{\tau}}\right)$
$\displaystyle\propto\sum_{i_{t}}{\rm
Pr}(i_{t+1}|i_{t})\exp\left(-l_{t|t-1}^{i_{t}}\right)\bar{w}_{t}^{i_{t}}$
$\displaystyle\overset{(a)}{=}\sum_{i_{t}}{\rm
Pr}(i_{t+1}|i_{t})\exp\left(-l_{t|t-1}^{i_{t}}\right){w}_{t|t-1}^{i_{t}}$
$\displaystyle\overset{(b)}{\propto}\sum_{i_{t}}{\rm
Pr}(i_{t+1}|i_{t})w_{t|t}^{i_{t}}\overset{(c)}{=}w_{t+1|t}^{i_{t+1}}$ (73)

where $(a)$ follows due to the induction assumption that
$\bar{w}_{t}^{i_{t}}=w_{t|t-1}^{i_{t}}$; and, $(b)$ and $(c)$ are based on
(38) and (35), respectively.

With $\bar{\ell}_{t+1|t}={\ell}_{t+1|t}^{\rm SW}$ being established according
to (73), multiplying (70) from $t=1$ to $T$ yields

$\displaystyle\sum_{\tau=1}^{T}\ell_{\tau+1|\tau}^{\rm
SW}-\sum_{\tau=1}^{T}l_{\tau+1|\tau}^{i_{\tau+1}}$
$\displaystyle=\log\bar{w}_{T}(\mathbf{i}_{T})-\log\bar{w}_{0}(\mathbf{i}_{T})$
$\displaystyle\leq-\log\bar{w}_{0}(\mathbf{i}_{T})\;.$ (74)

Since $\bar{w}_{0}(\mathbf{i}_{T})\\!=\\!{\rm Pr}(i_{1})\prod_{\tau=2}^{T}{\rm
Pr}(i_{\tau}|i_{\tau-1})\\!=\\!\frac{1}{M}q_{0}^{T-s}q_{1}^{s}$ with $s$
denoting the number of switches in $\mathbf{i}_{T}$, (74) can be rewritten as

$\displaystyle\sum_{\tau=1}^{T}\ell_{\tau+1|\tau}^{\rm
SW}\\!\\!-\\!\\!\sum_{\tau=1}^{T}l_{\tau+1|\tau}^{i_{\tau+1}}$
$\displaystyle\leq\log\\!M\\!-\\!(T\\!-\\!s)\log\\!q_{0}\\!-\\!s\log q_{1}$
$\displaystyle\overset{(b)}{\leq}\log\\!M\\!-\\!T\log\\!q_{0}\\!+\\!S\\!\log\frac{q_{0}}{1-q_{0}}$
(75)

where the inequality in $(a)$ results from $s\leq S$ and $q_{0}\geq
q_{1}=1-q_{0}$ based on (as4)-(as5). As the RHS of (75) is a convex function
of $q_{0}$, the following holds true upon setting it to its minimal value
taken at $q_{0}^{*}=(T-S)/T$

$\displaystyle\sum_{\tau=1}^{T}\ell_{\tau+1|\tau}^{\rm
SW}\\!-\\!\sum_{\tau=1}^{T}l_{\tau+1|\tau}^{i_{\tau+1}}$
$\displaystyle\leq\log
M\\!-\\!S\log\frac{S}{T}\\!+\\!(T\\!-\\!S)\log\frac{T}{T-S}$

which, upon leveraging $(T-S)\log\frac{T}{T-S}\leq(T-S)\frac{S}{T-S}=S$ for
$S\ll T$, yields Lemma 2.

###  9.4 Proof of Lemma 3

For any sequence $\mathbf{i}_{T}$, the cumulative loss over $T$ slots measured
by negative log-likelihood for fixed parameter set
${\mbox{\boldmath$\Theta$}}:=\\{{\mbox{\boldmath$\Theta$}}^{m}\\}_{m=1}^{M}$
is given by $\mathcal{L}_{\Theta}^{\mathbf{i}_{T}}:=-\log
p(\mathbf{y}_{T}|{\mbox{\boldmath$\Theta$}},\mathbf{i}_{T},\mathbf{X}_{T})=\sum_{\tau=1}^{T}\mathcal{L}({\mbox{\boldmath$\phi$}}^{i_{\tau}\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}^{i_{\tau}};y_{\tau})$,
whose expected value over factorized pdf
$q({\mbox{\boldmath$\Theta$}})=\prod_{m=1}^{M}q({\mbox{\boldmath$\theta$}}^{m})=\prod_{m=1}^{M}\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};{\mbox{\boldmath$\theta$}}^{m}_{*},\xi_{m}^{2}\mathbf{I}_{2D},)$
is

$\displaystyle\bar{\mathcal{L}}_{q_{\Theta}}^{\mathbf{i}_{T}}$
$\displaystyle:=\\!\\!\int\\!\\!\mathcal{L}_{\Theta}^{\mathbf{i}_{T}}q({\mbox{\boldmath$\Theta$}})d{\mbox{\boldmath$\Theta$}}\\!=\\!\sum_{\tau=1}^{T}\\!\int\\!\\!\mathcal{L}({\mbox{\boldmath$\phi$}}^{i_{\tau}\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}^{i_{\tau}};y_{\tau})q({\mbox{\boldmath$\theta$}}^{i_{\tau}})d{\mbox{\boldmath$\theta$}}^{i_{\tau}}$

which can be re-expressed as
$\bar{\mathcal{L}}_{q_{\Theta}}^{\mathbf{i}_{T}}=\sum_{m=1}^{M}\bar{\mathcal{L}}_{q_{\Theta}}^{m}$,
where

$\displaystyle\bar{\mathcal{L}}_{q_{\Theta}}^{m}:=\sum_{\tau\in\mathcal{T}_{m}}\int\mathcal{L}({\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}^{m};y_{\tau})q({\mbox{\boldmath$\theta$}}^{m})d{\mbox{\boldmath$\theta$}}^{m}$

with $\mathcal{T}_{m}$ collecting the $T_{m}$ slot indices when GP model from
expert $m$ is in action, that is,
$\mathcal{T}_{m}:=\\{\tau|i_{\tau}\\!=\\!m,1\\!\leq\\!\tau\\!\leq\\!T\\}$.

On the other hand, the cumulative online loss for any $\mathbf{i}_{T}$ is

$\displaystyle\sum_{\tau=1}^{T}l^{i_{\tau}}_{\tau|\tau-1}\\!\\!=\\!\sum_{\tau=1}^{T}\\!-\\!\log
p(y_{\tau}|\mathbf{y}_{\tau-1},i_{\tau},\\!\mathbf{X}_{\tau})\\!=\\!-\log
p(\mathbf{y}_{T}|\mathbf{i}_{T},\\!\mathbf{X}_{T}).$

With
$p({\mbox{\boldmath$\Theta$}})\\!=\\!\prod_{m=1}^{M}p({\mbox{\boldmath$\theta$}}^{m})\\!=\\!\prod_{m=1}^{M}\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\mathbf{0}_{2n_{\text{RF}}},\sigma_{\theta^{m}}^{2}\mathbf{I}_{2n_{\text{RF}}})$,
the following inequality can be proved in accordance with (57)

$\displaystyle\sum_{\tau=1}^{T}l_{\tau|\tau-1}^{i_{\tau}}-\bar{\mathcal{L}}_{q_{\Theta}}^{\mathbf{i}_{T}}=\int
q({\mbox{\boldmath$\Theta$}})\log\frac{p(\mathbf{y}_{T}|{\mbox{\boldmath$\Theta$}},\mathbf{i}_{T},\mathbf{X}_{T})}{p(\mathbf{y}_{T}|\mathbf{i}_{T},\mathbf{X}_{T})}d{\mbox{\boldmath$\Theta$}}$
$\displaystyle=\int
q({\mbox{\boldmath$\Theta$}})\log\frac{p({\mbox{\boldmath$\Theta$}}|\mathbf{y}_{T},\mathbf{i}_{T},\mathbf{X}_{T})}{p({\mbox{\boldmath$\Theta$}})}d{\mbox{\boldmath$\Theta$}}$
$\displaystyle=\sum_{m=1}^{M}\int
q({\mbox{\boldmath$\theta$}}^{m})\left(\log\frac{q({\mbox{\boldmath$\theta$}}^{m})}{p({\mbox{\boldmath$\theta$}}^{m})}-\log\frac{q({\mbox{\boldmath$\theta$}}^{m})}{p({\mbox{\boldmath$\theta$}}^{m}|\mathbf{y}^{m},\mathbf{X}^{m})}\right)d{\mbox{\boldmath$\theta$}}^{m}$
$\displaystyle\leq\sum_{m=1}^{M}KL(q({\mbox{\boldmath$\theta$}}^{m})\|p({\mbox{\boldmath$\theta$}}^{m}))$
$\displaystyle=\sum_{m=1}^{M}\\!\\!\left(\\!2n_{\text{RF}}\\!\log\frac{\sigma_{\theta^{m}}}{\xi_{m}}\\!+\\!\frac{1}{2\sigma_{\theta}^{2}}\left(\|{\mbox{\boldmath$\theta$}}_{*}^{m}\|^{2}\\!+\\!2n_{\text{RF}}\xi_{m}^{2}\right)\\!-\\!n_{\text{RF}}\\!\right).$
(76)

Next, following the derivations in (58)–(60) yields for any $m$

$\displaystyle\bar{\mathcal{L}}_{q_{\Theta}}^{m}\leq\sum_{t\in\mathcal{T}_{m}}\mathcal{L}({\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{t}){\mbox{\boldmath$\theta$}}^{m}_{*};y_{t})+\frac{T_{m}c\xi_{m}^{2}}{2}$

the sum of which over $m$ in conjunction with (76) results in

$\displaystyle\sum_{\tau=1}^{T}l_{\tau|\tau-1}^{i_{\tau}}-\sum_{\tau=1}^{T}\mathcal{L}({\mbox{\boldmath$\phi$}}^{i_{\tau}\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}^{i_{\tau}}_{*};y_{\tau})\leq$
(77)
$\displaystyle\sum_{m=1}^{M}\\!\\!\\!\left(\\!2n_{\text{RF}}\\!\log\frac{\sigma_{\theta^{m}}}{\xi_{m}}\\!+\\!\frac{1}{2\sigma_{\theta}^{2}}\left(\|{\mbox{\boldmath$\theta$}}_{*}^{m}\|^{2}\\!+\\!2n_{\text{RF}}\xi_{m}^{2}\right)\\!-\\!n_{\text{RF}}+\\!\frac{T_{m}c\xi_{m}^{2}}{2}\\!\right).$

It is evident that the RHS of (77) is a sum of $M$ convex functions of
$\xi_{m}$, each of which takes minimal value at
$\xi_{m}^{*}=\sqrt{(2n_{\text{RF}}\sigma_{\theta^{m}}^{2})/(2n_{\text{RF}}+T_{m}c\sigma_{\theta^{m}}^{2})}$.
Setting the RHS of (77) to its minimal value yields

$\displaystyle\sum_{\tau=1}^{T}l^{i_{\tau}}_{\tau|\tau-1}-\sum_{\tau=1}^{T}\mathcal{L}({\mbox{\boldmath$\phi$}}^{i_{\tau}\\!\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}^{i_{\tau}}_{*};y_{\tau})$
$\displaystyle\leq\sum_{m=1}^{M}\left(\frac{\|{\mbox{\boldmath$\theta$}}^{m}_{*}\|^{2}}{2\sigma_{\theta^{m}}^{2}}+n_{\text{RF}}\log\left(1+\frac{T_{m}c\sigma_{\theta^{m}}^{2}}{2n_{\text{RF}}}\right)\right)$
$\displaystyle\leq\sum_{m=1}^{M}\left(\frac{\|{\mbox{\boldmath$\theta$}}^{m}_{*}\|^{2}}{2\sigma_{\theta^{m}}^{2}}+n_{\text{RF}}\log\left(1+\frac{T_{m}c\sigma_{\theta^{*}}^{2}}{2n_{\text{RF}}}\right)\right)$
(78)

where $\sigma_{\theta^{*}}^{2}=\max_{m\in\mathcal{M}}\sigma_{\theta^{m}}^{2}$.
Since $\log(\cdot)$ is a concave function, it follows with
$T=\sum_{m=1}^{M}T_{m}$ that

$\displaystyle\sum_{m=1}^{M}\log\left(1+\frac{T_{m}c\sigma_{\theta^{*}}^{2}}{2n_{\text{RF}}}\right)\leq
M\log\left(1+\frac{Tc\sigma_{\theta^{*}}^{2}}{2n_{\text{RF}}}\right)$

which, upon plugging into (78), finalizes the proof of Lemma 3.

###  9.5 Proof of Theorem 2.

For a given $\mathbf{i}_{T}$, expert $m$ relies on batch data
$\\{\mathbf{x}_{\tau},y_{\tau},\tau\in\mathcal{T}_{m}\\}$ in hindsight to
learn the benchmark function $\hat{f}^{m}(\mathbf{x})$ in the RKHS and the RF-
based one $\check{f}_{*}^{m}(\mathbf{x})$ with parameter vector
${\mbox{\boldmath$\theta$}}^{m}_{*}$. Following (62)–(66) in Sec. 9.2, it
holds true with probability at least
$1-2^{8}(\frac{\sigma_{m}}{\epsilon})^{2}\exp\left(\frac{-n_{\text{RF}}\epsilon^{2}}{4d+8}\right)$
that

$\displaystyle\sum_{\tau\in\mathcal{T}_{m}}\\!\\!\\!\mathcal{L}\\!\left({\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){{\mbox{\boldmath$\theta$}}^{m}_{*}};y_{\tau}\\!\right)\\!-\\!\\!\\!\\!\sum_{\tau\in\mathcal{T}_{m}}\\!\\!\\!\mathcal{L}\\!\left(\hat{f}^{m}(\mathbf{x}_{\tau});y_{\tau}\\!\right)\\!\leq\\!\epsilon
LT_{m}C^{\prime}$

where
$C^{\prime}:=\underset{m\in\mathcal{M}}{\max}\sum_{\tau\in\mathcal{T}_{m}}\sigma_{\theta^{m}}^{2}|\hat{\alpha}_{\tau}^{m}|$.
Summing the above inequality over $m\in\mathcal{M}$ leads to

$\displaystyle\sum_{\tau=1}^{T}\\!\mathcal{L}\\!\left({\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){{\mbox{\boldmath$\theta$}}^{m}_{*}};y_{\tau}\\!\right)\\!-\\!\\!\\!\sum_{\tau=1}^{T}\\!\mathcal{L}\\!\left(\hat{f}^{m}(\mathbf{x}_{\tau});y_{\tau}\\!\right)\\!\leq\\!\epsilon
LTC^{\prime}\;$ (79)

holding true with probability at least
$1-2^{8}(\frac{\sigma_{*}}{\epsilon})^{2}\exp\left(\frac{-n_{\text{RF}}\epsilon^{2}}{4d+8}\right)$,
where $\sigma_{*}^{2}:=\max_{m\in\mathcal{M}}\sigma_{m}^{2}$.

Meanwhile, adapting the result in (68) to expert $m$ possessing data
$\mathcal{D}_{m}$, yields the ensuing inequality concerning
${\mbox{\boldmath$\theta$}}^{m}_{*}$

$\displaystyle\|{\mbox{\boldmath$\theta$}}^{m}_{*}\|^{2}:=\left\|\sigma_{\theta^{m}}^{2}\sum_{\tau\in\mathcal{T}_{m}}\hat{\alpha}^{m}_{\tau}{\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{\tau})\right\|^{2}$
(80)
$\displaystyle=\sigma_{\theta^{m}}^{4}\sum_{\tau\in\mathcal{T}_{m}}\sum_{\tau^{\prime}\in\mathcal{T}_{m}}\hat{\alpha}^{m}_{\tau}\hat{\alpha}^{m}_{\tau^{\prime}}{\mbox{\boldmath$\phi$}}^{m\top}_{\mathbf{v}}(\mathbf{x}_{\tau}){\mbox{\boldmath$\phi$}}^{m}_{\mathbf{v}}(\mathbf{x}_{\tau^{\prime}})\leq(1+\epsilon){C^{\prime}}^{2}$

which, in conjunction with Lemmas 2-3 and (79), readily proves Theorem 2.

## Acknowledgments

The authors would like to thank the anonymous reviewers for their constructive
feedback. We also gratefully acknowledge the support from NSF grant 1901134.

## References

  * [1] C. E. Rasmussen and C. K. Williams, _Gaussian processes for machine learning_. MIT press Cambridge, MA, 2006\.
  * [2] I. Kononenko, “Machine learning for medical diagnosis: history, state of the art and perspective,” _Artif. Intel. in Medicine_ , vol. 23, no. 1, pp. 89–109, 2001.
  * [3] K. Cutajar, M. Osborne, J. Cunningham, and M. Filippone, “Preconditioning kernel matrices,” _Proc. Int. Conf. Mach. Learn._ , pp. 2529–2538, 2016\.
  * [4] K. Wang, G. Pleiss, J. Gardner, S. Tyree, K. Q. Weinberger, and A. G. Wilson, “Exact Gaussian processes on a million data points,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 14 648–14 659, 2019.
  * [5] E. Snelson and Z. Ghahramani, “Sparse Gaussian processes using pseudo-inputs,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 1257–1264, 2006.
  * [6] M. Titsias, “Variational learning of inducing variables in sparse Gaussian processes,” _Proc. Int. Conf. Artif. Intel. and Stats._ , pp. 567–574, 2009\.
  * [7] M. Lázaro-Gredilla, J. Quiñonero Candela, C. E. Rasmussen, and A. Figueiras-Vidal, “Sparse spectrum Gaussian process regression,” _J. Mach. Learn. Res._ , vol. 11, no. Jun, pp. 1865–1881, 2010.
  * [8] A. Rahimi and B. Recht, “Random features for large-scale kernel machines,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 1177–1184, 2008.
  * [9] A. Gijsberts and G. Metta, “Real-time model learning using incremental sparse spectrum Gaussian process regression,” _Neural Networks_ , vol. 41, pp. 59–69, 2013.
  * [10] E. Hazan, “Introduction to online convex optimization,” _Foundations and Trends® in Optimization_ , vol. 2, no. 3-4, pp. 157–325, 2016\.
  * [11] N. Cesa-Bianchi and G. Lugosi, _Prediction, learning, and games_. Cambridge University press, 2006.
  * [12] C.-A. Cheng and B. Boots, “Incremental variational sparse Gaussian process regression,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 4410–4418, 2016\.
  * [13] T. D. Bui, C. Nguyen, and R. E. Turner, “Streaming sparse Gaussian process approximations,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 3299–3307, 2017.
  * [14] S. M. Kakade, M. W. Seeger, and D. P. Foster, “Worst-case bounds for Gaussian process models,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 619–626, 2006.
  * [15] E. Gilboa, Y. Saatçi, and J. P. Cunningham, “Scaling multidimensional inference for structured Gaussian processes,” _IEEE Trans. Pattern Anal. Mach. Intel._ , vol. 37, no. 2, pp. 424–436, 2013.
  * [16] J. P. Cunningham, K. V. Shenoy, and M. Sahani, “Fast Gaussian process methods for point process intensity estimation,” _Proc. Int. Conf. Mach. Learn._ , pp. 192–199, 2008.
  * [17] J. Quiñonero-Candela and C. E. Rasmussen, “A unifying view of sparse approximate Gaussian process regression,” _J. Mach. Learn. Res._ , vol. 6, no. Dec, pp. 1939–1959, 2005.
  * [18] Y. Gal and R. Turner, “Improving the Gaussian process sparse spectrum approximation by representing uncertainty in frequency inputs,” _Proc. Int. Conf. Mach. Learn._ , 2015.
  * [19] T. N. Hoang, Q. M. Hoang, and B. K. H. Low, “A unifying framework of anytime sparse Gaussian process regression models with stochastic variational inference for big data.” _Proc. Int. Conf. Mach. Learn._ , pp. 569–578, 2015\.
  * [20] H.-M. Kim, B. K. Mallick, and C. Holmes, “Analyzing nonstationary spatial data using piecewise Gaussian processes,” _J. of the Amer. Stat. Assoc._ , vol. 100, no. 470, pp. 653–668, 2005.
  * [21] V. Tresp, “A Bayesian committee machine,” _Neural Computation_ , vol. 12, no. 11, pp. 2719–2741, 2000.
  * [22] M. P. Deisenroth and J. W. Ng, “Distributed Gaussian processes,” in _Proc. of Intl. Conf. on Machine Learning_ , 2015.
  * [23] C. E. Rasmussen and Z. Ghahramani, “Infinite mixtures of Gaussian process experts,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 881–888, 2002.
  * [24] E. Meeds and S. Osindero, “An alternative infinite mixture of Gaussian process experts,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 883–890, 2006.
  * [25] M. Trapp, R. Peharz, F. Pernkopf, and C. E. Rasmussen, “Deep structured mixtures of Gaussian processes,” _Proc. Int. Conf. Artif. Intel. and Stats._ , pp. 2251–2261, 2020.
  * [26] C. A. Micchelli and M. Pontil, “Learning the kernel function via regularization,” _J. Mach. Learn. Res._ , vol. 6, no. Jul, pp. 1099–1125, 2005.
  * [27] M. A. Alvarez, L. Rosasco, N. D. Lawrence _et al._ , “Kernels for vector-valued functions: A review,” _Foundations and Trends® in Machine Learning_ , vol. 4, no. 3, pp. 195–266, 2012\.
  * [28] J. Lu, S. C. Hoi, J. Wang, P. Zhao, and Z.-Y. Liu, “Large scale online kernel learning,” _J. Mach. Learn. Res._ , vol. 17, no. 1, pp. 1613–1655, 2016\.
  * [29] R. Jin, S. C. Hoi, and T. Yang, “Online multiple kernel learning: Algorithms and mistake bounds,” in _Proc. of Intl. Conf. on Algorithmic Learning Theory_ , 2010, pp. 390–404.
  * [30] Y. Shen, T. Chen, and G. B. Giannakis, “Random feature-based online multi-kernel learning in environments with unknown dynamics,” _J. Mach. Learn. Res._ , vol. 20, no. 1, pp. 773–808, 2019.
  * [31] N. Lawrence, “Probabilistic non-linear principal component analysis with Gaussian process latent variable models,” _J. Mach. Learn. Res._ , vol. 6, no. Nov, pp. 1783–1816, 2005.
  * [32] N. D. Lawrence, “Learning for larger datasets with the Gaussian process latent variable model,” _Proc. Int. Conf. Artif. Intel. and Stats._ , pp. 243–250, 2007.
  * [33] A. C. Damianou, M. K. Titsias, and N. D. Lawrence, “Variational inference for latent variables and uncertain inputs in Gaussian processes,” _J. Mach. Learn. Res._ , vol. 17, no. 1, pp. 1425–1486, 2016.
  * [34] A. Yao, J. Gall, L. V. Gool, and R. Urtasun, “Learning probabilistic non-linear latent variable models for tracking complex activities,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 1359–1367, 2011.
  * [35] Y. Wang, M. A. Brubaker, B. Chaib-draa, and R. Urtasun, “Bayesian filtering with online Gaussian process latent variable models.” in _Uncertainty in Artif. Intel._ , 2014, pp. 849–857.
  * [36] X. Qin, P. Blomstedt, and S. Kaski, “Scalable Bayesian non-linear matrix completion,” _arXiv preprint arXiv:1908.01009_ , 2019.
  * [37] M. Ghashami, D. J. Perry, and J. Phillips, “Streaming kernel principal component analysis,” _Proc. Int. Conf. Artif. Intel. and Stats._ , pp. 1365–1374, 2016.
  * [38] M. E. Tipping and C. M. Bishop, “Mixtures of probabilistic principal component analyzers,” _Neural Computation_ , vol. 11, no. 2, pp. 443–482, 1999.
  * [39] A. Bellas, C. Bouveyron, M. Cottrell, and J. Lacaille, “Model-based clustering of high-dimensional data streams with online mixture of probabilistic PCA,” _Advances in Data Analysis and Classification_ , vol. 7, no. 3, pp. 281–300, 2013.
  * [40] Q. Lu, G. Karanikolas, Y. Shen, and G. B. Giannakis, “Ensemble Gaussian processes with spectral features for online interactive learning with scalability,” _Proc. Int. Conf. Artif. Intel. and Stats._ , pp. 1910–1920, 2020.
  * [41] C. Richard, J. C. M. Bermudez, and P. Honeine, “Online prediction of time series data with kernels,” _IEEE Trans. Sig. Process._ , vol. 57, no. 3, pp. 1058–1067, 2008.
  * [42] N. Xu, K. H. Low, J. Chen, K. K. Lim, and E. B. Ozgul, “GP-localize: Persistent mobile robot localization using online sparse Gaussian process observation model,” in _AAAI Conf. on Artif. Intel._ , 2014.
  * [43] S. M. Kakade and A. Y. Ng, “Online bounds for Bayesian algorithms,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 641–648, 2005.
  * [44] C. K. Williams and D. Barber, “Bayesian classification with Gaussian processes,” _IEEE Trans. Pattern Anal. Mach. Intel._ , vol. 20, no. 12, pp. 1342–1351, 1998.
  * [45] K. P. Murphy, _Machine Learning: A Probabilistic Perspective_. MIT press, 2012.
  * [46] Y. A. Malkov and D. A. Yashunin, “Efficient and robust approximate nearest neighbor search using hierarchical navigable small world graphs,” _IEEE Trans. Pattern Anal. Mach. Intel._ , vol. 42, no. 4, pp. 824–836, 2020.
  * [47] S. De Vito, E. Massera, M. Piga, L. Martinotto, and G. Di Francia, “On field calibration of an electronic nose for benzene estimation in an urban pollution monitoring scenario,” _Sensors and Actuators B: Chemical_ , vol. 129, no. 2, pp. 750–757, 2008.
  * [48] G. Rätsch, “Ensemble learning methods for classification,” April 1998, diploma thesis (in german). [Online]. Available: http://www.first.gmd.de/ raetsch/diplom.ps.gz
  * [49] T. G. Dietterich, A. N. Jain, R. H. Lathrop, and T. Lozano-Perez, “A comparison of dynamic reposing and tangent distance for drug activity prediction,” _Proc. Adv. Neural Inf. Process. Syst._ , pp. 216–223, 1994\.
  * [50] V. G. Sigillito, S. P. Wing, L. V. Hutton, and K. B. Baker, “Classification of radar returns from the ionosphere using neural networks,” _Johns Hopkins APL Technical Digest_ , vol. 10, no. 3, pp. 262–266, 1989.
  * [51] J. J. Hull, “A database for handwritten text recognition research,” _IEEE Trans. Pattern Anal. Mach. Intel._ , vol. 16, no. 5, pp. 550–554, 1994\.
  * [52] C. M. Bishop and G. D. James, “Analysis of multiphase flows using dual-energy gamma densitometry and neural networks,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 327, no. 2-3, pp. 580–593, 1993.
  * [53] F. Kawala, A. Douzal-Chouakria, E. Gaussier, and E. Dimert, “Prédictions d’activité dans les réseaux sociaux en ligne,” in _4ième conférence sur les modèles et l’analyse des réseaux : Approches mathématiques et informatiques_ , France, Oct. 2013, p. 16.
  * [54] D. Dua and C. Graff, “UCI machine learning repository,” 2017. [Online]. Available: http://archive.ics.uci.edu/ml
  * [55] C. E. Rasmussen and H. Nickisch, “Gaussian processes for machine learning (gpml) toolbox,” _J. Mach. Learn. Res._ , vol. 11, no. Nov, pp. 3011–3015, 2010.
  * [56] D. A. Ross, J. Lim, R.-S. Lin, and M.-H. Yang, “Incremental learning for robust visual tracking,” _Intl. J. of Comput. Vis._ , vol. 77, no. 1-3, pp. 125–141, 2008.
  * [57] A. Levey and M. Lindenbaum, “Sequential Karhunen-Loeve basis extraction and its application to images,” _IEEE Trans. Image Process._ , vol. 9, no. 8, pp. 1371–1374, 2000.
  * [58] B. Schölkopf, A. Smola, and K.-R. Müller, “Nonlinear component analysis as a kernel eigenvalue problem,” _Neural Computation_ , vol. 10, no. 5, pp. 1299–1319, 1998.

Figure 9: IE-GP inference on a synthetic dataset. Shaded regions indicate
$95\%$ confidence intervals. Figure 10: IE-GP cumulative static regret (34) on
a synthetic dataset. Figure 11: SIE-GP average switching regret,
$\frac{\mathcal{R}^{\rm SW}(T)}{T}$, versus $T$.

##  Appendix A Synthetic tests

This section contains tests with synthetic data to validate the expected
uncertainty quantification performance of IE-GP and the regret bounds of
(S)IE-GP. To start with, scalar input data $\\{x_{t}\\}_{t=1}^{30}$ were
randomly drawn from a standardized normal distribution, and outputs were
generated as $y_{t}=\sin(2x_{t})+\sin(3x_{t})+\epsilon_{t}$, where
$\epsilon_{t}\sim\mathcal{N}(0,0.01)$. The inferred mean function as well as
(approximate) $95\%$ confidence intervals are shown in Fig. 9. As expected,
regions populated with training examples correspond to tighter confidence
bands relative to unpopulated ones.

To validate the regret bound for IE-GP (cf. (34)), datasets of increasing size
$T$ were generated from the aforementioned synthetic model, albeit with
$x_{t}\sim\mathcal{N}(0,100)$. It is evident from Fig. 10 that the regret can
be upper bounded by $\mathcal{O}(\log T)$, as predicted by the Theorem 1. To
further corroborate the regret bound of SIE-GP, synthetic switching datasets
with growing $T$ were generated, where the first $T/2$ noiseless outputs were
drawn from $\mathcal{GP}(0,\kappa_{1})$ and the second half from
$\mathcal{GP}(0,\kappa_{2})$. Here, $\kappa_{1}$ and $\kappa_{2}$ were two RBF
kernels whose magnitudes are $1$ and characteristic lengthscales are $0.01$
and $100$, respectively. The observed outputs were obtained by adding noises
sampled from zero-mean Gaussian distribution with variance $1$. After running
SIE-GP with these two GP models, the average switching regret (39) was plotted
as a function of $T$ in Fig. 11, where the convergence behavior of the curve
agrees well with Theorem 2.