# Weighted Ensembles for Active Learning with Adaptivity

Konstantinos D. Polyzos , Qin Lu, and Georgios B. Giannakis
University of Minnesota
{polyz003, qlu<EMAIL_ADDRESS>The work of Konstantinos D. Polyzos was also
supported by the Onassis Foundation Scholarship.

###### Abstract

Labeled data can be expensive to acquire in several application domains,
including medical imaging, robotics, and computer vision. To efficiently train
machine learning models under such high labeling costs, active learning (AL)
judiciously selects the most informative data instances to label on-the-fly.
This active sampling process can benefit from a statistical function model,
that is typically captured by a Gaussian process (GP). While most GP-based AL
approaches rely on a single kernel function, the present contribution
advocates an _ensemble_ of _GP models_ with weights adapted to the labeled
data collected incrementally. Building on this novel EGP model, a suite of
acquisition functions emerges based on the uncertainty and disagreement rules.
An adaptively weighted _ensemble_ of EGP-based _acquisition functions_ is also
introduced to further robustify performance. Extensive tests on synthetic and
real datasets showcase the merits of the proposed EGP-based approaches with
respect to the single GP-based AL alternatives.

##  1 Introduction

Identifying machine learning models usually relies on a sufficient number of
labeled input-output data pairs, which may not be feasible due to labeling
costs or privacy concerns in a number of application domains, ranging from
healthcare [7], computer vision [9], to robotics [35]. To cope with sampling
cost constraints, active learning (AL) selects a set of most informative data
to label incrementally [32, 37]. This selection process can benefit from a
statistical model for the learning function $f({\bf x})$ that maps each input
feature vector ${\bf x}_{\tau}$ to the output/label $y_{\tau}$ [2]. Capable of
learning nonlinear functions with uncertainty quantification in a sample-
efficient fashion, Gaussian processes (GPs) are widely adopted to model the
aforementioned $f$ in AL; see e.g., [12, 9]. Given labeled samples
$\\{\mathbf{x}_{\tau},y_{\tau}\\}_{\tau=1}^{t}$ collected in the set ${\cal
L}_{t}$, GP modeling yields a function posterior probability density function
(pdf) $p(f({\bf x})|{\cal L}_{t})$. For regression, the latter is Gaussian
with mean and variance available in closed form, while the model uncertainty
captured by the variance is essential to guide the selection of subsequent
instances to be queried.

The expressiveness of GP models depends on how well the chosen covariance
(kernel) is adapted to the data at hand. Typically, GP-based AL approaches
rely on a single kernel with preselected form, which may exhibit limited
expressiveness for AL, where labeled data are collected incrementally. On par
with this online interactive operation, a key desideratum is a more expressive
function space to improve adaptation to the labeled instances on-the-fly.
Besides expressiveness of the function model, the AL performance is critically
affected by the data acquisition rule. For regression, GP-AL typically
leverages the function variance to select the next instance to be queried.
Devising alternative acquisition rules for more expressive GP models is an
open problem. Clearly, with multiple acquisition rules available, devising a
strategy to combine and adapt them will be conducive to robust performance
across tasks.

To address the aforementioned challenges, our contributions here can be
summarized in the following three aspects.

  * c1)

Relative to GP-based AL that relies on a single GP with a preselected kernel,
we introduce _weighted ensembles_ (E) of GPs for enhanced expressiveness with
weights capturing the contributions of individual GP models, and adapting to
labeled data collected on-the-fly.

  * c2)

Utilizing the novel EGP model, we devise a suite of acquisition functions
(AFs), including _weighted ensembles_ of AFs that further robustify
performance.

  * c3)

Our thorough tests on synthetic and real data corroborate the impressive
merits of GP and AF ensembles.

##  2 Related works

This section outlines the context of the present work.

Statistical models for AL. The performance of model-based AL depends
critically on the chosen statistical model. GPs are widely used because they
come with uncertainty quantification and sample efficiency [12, 9]. However,
most of the existing approaches rely on a preselected GP kernel with limited
expressiveness, which may fall short in characterizing the incrementally
collected data in AL. Broadening the scope of a single GP, a mixture of GPs is
viable by training each GP component on a subset of labeled data [40]. This GP
mixture model can account for multi-modalities in the function space, but
needs to be refitted per iteration using nontrivial variational methods. On
the other hand, the proposed EGP model requires minimal refitting efforts, and
trains each GP on the whole labeled set. Besides GP-based statistical models,
Bayesian convolutional networks have been leveraged for image data [5].

Acquisition rules for AL. The acquisition function (AF) also plays a
performance - critical role in AL. Different AFs have been devised based on
distinct selection criteria, learning tasks, as well as whether or not a
statistical model is capitalized on; see e.g., [32]. How to appropriately
select the AF from the available choices usually calls for domain expertise.
It has been shown empirically that no acquisition rule excels in all tasks.
This motivates a strategy that can combine and adapt candidate selection
rules. In the closely related context of Bayesian optimization, multiple AFs
have been combined to robustify performance across tasks [6, 34]. Such AF
ensembles however, have not been explored for AL.

Kernel selection for GPs. Adaptively choosing the form of the kernel from
training data has been reported for conventional GP learning; see, e.g., [36,
3, 11, 22]. These approaches usually operate in a batch offline mode, and rely
on a large number of samples – what discourages their use for AL, where data
are not only acquired online, but are also scarce due to the potentially high
labeling cost. Recently, an online scalable kernel selection approach has been
introduced that leverages a set of GP experts in an online Bayesian model
averaging fashion [17, 18]. Although the notion of ensembles of GP experts has
been explored in unsupervised learning [10], graph-guided semi-supervised
learning [25, 26, 27], reinforcement learning [16, 15, 28] and Bayesian
optimization [19], to the best of our knowledge it has not yet been utilized
in the AL context, which also entails the AF design.

##  3 Preliminaries

This section outlines the motivation and preliminaries for the AL approach of
interest.

Typical learning approaches boil down to estimating the mapping $f(\cdot)$
that relates the $d\times 1$ input feature vector $\mathbf{x}_{\tau}$ to the
output $y_{\tau}$ (that is either a real number in regression or it belongs to
a finite alphabet in classification) as $\mathbf{x}_{\tau}\rightarrow
f(\mathbf{x}_{\tau})\rightarrow y_{\tau}$. This estimation task relies on a
sufficient number of labeled training samples
$\\{\mathbf{x}_{\tau},y_{\tau}\\}_{\tau=1}^{T}$. In several applications
however, the input can be readily obtained whereas the corresponding label can
be expensive to acquire due to sampling costs or privacy concerns. In
healthcare for instance, many labels describing the medical condition of
patients may not be revealed to preserve confidentiality. Faced with this
challenge, one can resort to the AL paradigm, which judiciously and
proactively selects the most informative instances to label so that the sought
mapping can be inferred in a sample-efficient manner.

AL starts with a small-size set
$\mathcal{L}_{0}:=\\{(\mathbf{x}_{\tau},y_{\tau})\\}_{\tau=-L_{0}+1}^{0}$ of
labeled samples,111The negative instance index here is used for notational
brevity as more labeled data are included next. and a larger collection of
unlabeled features
$\mathcal{U}_{0}:=\\{\mathbf{x}_{\tau}^{u}\\}_{\tau=1}^{U_{0}}$ ($L_{0}\ll
U_{0}$). Given corresponding sets $\mathcal{L}_{t}$ and $\mathcal{U}_{t}$ up
to index $t>0$, model-based AL entails a statistical function model, namely
the pdf $p(f({\bf x})|{\cal L}_{t})$. The latter is utilized by the so-termed
acquisition function (AF) $\alpha(\cdot)$ to select the instance
$\mathbf{x}_{t+1}$ from the corresponding unlabeled set $\mathcal{U}_{t}$, as
[2]

$\displaystyle{\bf
x}_{t+1}=\underset{\mathbf{x}\in\mathcal{U}_{t}}{\arg\max}\;\;\alpha(\mathbf{x};{\cal
L}_{t})\;.$ (1)

Intuitively, $\alpha$ is chosen to guide exploration of the space $f(\cdot)$
belongs to, which hinges on quantifying the uncertainty of the belief model
$p(f({\bf x})|{\cal L}_{t})$. Upon querying an oracle for the associated label
$y_{t+1}$, the labeled set is then augmented with the new pair and the feature
vector is removed from the unlabeled set, that is ${\cal L}_{t+1}:={\cal
L}_{t}\cup\\{(\mathbf{x}_{t+1},y_{t+1})\\}$ and ${\cal U}_{t+1}:={\cal
U}_{t}\setminus\\{\mathbf{x}_{t+1}\\}$. Apparently, the two critical choices
are the model for $f$, and the AF design for $\alpha$. Focusing on the
regression task, we will outline the GP-based model for $f$, and the
associated acquisition rules next.

###  3.1 GP-based active learning

GPs estimate a nonparametric function model in a sample-efficient manner,
while also offering quantification of the associated model uncertainty [29].
Sample efficiency justifies their wide adoption in AL. The GP model postulates
$f$ as being randomly drawn from a GP prior; that is
$f\sim\mathcal{GP}(0,\kappa(\mathbf{x},\mathbf{x}^{\prime}))$ with
$\kappa(\mathbf{x},\mathbf{x}^{\prime})$ being a positive-definite kernel
function that measures the pairwise similarity between two distinct inputs
$\mathbf{x}$ and $\mathbf{x}^{\prime}$. With ⊤ denoting transposition, this
means that the random vector $\mathbf{f}_{t}:=[f(\mathbf{x}_{1})\ldots
f(\mathbf{x}_{t})]^{\top}$ consisting of the function evaluations at instances
$\mathbf{X}_{t}:=\left[\mathbf{x}_{1}\ldots\mathbf{x}_{t}\right]^{\top}$ is
Gaussian distributed as [29]

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

where $\mathbf{K}_{t}$ denotes the $t\times t$ covariance matrix whose $(i,j)$
entry is $[{\bf K}_{t}]_{i,j}={\rm
cov}(f(\mathbf{x}_{i}),f(\mathbf{x}_{j})):=\kappa(\mathbf{x}_{i},\mathbf{x}_{j})$.

Next, the (possibly noisy) output data $\mathbf{y}_{t}:=[y_{1}\cdots
y_{t}]^{\top}$ are related to the function evaluations $\mathbf{f}_{t}$
through the likelihood $p(\mathbf{y}_{t}|\mathbf{f}_{t};\mathbf{X}_{t})$ which
is assumed to be factored as
$p(\mathbf{y}_{t}|\mathbf{f}_{t};\mathbf{X}_{t})=\prod_{\tau=1}^{t}p(y_{\tau}|f(\mathbf{x}_{\tau}))$
with $p(y_{\tau}|f(\mathbf{x}_{\tau}))$ denoting the per-datum known
likelihood factor. The latter holds in the regression setup, where
$p(y_{\tau}|f(\mathbf{x}_{\tau}))=\mathcal{N}(y_{\tau};f(\mathbf{x}_{t}),\sigma_{n}^{2})$,
so long as $y_{\tau}$ is expressed as $y_{\tau}=f_{\tau}+n_{\tau}$ with
$n_{\tau}\sim\mathcal{N}(n_{\tau};0,\sigma_{n}^{2})$ being white Gaussian
noise uncorrelated across instances. For the regression task, predicting the
function value of a generic input $\mathbf{x}$ is carried out after writing
the joint pdf of $\mathbf{y}_{t}$ and $f({\bf x})$ as

$\displaystyle\begin{bmatrix}\mathbf{y}_{t}\\\ f({\bf
x})\end{bmatrix}\\!\sim\\!\mathcal{N}\\!\left(\\!\mathbf{0}_{t+1},\\!\\!\begin{bmatrix}\mathbf{K}_{t}\\!+\\!\sigma_{n}^{2}\mathbf{I}_{t}&\mathbf{k}_{t}({\bf
x})\\\ \mathbf{k}_{t}^{\top}({\bf
x})&\kappa(\mathbf{x},\mathbf{x})\\!+\\!\sigma_{n}^{2}\end{bmatrix}\\!\right)$

where $\mathbf{k}_{t}({\bf
x}):=[\kappa(\mathbf{x}_{1},\mathbf{x}),\ldots,\kappa(\mathbf{x}_{t},\mathbf{x})]^{\top}$.

With the joint pdf at hand, the posterior pdf of $f({\bf x})$ is

$\displaystyle p(f({\bf x})|\mathbf{y}_{t};\mathbf{X}_{t})=\mathcal{N}(f({\bf
x});\mu_{t}({\bf x}),\sigma_{t}^{2}({\bf x}))$ (2)

where

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

Note that the mean in (3a) is a point prediction of $f({\bf x})$, while the
variance in (3b) quantifies the associated uncertainty. In the AL context,
this uncertainty is used by the acquisition function that selects the next
query instance as

$\displaystyle{\bf
x}_{t+1}=\underset{\mathbf{x}\in\mathcal{U}_{t}}{\arg\max}\;\;\sigma_{t}^{2}({\bf
x})\;.$ (4)

It is worth mentioning that for a Gaussian pdf, (4) is tantamount to
maximizing the entropy [20].

The posterior mean and variance in (3) rely on all $t$ instances in
$\mathbf{X}_{t}$ to form $\mathbf{K}_{t}$, and the associated complexity for
its inversion is thus $\mathcal{O}(t^{3})$. Although this complexity can be
affordable in AL where $t$ is small, it can be further reduced. In addition,
$\mu_{t}({\bf x})$ and $\sigma_{t}^{2}({\bf x})$ require direct access to
$\\{\mathbf{x}_{\tau}\\}_{\tau=1}^{t}$, which may be discouraged due to
privacy concerns as in e.g medical records and financial statements. Further,
GP-based AL relies on a preselected kernel function, which may exhibit limited
expressiveness. These limitations can be ameliorated through our novel
ensemble approach that leverages also random spectral features, as delineated
next.

Algorithm 1 RF-based EGP-AL.

1:Initialization: $\mathcal{U}_{0}$, $\mathcal{L}_{0}$, $\mathcal{K}$;

2:for $t=0,\ldots,T$ do

3: Obtain $p(f({\bf x})|{\cal L}_{t})$ via ${\mbox{\boldmath$\Xi$}}_{t}$;

4: Select ${\bf x}_{t+1}$ based on one from (18)- (20), (22)-(23);

5: Query the oracle to obtain $y_{t+1}$;

6: ${\cal L}_{t+1}\\!\\!:=\\!{\cal
L}_{t}\cup\\{(\mathbf{x}_{t+1},y_{t+1})\\}$, ${\cal U}_{t+1}\\!\\!:=\\!{\cal
U}_{t}\setminus\\{\mathbf{x}_{t+1}\\}$;

7:end for

##  4 Ensemble GPs for AL

The chosen function model affects critically the performance of AL approaches.
Unlike most existing works that rely on a single GP with a preselected kernel,
we advocate an ensemble (E) of $M$ GPs to enhance expressiveness. Each GP has
a distinct kernel function selected from a given dictionary
$\mathcal{K}:=\\{\kappa^{1},\ldots,\kappa^{M}\\}$, that is formed using
kernels of different types and with different hyperparameters. Specifically,
each GP $m\in\mathcal{M}:=\\{1,\ldots,M\\}$ places a unique prior on $f$ as
$f|m\sim\mathcal{GP}(0,\kappa^{m}(\mathbf{x},\mathbf{x}^{\prime}))$. The EGP
prior of $f(\mathbf{x})$ is then a weighted ensemble of the individual GP
priors as

$\displaystyle f(\mathbf{x})\sim\sum_{m=1}^{M}w^{m}_{0}{\cal
GP}(0,\kappa^{m}(\mathbf{x},\mathbf{x}^{\prime})),\;\;\;\;\sum_{m=1}^{M}w^{m}_{0}=1$
(5)

where $w^{m}_{0}:={\rm Pr}(i=m)$ is the prior probability that measures the
contribution of GP model $m$. With labeled data collected on-the-fly, the sum-
product rule allows one to express the EGP-based function posterior pdf as

$\displaystyle{p}(f(\mathbf{x})|\mathcal{L}_{t})\\!$
$\displaystyle=\\!\sum_{m=1}^{M}\\!{\rm
Pr}(i\\!=\\!m|\mathcal{L}_{t}){p}(f(\mathbf{x})|i\\!=\\!m,\mathcal{L}_{t})$
(6)

which is a mixture of posterior GPs with weights $w_{t}^{m}:={\rm
Pr}(i=m|{\cal L}_{t})$ that signify the significance of the GP experts. These
weights thus enable online model adaptation.

To efficiently update this EGP function model across $t$, we will leverage a
parametric function approximant, formed by the so-termed random features
(RFs), as outlined next.

###  4.1 RF-based approximation per GP

The RF-based approximation relies on a shift-invariant kernel
$\bar{\kappa}(\mathbf{x},\mathbf{x}^{\prime})=\bar{\kappa}(\mathbf{x}-\mathbf{x}^{\prime})$
satisfying $\kappa=\sigma_{\theta}^{2}\bar{\kappa}$, which can be expressed as
the inverse Fourier transform of a spectral density
$\pi_{\bar{\kappa}}({\mbox{\boldmath$\zeta$}})$ as [30]

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

where
$\int\pi_{\bar{\kappa}}({\mbox{\boldmath$\zeta$}})d{\mbox{\boldmath$\zeta$}}=1$,
allowing $\pi_{\bar{\kappa}}$ to be deemed as a pdf. Since $\bar{\kappa}$ is
real, the last expectation is equal to
$\mathbb{E}_{\pi_{\bar{\kappa}}}\left[\cos({\mbox{\boldmath$\zeta$}}^{\top}(\mathbf{x}-\mathbf{x}^{\prime}))\right]$;
and after drawing a sufficient number $D$ of independent and identically
distributed (i.i.d.) samples $\\{{\mbox{\boldmath$\zeta$}}_{j}\\}_{j=1}^{D}$
from $\pi_{\bar{\kappa}}({\mbox{\boldmath$\zeta$}})$, kernel $\bar{\kappa}$ is
approximated by

$\displaystyle\check{\bar{\kappa}}(\mathbf{x},\mathbf{x}^{\prime}):=\frac{1}{D}\sum_{j=1}^{D}\cos\left({\mbox{\boldmath$\zeta$}}_{j}^{\top}(\mathbf{x}-\mathbf{x}^{\prime})\right)\;.$
(7)

Defining the $2D\\!\times\\!1$ RF vector as [14]

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

the sample average (7) can be re-written as
$\check{\bar{\kappa}}(\mathbf{x},\mathbf{x}^{\prime})={\mbox{\boldmath$\phi$}}_{{\mbox{\boldmath$\zeta$}}}^{\top}(\mathbf{x}){\mbox{\boldmath$\phi$}}_{{\mbox{\boldmath$\zeta$}}}(\mathbf{x}^{\prime})$.
This allows the _parametric linear_ function

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

to have an approximate GP prior. Such an RF-based parametric function readily
yields an efficient model update by propagating the posterior pdf
$p({\mbox{\boldmath$\theta$}}|\mathbf{y}_{t};\mathbf{X}_{t})=\mathcal{N}({\mbox{\boldmath$\theta$}};\hat{{\mbox{\boldmath$\theta$}}}_{t},{\mbox{\boldmath$\Sigma$}}_{t})$
per slot $t$ in a recursive Bayes fashion. It is also worth stressing that the
model learning step does not require direct access to $\mathbf{x}_{t}$, but
relies only on the RF vector
$\boldsymbol{\phi}_{{\mbox{\boldmath$\zeta$}}}(\mathbf{x}_{t})$, which can be
viewed as an encrypted version of $\mathbf{x}_{t}$. This may be appealing if
$\mathbf{x}_{t}$, which may e.g., comprise private medical data, cannot be
revealed during model learning.

Having outlined the RF-based approximation per GP, we can proceed with
updating our RF-based EGP function model, as labeled instances become
available incrementally.

###  4.2 EGP parametric model updates

When kernels in the dictionary are shift-invariant, the RF vector
${\mbox{\boldmath$\phi$}}_{{\mbox{\boldmath$\zeta$}}}^{m}(\mathbf{x})$ per GP
$m$ can be formed via (8) by first drawing i.i.d. random vectors
$\\{{\mbox{\boldmath$\zeta$}}_{j}^{m}\\}_{j=1}^{D}$ from
$\pi_{\bar{\kappa}}^{m}({\mbox{\boldmath$\zeta$}})$, which is the spectral
density of the standardized kernel $\bar{\kappa}^{m}$. Let
$\sigma_{\theta^{m}}^{2}$ be the kernel magnitude, so that
$\kappa^{m}=\sigma_{\theta^{m}}^{2}\bar{\kappa}^{m}$. The generative model for
the sought function and the noisy output $y$ per GP $m$ are characterized by
the $2D\times 1$ vector ${\mbox{\boldmath$\theta$}}^{m}$ as

$\displaystyle{p}({\mbox{\boldmath$\theta$}}^{m})$
$\displaystyle=\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\mathbf{0}_{2D},\sigma_{\theta^{m}}^{2}\mathbf{I}_{2D})$
$\displaystyle p(f(\mathbf{x}_{\tau})|i=m,{\mbox{\boldmath$\theta$}}^{m})$
$\displaystyle=\delta(f(\mathbf{x}_{\tau})-{\mbox{\boldmath$\phi$}}_{{\mbox{\boldmath$\zeta$}}}^{m\top}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}^{m})$
$\displaystyle p(y_{\tau}|{\mbox{\boldmath$\theta$}}^{m},\mathbf{x}_{\tau})$
$\displaystyle=\mathcal{N}(y_{\tau};{\mbox{\boldmath$\phi$}}_{{\mbox{\boldmath$\zeta$}}}^{m\top}(\mathbf{x}_{\tau}){\mbox{\boldmath$\theta$}}^{m},\sigma_{n}^{2})\;.$
(10)

This parametric form allows one to capture the function posterior pdf per GP
$m$ via
$p({\mbox{\boldmath$\theta$}}^{m}|\mathcal{L}_{t})=\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m},{\mbox{\boldmath$\Sigma$}}_{t}^{m})$,
which together with the weight $w_{t}^{m}$, approximates the EGP function
posterior (6) via

$\displaystyle p(\check{f}({\bf x})|{\cal L}_{t})=\sum_{m=1}^{M}w_{t}^{m}{\cal
N}(\check{f}({\bf x});\check{\mu}_{t}^{m}({\bf
x}),(\check{\sigma}_{t}^{m}({\bf x}))^{2})$ (11)

with

$\displaystyle\check{\mu}_{t}^{m}({\bf x})$
$\displaystyle={\mbox{\boldmath$\phi$}}^{m\top}_{{\mbox{\boldmath$\zeta$}}}(\mathbf{x})\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m}$
(12a) $\displaystyle(\check{\sigma}_{t}^{m}({\bf x}))^{2}$
$\displaystyle={\mbox{\boldmath$\phi$}}^{m\top}_{{\mbox{\boldmath$\zeta$}}}(\mathbf{x}){\mbox{\boldmath$\Sigma$}}^{m}_{t}{\mbox{\boldmath$\phi$}}^{m}_{{\mbox{\boldmath$\zeta$}}}(\mathbf{x})\;.$
(12b)

Next, we will see how RF-based EGP propagates the function model by updating
across $t$ the parameter set

$\displaystyle{\mbox{\boldmath$\Xi$}}_{t}:=\\{w_{t}^{m},\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m},{\mbox{\boldmath$\Sigma$}}_{t}^{m},m\in\mathcal{M}\\}\;.$
(13)

Upon acquiring the newly labeled pair $\\{\mathbf{x}_{t+1},y_{t+1}\\}$, the
updated weight $w_{t+1}^{m}:={\rm Pr}(i=m|\mathcal{L}_{t+1})$ can be obtained
per GP $m$ via Bayes’ rule as

$\displaystyle w_{t+1}^{m}$ $\displaystyle=\frac{{\rm
Pr}(i=m|\mathcal{L}_{t}){p}(y_{t+1}|\mathbf{x}_{t+1},i=m,\mathcal{L}_{t})}{{p}(y_{t+1}|\mathbf{x}_{t+1},\mathcal{L}_{t})}.$

Since the per-model predictive likelihood is given by

$\displaystyle{p}(y_{t+1}|i\\!=\\!m,\mathcal{L}_{t},\mathbf{x}_{t+1})\\!\\!$
$\displaystyle=\\!\\!\int\\!\\!p(y_{t+1}|{\mbox{\boldmath$\theta$}}^{m}\\!,\mathbf{x}_{t+1})p({\mbox{\boldmath$\theta$}}^{m}|\mathcal{L}_{t})d{\mbox{\boldmath$\theta$}}^{m}$
$\displaystyle=\mathcal{N}(y_{t+1};\hat{y}_{t+1|t}^{m},(\sigma_{t+1|t}^{m})^{2})$

with

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

the updated weight can thus be expressed as

$\displaystyle w_{t+1}^{m}$
$\displaystyle=\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)}\;.$
(15)

Bayes’ rule further allows updating the posterior of
${\mbox{\boldmath$\theta$}}^{m}$ as

$\displaystyle{p}({\mbox{\boldmath$\theta$}}^{m}|\mathcal{L}_{t+1})$
$\displaystyle=\frac{{p}({\mbox{\boldmath$\theta$}}^{m}|\mathcal{L}_{t}\\!){p}(y_{t+1}|{\mbox{\boldmath$\theta$}}^{m},\mathbf{x}_{t+1})}{{p}(y_{t+1}|\mathbf{x}_{t+1},i=m,\mathcal{L}_{t})}$
$\displaystyle=\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\hat{{\mbox{\boldmath$\theta$}}}_{t+1}^{m},{\mbox{\boldmath$\Sigma$}}^{m}_{t+1})$
(16)

where the 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}_{{\mbox{\boldmath$\zeta$}}}(\mathbf{x}_{t+1})(y_{t+1}\\!-\\!\hat{y}_{t+1|t}^{m})$
$\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}_{{\mbox{\boldmath$\zeta$}}}\\!(\mathbf{x}_{t+1}){\mbox{\boldmath$\phi$}}^{m\top}_{{\mbox{\boldmath$\zeta$}}}\\!\\!(\mathbf{x}_{t+1}){\mbox{\boldmath$\Sigma$}}^{m}_{t}\\!.$

###  4.3 Acquisition rules for EGP-based AL

Using the EGP posterior in (11), we are ready to devise AFs that select the
next query point based on different rules.

####  4.3.1 Weighted variance

Motivated by (4), the first AF relies on the uncertainty expressed by the
variance. With GP expert $m$ forming the function posterior with variance
$(\check{\sigma}_{t}^{m}({\bf x}))^{2}$, a weighted combination over all the
$M$ experts yields the AF

$\displaystyle\alpha^{\rm wVar}({\bf x};{\cal
L}_{t}):=\sum_{m=1}^{M}w_{t}^{m}(\check{\sigma}_{t}^{m}({\bf x}))^{2}\;.$ (18)

####  4.3.2 Weighted entropy

Relying alternatively on the entropy as the uncertainty measure, one can take
a weighted sum of the entropy values given by the $M$ GP experts, yielding

$\displaystyle\alpha^{\rm wEnt}({\bf x};{\cal
L}_{t}):=\frac{1}{2}\sum_{m=1}^{M}w_{t}^{m}\ln(2\pi(\check{\sigma}_{t}^{m}({\bf
x}))^{2})\;.$ (19)

####  4.3.3 Query-by-Committee (QBC)

Besides capturing uncertainty by variance or entropy, an alternative
disagreement-based AF – QBC, has been reported for classification [33], and
regression using neural networks [13]. With the $M$ GP experts forming a
committee, the novel EGP-based QBC rule is (cf. (12a))

$\displaystyle\alpha^{\rm QBC}({\bf x};{\cal
L}_{t}):=\sum_{m=1}^{M}w_{t}^{m}(\check{\mu}_{t}^{m}({\bf
x})-\bar{\mu}_{t}({\bf x}))^{2}$ (20)

where $\bar{\mu}_{t}({\bf x})$ is the consensus of the committee given by

$\displaystyle\bar{\mu}_{t}({\bf
x})=\sum_{m=1}^{M}w_{t}^{m}\check{\mu}_{t}^{m}({\bf x})\;.$ (21)

Unlike previous QBC approaches that have equal weights per committee member,
the weights $w_{t}^{m}$ in (20) and (21) are generally different across $m$.

####  4.3.4 Variance of GP mixtures

Rather than directly summing per-GP weighted variances in (18), one can
alternatively obtain the variance based on the GP mixture of the function
posterior (cf. (11)) as

$\displaystyle\alpha^{\rm GPM-Var}({\bf x};{\cal L}_{t})$
$\displaystyle:=\sum_{m=1}^{M}w_{t}^{m}((\check{\sigma}_{t}^{m}({\bf
x}))^{2}+(\check{\mu}_{t}^{m}({\bf x})-\bar{\mu}_{t}({\bf x}))^{2}$ (22)

which, interestingly, is the sum of (18) and (20).

####  4.3.5 Entropy of GP mixtures

The last AF is given by the entropy of the GP mixture in (11), which
unfortunately, has no analytic expression. Aiming at a tractable form, we will
resort to its analytic lower bound [8], which can be expressed as

$\displaystyle-\\!\\!\sum_{m=1}^{M}\\!w_{t}^{m}\\!\\!\int\\!\\!{\cal
N}(\check{f}({\bf x});\check{\mu}_{t}^{m}({\bf
x}),\\!(\check{\sigma}_{t}^{m}({\bf x}))^{2})\log p(\check{f}({\bf
x})|\\!{\cal L}_{t})d\check{f}({\bf x})$
$\displaystyle\overset{(a)}{\geq}-\sum_{m=1}^{M}w_{t}^{m}\log\Big{(}\int{\cal
N}(\check{f}({\bf x});\check{\mu}_{t}^{m}({\bf
x}),(\check{\sigma}_{t}^{m}({\bf x}))^{2})$
$\displaystyle\qquad\qquad\qquad\qquad\quad\times p(\check{f}({\bf x})|{\cal
L}_{t})d\check{f}({\bf x})\Big{)}$

where $(a)$ holds due to Jensen’s inequality. Upon obtaining the analytic
expression for the term inside the logarithm, the last AF is then given by

$\displaystyle\alpha^{\rm GPM-Ent}({\bf x};{\cal
L}_{t}):=-\sum_{m=1}^{M}w_{t}^{m}\log\left(\sum_{m^{\prime}=1}^{M}w_{t}^{m^{\prime}}z^{m,m^{\prime}}_{t}\right)$
(23)

with $z^{m,m^{\prime}}$ accounting for the interaction of any two GP models as

$\displaystyle z^{m,m^{\prime}}_{t}$ $\displaystyle:=\int{\cal
N}(\check{f}({\bf x});\check{\mu}_{t}^{m}({\bf
x}),(\check{\sigma}_{t}^{m}({\bf x}))^{2})$ $\displaystyle\quad\times{\cal
N}(\check{f}({\bf x});\check{\mu}_{t}^{m^{\prime}}({\bf
x}),(\check{\sigma}_{t}^{m^{\prime}}({\bf x}))^{2})d\check{f}({\bf x})$
$\displaystyle={\cal N}(\check{\mu}_{t}^{m}({\bf
x});\check{\mu}_{t}^{m^{\prime}}({\bf x}),(\check{\sigma}_{t}^{m}({\bf
x}))^{2}+(\check{\sigma}_{t}^{m^{\prime}}({\bf x}))^{2})\;.$

Based on our novel EGP-based AFs, implementation of the proposed EGP-AL
approach is summarized in Alg. 1.

Algorithm 2 EGP-MultiAFs for AL.

1:Initialization: $\mathcal{L}_{0}$, $\mathcal{U}_{0}$, $\mathcal{V}$,
$\mathcal{K}$;

2:$\boldsymbol{\omega}_{0}=\frac{1}{K}[1,\ldots,1]^{\top}$;

3:for $t=0,1,\ldots,T$ do

4: Obtain EGP ${\mbox{\boldmath$\Xi$}}_{t}$ based on ${\cal L}_{t}$ using
(15)-(16);

5: for $k=1,\ldots,K$ do

6: Obtain instance $\tilde{\mathbf{x}}_{t+1}^{k}\in\mathcal{U}_{t}$ by (24);

7: Obtain pseudo-label $\tilde{y}_{t+1}^{k}$ using
${\mbox{\boldmath$\Xi$}}_{t}$ via (25);

8: Using pseudo pair $\\{\tilde{{\bf x}}_{t+1}^{k},\tilde{y}_{t+1}^{k}\\}$
obtain $\tilde{{\mbox{\boldmath$\Xi$}}}_{t+1}^{k}$;

9: Obtain error $\epsilon_{t+1}^{v,k}$ on $\mathcal{V}$ via (27);

10: end for

11: Update per AF weight using (29);

12: Obtain $\mathbf{x}_{t+1}\in\mathcal{U}_{t}$ by (30);

13: Query the oracle to obtain true label $y_{t+1}$;

14: $\mathcal{L}_{t+1}=\mathcal{L}_{t}\cup(\mathbf{x}_{t+1},y_{t+1})$;

15: $\mathcal{U}_{t+1}=\mathcal{U}_{t}\setminus\\{\mathbf{x}_{t+1}\\}$;

16:end for

##  5 Ensemble of EGP-based AFs

So far, we have introduced a novel EGP-based function model along with several
choices for the AF. In the context of Bayesian optimization though, it is
known that no single AF excels at all tasks [6]. Hence, combining candidate
AFs can intuitively offer robustness and improved performance. To this end, we
will rely on a validation set ${\cal
V}:=\\{(\mathbf{x}_{\tau}^{v},y_{\tau}^{v})\\}_{\tau=1}^{V}$ to evaluate the
performance of different AFs. Similar to EGP, each of the $K$ candidate AFs
will come with a weight (probability) $\omega_{t}^{k}\in[0,1]$ to capture its
contribution per slot $t$, such that $\sum_{k=1}^{K}\omega_{t}^{k}=1$.

Figure 1: NMSE and NPLL versus iterations for Ackley5D, Currin exponential,
Branin, Gramacy and Higdon functions

Upon identifying the RF-based EGP set ${\mbox{\boldmath$\Xi$}}_{t}$ in (13)
using the labeled set ${\cal L}_{t}$ at slot $t$, each AF $k$ selects its
query point $\tilde{{\bf x}}_{t+1}^{k}$ at slot $t+1$ by optimizing the
associated criterion as

$\displaystyle\tilde{{\bf
x}}_{t+1}^{k}=\underset{\mathbf{x}\in\mathcal{U}_{t}}{\arg\max}\;\;\alpha^{k}(\mathbf{x};{\cal
L}_{t})\;.$ (24)

Upon obtaining $\tilde{{\bf x}}_{t+1}^{k}$, AF $k$ constructs a ‘pseudo label’
$\tilde{y}_{t+1}^{k}$ using the EGP parameters in
${\mbox{\boldmath$\Xi$}}_{t}$, as

$\displaystyle\tilde{y}_{t+1}^{k}=\sum_{m=1}^{M}w_{t}^{m}{\mbox{\boldmath$\phi$}}^{m\top}_{{\mbox{\boldmath$\zeta$}}}(\tilde{{\bf
x}}_{t+1}^{k})\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m}\;.$ (25)

This pseudo pair $\\{\tilde{{\bf x}}_{t+1}^{k},\tilde{y}_{t+1}^{k}\\}$ allows
one to leverage (15)– (16) to find the updated EGP parameter vector as

$\displaystyle\tilde{{\mbox{\boldmath$\Xi$}}}_{t+1}^{k}=\\{\tilde{w}_{t+1}^{m,k},\tilde{{\mbox{\boldmath$\theta$}}}_{t+1}^{m,k},\tilde{{\mbox{\boldmath$\Sigma$}}}_{t+1}^{m,k},m\in\mathcal{M}\\}$
(26)

based on which the loss per AF can be evaluated.

To find this loss, AF $k$ capitalizes on
$\tilde{{\mbox{\boldmath$\Xi$}}}_{t+1}^{k}$ in order to obtain the prediction
error at the validation set

$\displaystyle\epsilon_{t+1}^{v,k}=V^{-1}\sum_{\tau=1}^{V}(y_{\tau}^{v}-\hat{y}_{\tau|t+1}^{v,k})^{2}$
(27)

where the predicted label per validation sample $\tau$ is

$\displaystyle\hat{y}_{\tau|t+1}^{v,k}=\sum_{m=1}^{M}\tilde{w}_{t+1}^{m,k}{\mbox{\boldmath$\phi$}}^{m\top}_{{\mbox{\boldmath$\zeta$}}}({\bf
x}_{\tau}^{v})\tilde{{\mbox{\boldmath$\theta$}}}_{t+1}^{m,k}\;.$ (28)

Having available the prediction error over the validation set per AF $k$, the
associated weight can then be updated as

$\displaystyle\omega_{t+1}^{k}=\frac{\omega_{t}^{k}\exp(-\eta\epsilon_{t+1}^{v,k})}{\sum_{k^{\prime}=1}^{K}\omega_{t}^{k^{\prime}}\exp(-\eta\epsilon_{t+1}^{v,k^{\prime}})}$
(29)

where $\eta$ denotes the learning rate. Here, the weight update rule is
similar to that in EGP (cf. (15)), and belongs to the exponential weight
update in online learning with expert advice; see e.g., [1].

Given the updated weights, the next query point is identified by maximizing
the weighted ensemble of AFs, that is

$\displaystyle{\bf
x}_{t+1}=\underset{\mathbf{x}\in\mathcal{U}_{t}}{\arg\max}\sum_{k=1}^{K}\omega_{t+1}^{k}\alpha^{k}(\mathbf{x};{\cal
L}_{t})\;.$ (30)

Upon querying the oracle for the label $y_{t+1}$ of instance ${\bf x}_{t+1}$,
the labeled and unlabeled sets are updated, thus completing one iteration of
the novel “EGP-MultiAFs" approach, that is implemented as listed in Alg. 2.

Figure 2: NMSE (top) and NPLL (bottom) versus iterations for California
housing, Malaria, and Diabetes datasets.

##  6 Numerical tests

In this section, the performance of the proposed EGP-based AL models will be
compared against several benchmark synthetic functions, and it will be tested
with real datasets ranging from biomedical to robotic based ones. Based on the
novel EGP model, the innovative acquisition criteria to be tested are the ones
described in (18) - (23) and (30), which henceforth will be abbreviated as
“EGP-wVar,” “EGP-wEnt,” “EGP-QBC,” “EGP-GPM-Var,” “EGP-GPM-Ent,” and “EGP-
multiAFs,” respectively. The competing baseline is a single GP model that
utilizes the maximum variance (entropy) AF in (4) that has been extensively
used in AL; see e.g., [9, 24, 31]. For all approaches, the few initially
labeled data collected in $\mathcal{L}_{0}$ are utilized to obtain the kernel
hyperparameters per GP expert by maximizing the marginal likelihood using the
sklearn package. The RF-based GP approximate models rely on $D=50$ RFs. For
EGP-based approaches, the kernel dictionary $\mathcal{K}$ consists of radial
basis functions (RBFs) with lengthscales $\\{10^{c}\\}_{c=-4}^{6}$. For the
EGP-multiAFs approach, each $\alpha^{k}(\mathbf{x};{\cal L}_{t})$ in (30) is
divided by its maximum value so that to range between 0 and 1.

The performance of the competing methods is evaluated on a held-out test set
$\mathcal{T}^{e}:=(\mathbf{x}_{\tau}^{e},y_{\tau}^{e})_{\tau=1}^{T^{e}}$
(superscript e stands for evaluation) using two metrics. The first performance
metric is the normalized mean-square error (NMSE) that for iteration $t$, is
given by

$\displaystyle\text{NMSE}_{t}:=\frac{1}{T^{e}}\sum_{\tau=1}^{T^{e}}(\hat{y}_{\tau|t}^{e}-y_{\tau}^{e})^{2}/\sigma_{y}^{2}$

where $\hat{y}_{\tau|t}^{e}$ denotes the point prediction of test instance
$\tau$, and
$\sigma_{y}^{2}:=\mathbb{E}\|\mathbf{y}_{T}^{e}-\mathbb{E}\\{\mathbf{y}_{T}^{e}\\}\|^{2}$,
where $\mathbf{y}_{T}^{e}:=[y_{1}^{e}\ldots y_{T}^{e}]^{\top}$. A second
metric used to assess the associated uncertainty is the negative predictive
log-likelihood (NPLL) expressed as

$\displaystyle\text{NPLL}_{t}:=-\log\;p(\mathbf{y}_{T}^{e}|\mathcal{L}_{t},\mathbf{X}_{T^{e}})$

where the matrix
$\mathbf{X}_{T^{e}}:=[\mathbf{x}_{1}^{e}\dots\mathbf{x}_{T}^{e}]^{\top}$
collects the feature vectors of all $T^{e}$ test instances. All methods are
tested over 10 realizations, and their average performance is reported along
with the corresponding standard deviation. More details about the experimental
set up can be found in Table 1 in the supplementary file.

Figure 3: NMSE (top) and NPLL (bottom) versus iterations for Robot Pushing 3D
and 4D tasks.

###  6.1 Synthetic functions

The tests here are run for known synthetic functions, including Ackley5D,
Currin exponential, Branin, Gramacy and Higdon; see Table 2 in the
supplementary file for their analytic expression. Fig. 1 demonstrates that all
EGP-based approaches with a single AF outperform the GP-Var baseline in the
Currin exponential and Gramacy functions, in terms of NMSE and most have
superior performance in the remaining three datasets. In addition, all EGP-
based methods enjoy the lowest NPLL in four out of five datasets, which
corroborates the merits of having an ensemble of GPs and using them in the
corresponding acquisition criteria. Further exploiting an ensemble of AFs in
the adaptive EGP-multiAFs approach, significantly improves the prediction
performance, and also effectively quantifies the prediction uncertainty, thus
rendering it the best performing approach over all datasets.

###  6.2 Real datasets

All approaches here are tested on California housing, Malaria, and Diabetes
real datasets. The last two deal with real medical data that are well
motivated for AL because of the scarcity of labeled instances emanating from
medical confidentiality. These datasets are described next.

California housing dataset. This set considers 8 features of districts in
California, including demographic and location data, but also more general
features such as average number of rooms and bedrooms per household, based on
which a regression task is formed where the target variable is the median
house price in these districts [23].

Malaria dataset. In this dataset, each instance corresponds to a specific
location where the 8-dimensional feature vector consists of the latitude,
longitude and some bioclimatic characteristics of this location, while the
target variable is the infection rate of Plasmodium falciparum, the parasite
that causes malaria [39, 21].

Diabetes dataset. This dataset considers 10 characteristics of diabetes
patients, including age, sex, body mass index, average blood pressure, and six
blood-related measurements. The target variable is a metric that quantifies
the disease progression in a single year [4].

It is evident in Fig. 2 that all EGP-based approaches markedly outperform the
GP-Var baseline in terms of NMSE and NPLL in the California housing and
Diabetes datasets, showcasing the significance of adopting the EGP model to
estimate the learning function, along with the corresponding AFs. It is worth
mentioning that although most EGP-based approaches are comparable in terms of
prediction error in the California housing dataset, EGP-multiAFs outperforms
all other methods in terms of NPLL. Although the single GP approach is more
effective in the Malaria dataset compared to the EGP-based ones with a single
AF, EGP-multiAFs that judiciously combines the benefits of all AFs, enjoys the
lowest NMSE. This illustrates the significance of properly adjusting AF
weights in an online adaptive fashion. Regarding the NPLL metric, although it
increases per iteration in the Malaria dataset, five out of six EGP-based
approaches achieve lower NPLL compared to the GP-Var baseline.

###  6.3 Robotic tasks

The last experiments focus on a practical robotic task, where a robot pushes
an object to a specific location [38]. Specifically, given as input the robot
location $(r_{\tau}^{x},r_{\tau}^{y})$ and pushing duration $t_{\tau}^{p}$ at
slot $\tau$, the object ends up in a location
$\mathbf{o}_{\tau}:=(o_{\tau}^{x},o_{\tau}^{y})$. We form a regression task
where the goal per slot $\tau$ is to map the $3\times 1$ feature vector
$\mathbf{x}_{\tau}:=[r_{\tau}^{x},r_{\tau}^{y},t_{\tau}^{p}]^{\top}$ to the
target variable $y_{\tau}:=||\mathbf{o}_{\tau}-\mathbf{d}||_{2}$, with
$\mathbf{d}:=[d^{x},d^{y}]$ denoting a pre-defined constant position vector,
yielding the Robot pushing 3D dataset. This is of practical interest in
various robotic problems such as obstacle avoidance, where $y_{\tau}$ is
desired to be greater than a pre-defined threshold $y^{\text{th}}$. Augmenting
the feature vector $\mathbf{x}_{\tau}$ with the additional feature
$r_{\tau}^{\theta}$ that represents the pushing angle, yields the Robot
pushing 4D dataset.

Fig. 3 depicts the NMSE and NPLL at each iteration of all competing AL
approaches for the Robot pushing 3D and 4D tasks respectively. It is evident
that all EGP-based approaches enjoy remarkably lower NMSE and NPLL compared to
the single GP based AL counterpart in both datasets, with the EGP-MultiAFs
consistently being the best-performing one. This implies that in these
practical robotic tasks, the function expressiveness offered by the advocated
EGP model and the ensuing innovative acquisition criteria considerably improve
the prediction performance providing also quantifiable prediction uncertainty.

##  7 Conclusions

This work advocated a weighted ensemble of GPs as the statistical model in AL.
By adapting the weights of individual GPs, the EGP model selects the
appropriate kernel on-the-fly as new labeled data are included incrementally.
Building on the novel EGP model, several AFs have been devised based on
different criteria. Combining the candidate EGP-based AFs with weights being
adjusted in an adaptive manner, the AL performance can be further robustified.
Tests on synthetic functions and real datasets showcase the merits of weighted
ensembles of GPs and AFs in AL.

## References

  * [1] N. Cesa-Bianchi and G. Lugosi.  Prediction, Learning, and Games.  Cambridge University Press, 2006.
  * [2] D. A. Cohn, Z. Ghahramani, and M. I. Jordan.  Active learning with statistical models.  Journal of Artificial Intelligence Research, 4:129–145, 1996.
  * [3] D. Duvenaud, J. Lloyd, R. Grosse, J. Tenenbaum, and G. Zoubin.  Structure discovery in nonparametric regression through compositional kernel search.  Proc. Intl. Conf. Mach. Learn., pages 1166–1174, 2013.
  * [4] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani.  Least angle regression.  The Annals of Statistics, 32(2):407–499, 2004.
  * [5] Y. Gal, R. Islam, and Z. Ghahramani.  Deep Bayesian active learning with image data.  Proc. Intl. Conf. Mach. Learn., pages 1183–1192, 2017.
  * [6] M. Hoffman, E. Brochu, N. de Freitas, et al.  Portfolio allocation for Bayesian optimization.  Proc. Conf. Uncertainty in Artif. Intel., pages 327–336, 2011.
  * [7] S. C. Hoi, R. Jin, J. Zhu, and M. R. Lyu.  Batch mode active learning and its application to medical image classification.  In Proc. Intl. Conf. Mach. Learn., pages 417–424, 2006.
  * [8] M. F. Huber, T. Bailey, H. Durrant-Whyte, and U. D. Hanebeck.  On entropy approximation for Gaussian mixture random vectors.  In IEEE Intl. Conf. on Multisensor Fusion and Integ. for Intell. Syst., pages 181–188, 2008.
  * [9] A. Kapoor, K. Grauman, R. Urtasun, and T. Darrell.  Active learning with Gaussian processes for object categorization.  Proc. Intl. Conf. Comp. Vision, 2007.
  * [10] G. V. Karanikolas, Q. Lu, and G. B. Giannakis.  Online unsupervised learning using ensemble Gaussian processes with random features.  Proc. IEEE Int. Conf. Acoust., Speech, Sig. Process., pages 3190–3194, 2021.
  * [11] H. Kim and Y. W. Teh.  Scaling up the automatic statistician: Scalable structure discovery using Gaussian processes.  Proc. Intl. Conf. Artif. Intel. and Stats., pages 575–584, 2018\.
  * [12] A. Krause, A. Singh, and C. Guestrin.  Near-optimal sensor placements in Gaussian processes: Theory, efficient algorithms and empirical studies.  J. Mach. Learn. Res., 9(2), 2008.
  * [13] A. Krogh and J. Vedelsby.  Neural network ensembles, cross validation, and active learning.  Proc. Adv. Neural Inf. Process. Syst., 7, 1994.
  * [14] M. Lázaro-Gredilla, J. Quiñonero Candela, C. E. Rasmussen, and A. Figueiras-Vidal.  Sparse spectrum Gaussian process regression.  J. Mach. Learn. Res., 11(Jun):1865–1881, 2010.
  * [15] Q. Lu and G. B. Giannakis.  Gaussian process temporal-difference learning with scalability and worst-case performance guarantees.  Proc. IEEE Int. Conf. Acoust., Speech, Sig. Process., pages 3485–3489, 2021.
  * [16] Q. Lu and G. B. Giannakis.  Robust and adaptive temporal-difference learning using an ensemble of Gaussian processes.  arXiv preprint arXiv:2112.00882, 2021.
  * [17] Q. Lu, G. Karanikolas, Y. Shen, and G. B. Giannakis.  Ensemble Gaussian processes with spectral features for online interactive learning with scalability.  Proc. Intl. Conf. Artif. Intel. and Stats., pages 1910–1920, 2020\.
  * [18] Q. Lu, G. V. Karanikolas, and G. B. Giannakis.  Incremental ensemble Gaussian processes.  IEEE Trans. Pattern Anal. Mach. Intel., 2022.
  * [19] Q. Lu, K. D. Polyzos, B. Li, and G. B. Giannakis.  Surrogate modeling for bayesian optimization beyond a single Gaussian process.  arXiv preprint arXiv:2205.14090, 2022.
  * [20] D. J. MacKay.  Information-based objective functions for active data selection.  Neural Computation, 4(4):590–604, 1992.
  * [21] W. J. Maddox, S. Stanton, and A. G. Wilson.  Conditioning sparse variational Gaussian processes for online decision-making.  Proc. Adv. Neural Inf. Process. Syst., 34, 2021.
  * [22] G. Malkomes, C. Schaff, and R. Garnett.  Bayesian optimization for automated model selection.  Proc. Adv. Neural Inf. Process. Syst., 2016.
  * [23] R. K. Pace and R. Barry.  Sparse spatial autoregressions.  Statistics and Probability Letters, 33(3):291–297, 1997.
  * [24] E. Pasolli and F. Melgani.  Gaussian process regression within an active learning scheme.  In IEEE Intl. Geoscience and Remote Sensing Symp., pages 3574–3577, 2011.
  * [25] K. D. Polyzos, Q. Lu, and G. B. Giannakis.  Ensemble gaussian processes for online learning over graphs with adaptivity and scalability.  IEEE Trans. Sig. Process., 2021.
  * [26] K. D. Polyzos, Q. Lu, and G. B. Giannakis.  Graph-adaptive incremental learning using an ensemble of Gaussian process experts.  Proc. IEEE Int. Conf. Acoust., Speech, Sig. Process., pages 5220–5224, 2021.
  * [27] K. D. Polyzos, Q. Lu, and G. B. Giannakis.  Online graph-guided inference using ensemble Gaussian processes of egonet features.  Proc. Asilomar Conf. Sig., Syst., Comput., pages 182–186, 2021\.
  * [28] K. D. Polyzos, Q. Lu, A. Sadeghi, and G. B. Giannakis.  On-policy reinforcement learning via ensemble Gaussian processes with application to resource allocation.  Proc. Asilomar Conf. Sig., Syst., Comput., pages 1018–1022, 2021\.
  * [29] C. E. Rasmussen and C. K. Williams.  Gaussian processes for machine learning.  MIT press Cambridge, MA, 2006.
  * [30] W. Rudin.  Principles of Mathematical Analysis, volume 3.  McGraw-hill New York, 1964.
  * [31] J. Schreiter, D. Nguyen-Tuong, M. Eberts, B. Bischoff, H. Markert, and M. Toussaint.  Safe exploration for active learning with Gaussian processes.  In European Conf. Mach. Learn. and Princ. and Practice of Knowledge Disc. in Databases, 2015.
  * [32] B. Settles.  Active learning.  Synthesis Lectures on Artif. Intel. and Mach. Learn., 6(1):1–114, 2012.
  * [33] H. S. Seung, M. Opper, and H. Sompolinsky.  Query by committee.  In Proc. of the Workshop on Computational Learning Theory, pages 287–294, 1992.
  * [34] B. Shahriari, Z. Wang, M. W. Hoffman, A. Bouchard-Côté, and N. de Freitas.  An entropy search portfolio for Bayesian optimization.  arXiv preprint arXiv:1406.4625, 2014.
  * [35] A. T. Taylor, T. A. Berrueta, and T. D. Murphey.  Active learning in robotics: A review of control principles.  Mechatronics, 77:102576, 2021.
  * [36] T. Teng, J. Chen, Y. Zhang, and B. K. H. Low.  Scalable variational Bayesian kernel selection for sparse Gaussian process regression.  Proc. AAAI Conf. Artif. Intel., 34(04):5997–6004, 2020.
  * [37] S. Tong.  Active Learning: Theory and Applications.  Stanford University, 2001.
  * [38] Z. Wang and S. Jegelka.  Max-value entropy search for efficient Bayesian optimization.  Proc. Intl. Conf. Mach. Learn., pages 3627–3635, 2017.
  * [39] D. J. Weiss, T. C. Lucas, M. Nguyen, A. K. Nandi, D. Bisanzio, K. E. Battle, E. Cameron, K. A. Twohig, D. A. Pfeffer, J. A. Rozier, et al.  Mapping the global prevalence, incidence, and mortality of plasmodium falciparum, 2000–17: a spatial and temporal modelling study.  The Lancet, 394(10195):322–331, 2019.
  * [40] J. Zhao, S. Sun, H. Wang, and Z. Cao.  Promoting active learning with mixtures of Gaussian processes.  Knowledge-Based Systems, 188:105044, 2020.

Supplementary file

Table 1: Additional experimental details Dataset |  $\mathcal{L}_{0}$ size |  $\mathcal{V}$ size |  $\mathcal{U}_{0}$ size |  $\mathcal{T}$ size | $\eta$
---|---|---|---|---|---
Ackley-5D | 10 | 50 | 500 | 100 | 1
Branin | 10 | 50 | 500 | 100 | 100
Currin exponential | 10 | 50 | 500 | 100 | 100
Gramacy | 10 | 50 | 500 | 100 | 100
Higdon | 10 | 50 | 500 | 100 | 100
Diabetes | 15 | 55 | 261 | 111 | 100
Malaria | 40 | 60 | 500 | 680 | 20
Robot pushing 3D | 20 | 50 | 500 | 200 | 20
Robot pushing 4D | 20 | 50 | 500 | 200 | 20
California housing | 50 | 70 | 1000 | 1032 | 0.05
Table 2: Analytical expression of all synthetic functions Function |  Analytical expression
---|---
Ackley-5D |  $-20e^{-0.2\sqrt{(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}+x_{5}^{2})/5}}-e^{(\text{cos}(2\pi x_{1})+\text{cos}(2\pi x_{2})+\text{cos}(2\pi x_{3})+\text{cos}(2\pi x_{4})+\text{cos}(2\pi x_{5}))/5}+20+e^{1}$
Branin |  $(x_{2}-5.1/(4\pi^{2})x_{1}^{2}+5x_{1}/\pi-6)^{2}+10(1-1/(8\pi))\text{cos}(x_{1})+10$
Currin exponential |  $(1-e^{-1/(2x_{2})})(2300x_{1}^{3}+1900x_{1}^{2}+2092x_{1}+60)/(100x_{1}^{3}+500x_{1}^{2}+4x_{1}+20)$
Gramacy |  $\text{sin}(10\pi x)/(2x)+(x-1)^{4}$
Higdon |  $\text{sin}(2\pi x/10)+0.2\text{sin}(2\pi x/2.5)$