# Sparse within Sparse Gaussian Processes using Neighbor Information

Gia-Lac Tran  Dimitrios Milios  Pietro Michiardi  Maurizio Filippone

###### Abstract

Approximations to Gaussian processes (GPs) based on inducing variables,
combined with variational inference techniques, enable state-of-the-art sparse
approaches to infer GPs at scale through mini-batch based learning. In this
work, we further push the limits of scalability of sparse GPs by allowing
large number of inducing variables without imposing a special structure on the
inducing inputs. In particular, we introduce a novel hierarchical prior, which
imposes sparsity on the set of inducing variables. We treat our model
variationally, and we experimentally show considerable computational gains
compared to standard sparse GPs when sparsity on the inducing variables is
realized considering the nearest inducing inputs of a random mini-batch of the
data. We perform an extensive experimental validation that demonstrates the
effectiveness of our approach compared to the state-of-the-art. Our approach
enables the possibility to use sparse GPs using a large number of inducing
points without incurring a prohibitive computational cost.

Machine Learning, ICML

##  1 Introduction

Gaussian Processes (gps) (Rasmussen & Williams, 2006) offer a powerful
framework to perform inference over functions; being Bayesian, gps provide
rigorous uncertainty quantification and prevent overfitting. However, the
applicability of gps on large datasets is hindered by their computational
complexity of $\mathcal{O}\left(N^{3}\right)$, where $N$ is the training size.
This issue has fuelled a considerable amount of research towards scalable gp
methodologies that operate on a set of _inducing variables_ (Quiñonero Candela
& Rasmussen, 2005). In the literature, there is a plethora of approaches that
offer different treatments of the inducing variables (Lawrence et al., 2002;
Seeger et al., 2003; Snelson & Ghahramani, 2005; Naish-Guzman & Holden, 2007;
Titsias, 2009; Hensman et al., 2013; Wilson & Nickisch, 2015; Hensman et al.,
2015). Some of the more recent approaches, such as Scalable Variational
Gaussian Processes (svgps) (Hensman et al., 2015), allow for the application
of gps to problems with millions of data-points. In most applications of
scalable gps, these are approximated using $M$ inducing points, which results
in a complexity of $\mathcal{O}\left(M^{3}\right)$. It has been shown recently
by Burt et al. (2019) that it is possible to obtain an arbitrarily good
approximation for a certain class of gp models (i.e. conjugate likelihoods,
concentrated distribution for the training data) with $M$ growing more slowly
than $N$. However, the general case remains elusive and it is still possible
that the required value for $M$ may exceed a certain computational budget. Our
result contributes to strengthen our belief that sparsity does not only enjoy
desirable theoretical properties, but it also constitutes an extremely
computationally efficient method in practice.

In this work, we push the limits of scalability and effectiveness of sparse
gps enabling a further reduction in complexity, which can be translated to
higher accuracy by considering a larger set of inducing variables. The idea is
to operate on a subset of $H$ inducing points during training and prediction,
with $H\ll M$, while maintaining a sparse approximation with $M$ inducing
variables. We formalize our strategy by imposing a sparsity-inducing structure
on the prior over the inducing variables and by carrying out a variational
formulation of this model. This extends the original svgp framework and
enables mini-batch based optimization for the variational objective. We then
consider ways to select the set of $H$ inducing points based on neighbor
information; at training time, for a given mini-batch, we activate $H$ out of
$M$ inducing variables considering the nearest inducing inputs to the samples
in the mini-batch, whereas at test time we select inducing variables
corresponding to the inducing inputs which are nearest to the test data-
points. We name our proposal _Sparse within a Sparse gp_ (swsgp). swsgp is
characterized by a number of attractive features: (i) it improves
significantly the prediction quality using a small number of neighboring
inducing inputs, and (ii) it accelerates the training phase, especially when
the total number of inducing points becomes large. We extensively validate
these properties on a variety of regression and classification tasks. We also
showcase swsgp on a large scale classification problem with $M=100,000$; we
are not aware of other approaches that can handle such a large set of inducing
inputs without imposing some special structure on them (e.g., grid) or without
considering one-dimensional inputs.

Hierarchical priors are often applied in Bayesian modeling to achieve
compression and to improve flexibility (Molchanov et al., 2017; Louizos et
al., 2017). To the best of our knowledge, this work is the first to explore
these ideas as means to sparsify the inducing set in sparse gps.

##  2 Related Work and Background

Sparse gps that operate on inducing inputs have been extensively studied in
the last 20 years (Csató & Opper, 2002; Lawrence et al., 2002; Snelson &
Ghahramani, 2005; Quiñonero Candela & Rasmussen, 2005; Naish-Guzman & Holden,
2007). Many attempts on sparse gps specified inducing inputs by satisfying
certain criteria that produce an informative set of inducing variables (Csató
& Opper, 2002; Lawrence et al., 2002; Seeger et al., 2003). A different
treatment has been proposed by Titsias (2009), which involves formulating the
selection of inducing inputs as optimization of a variational lower bound to
the marginal likelihood. The variational framework was later extended so that
stochastic optimization can be admitted, thus improving scalability for
regression (Hensman et al., 2013) and classification (Hensman et al., 2015),
and to deal with large-dimensional inputs (Panos et al., 2018). All the
aforementioned methodologies share a computational complexity of
$\mathcal{O}\left(M^{3}\right)$. Although there have been some attempts in the
literature to infer the appropriate number of inducing points as well as the
inducing inputs (Pourhabib et al., 2014; Burt et al., 2019), a large number of
inducing variables is desirable to improve posterior approximation. In this
work, we present a methodology that builds on the svgp framework (Hensman et
al., 2015) and reduces its complexity, thus increasing the potential of sparse
gp application on even larger datasets and with a larger set of inducing
variables.

A different approach to scalable gps was introduced by Wilson & Nickisch
(2015), namely Kernel Interpolation for Scalable Structured gps (kiss-gp).
This line of work involves arranging a large number of inducing inputs into a
grid structure; this allows one to scale to very large datasets by means of
fast linear algebra. The applicability of kiss-gp on higher-dimensional
problems has been addressed by Wilson et al. (2015) by means of low-
dimensional projections. A more recent extension allows for a constant-time
variance prediction using Lanczos methods (Pleiss et al., 2018). Our work
takes a different approach by keeping the gp prior intact, and by imposing
sparsity on the set of inducing variables.

The local approximation of gps inspired by the the concept of divide-and-
conquer is also a practical solution to implement scalable gps (Kim et al.,
2005; Urtasun & Darrell, 2008; Datta et al., 2016; Park & Huang, 2016; Park &
Apley, 2017). More recently, Liu & Liu (2019) proposed an amortized
variational inference framework where the nearest training points are selected
for inference. In our work, we use neighbor information in a different way, by
incorporating it in a certain hierarchical structure of the auxiliary
variables through a variational scheme.

###  2.1 Scalable Variational Gaussian Processes

Consider a supervised learning problem with inputs
$\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{N})^{\top}$ associated with
labels $\mathbf{y}=(y_{1},\ldots,y_{N})^{\top}$. Given a set of latent
variables $\mathbf{f}=(f_{1},\ldots,f_{N})^{\top}$, gp models assume that
labels are stochastic realizations based on $\mathbf{f}$ and a likelihood
function $p(\mathbf{y}\mid\mathbf{f})$. In svgps, the set of inducing points
is characterized by inducing inputs
$\mathbf{Z}=(\mathbf{z}_{1},\ldots,\mathbf{z}_{M})^{\top}$ and inducing
variables $\mathbf{u}=(u_{1},\ldots,u_{M})^{\top}$. Regarding $\mathbf{f}$ and
$\mathbf{u}$, we have the following joint prior:

$p(\mathbf{f},\mathbf{u})=\mathcal{N}\left(0,\begin{bmatrix}\mathbf{K}_{\mathbf{X}}&\mathbf{K}_{\mathbf{X},\mathbf{Z}}\\\
\mathbf{K}_{\mathbf{Z},\mathbf{X}}&\mathbf{K}_{\mathbf{Z}}\end{bmatrix}\right)\text{,}$
(1)

where $\mathbf{K}_{\mathbf{X}}$, $\mathbf{K}_{\mathbf{Z}}$ and
$\mathbf{K}_{\mathbf{X},\mathbf{Z}}$ are covariance matrices evaluated at the
inputs indicated by the subscripts. The posterior over inducing variables is
approximated by a variational distribution
$q\left(\mathbf{u}\right)=\mathcal{N}\left(\mathbf{u}\mid\mathbf{m},\mathbf{S}\right)$,
while keeping the exact conditional $p(\mathbf{f}\mid\mathbf{u})$ intact, that
is $q(\mathbf{f},\mathbf{u})=p(\mathbf{f}\mid\mathbf{u})q(\mathbf{u})$. The
variational parameters $\mathbf{m}$ and $\mathbf{S}$, as well as the inputs
$\mathbf{Z}$, are optimized by maximizing a lower bound on the marginal
likelihood $p(\mathbf{y}\mid\mathbf{X})=\int
p(\mathbf{y}\mid\mathbf{f})p(\mathbf{f}\mid\mathbf{X})d\mathbf{f}$. The lower
bound on $\log p\left(\mathbf{y}\mid\mathbf{X}\right)$ can be obtained by
considering the form of $q(\mathbf{f},\mathbf{u})$ above and by applying
Jensen’s inequality:

$\begin{split}\mathrm{E}_{q\left(\mathbf{f}\right)}\left[\log
p\left(\mathbf{y}\mid\mathbf{f}\right)\right]-\mathrm{KL}\left(q\left(\mathbf{u}\right)\|\,p\left(\mathbf{u}\right)\right)\text{.}\end{split}$
(2)

The approximate posterior $q\left(\mathbf{f}\right)$ can be computed by
integrating out $\mathbf{u}$: $q\left(\mathbf{f}\right)=\int
q\left(\mathbf{u}\right)p\left(\mathbf{f}\left|\mathbf{u}\right.\right)d\mathbf{u}.$
Thanks to the Gaussian form of $q\left(\mathbf{u}\right)$,
$q\left(\mathbf{f}\right)$ can be computed analytically:

$q\left(\mathbf{f}\right)=\mathcal{N}\left(\mathbf{f}\ |\
\mathbf{A}\mathbf{m},\
\mathbf{K}_{\mathbf{X}}+\mathbf{A}\left(\mathbf{S}-\mathbf{K}_{\mathbf{Z}}\right)\mathbf{A}\right)\text{,}$
(3)

where
$\mathbf{A}=\mathbf{K}_{\mathbf{X},\mathbf{Z}}\mathbf{K}_{\mathbf{Z}}^{-1}$.
When the likelihood factorizes over training points, the lower bound can be
re-written as:

$\textstyle\sum_{i=1}^{N}\mathrm{E}_{q\left(f_{i}\right)}\left[\log
p\left(y_{i}\ |\
f_{i}\right)\right]-\mathrm{KL}\left(q\left(\mathbf{u}\right)\|\,p\left(\mathbf{u}\right)\right)\text{.}$
(4)

Each term of the one-dimensional expectation of the log-likelihood can be
computed by Gauss-Hermite quadrature for any likelihoods (and analytically for
the Gaussian likelihood). The $\mathrm{KL}\left(q\left(\mathbf{u}\right)\ \|\
p\left(\mathbf{u}\right)\right)$ term can be computed analytically given that
$q\left(\mathbf{u}\right)$ and $p\left(\mathbf{u}\right)$ are both Gaussian.
To maintain positive-definiteness of $\mathbf{S}$ and perform unconstrained
optimization, $\mathbf{S}$ is parametrized as
$\mathbf{S}=\mathbf{L}\mathbf{L}^{T}$, with $\mathbf{L}$ lower triangular.

##  3 Sparse Within Sparse GP

We present a novel formulation of sparse gps, which permits the use of a
random subset of the inducing points with little loss in performance. We
introduce a set of binary random variables $\mathbf{w}\in\\{0,1\\}^{M}$ to
govern the inclusion of inducing inputs $\mathbf{Z}$ and the corresponding
variables $\mathbf{u}$. We then employ these random variables to define a
hierarchical structure on the prior as follows:

$p\left(\mathbf{u}\mid\mathbf{w}\right)=\mathcal{N}\left(\mathbf{0},\mathbf{D}_{\mathbf{w}}\mathbf{K}_{\mathbf{Z}}\mathbf{D}_{\mathbf{w}}\right),$
(5)

where $\mathbf{D}_{\mathbf{w}}=\mathrm{diag}\left(\mathbf{w}\right)\text{, and
}\mathbf{w}\sim p\left(\mathbf{w}\right)$. Although the marginalized prior
$p(\mathbf{u})$ is not Gaussian, it is possible to use the joint
$p(\mathbf{u},\mathbf{w})=p(\mathbf{u}\mid\mathbf{w})\ p(\mathbf{w})$ within a
variational scheme. We thus consider a random subset of the inducing points
during the evaluation of the prior in the variational scheme that follows; no
inducing points are permanently removed. Regarding $p(\mathbf{w})$, we
consider an implicit distribution: its analytical form is unknown, but we can
draw samples from it. Later, we will consider $p(\mathbf{w})$ based on the
nearest inducing inputs to random mini-batches of data.

#### Remarks on the prior over $\mathbf{f}$

Our strategy simply assumes a certain structure on the auxiliary variables,
but it has no effect on the prior over $\mathbf{f}$; the latter remains
unchanged. Let $\mathcal{I}$ and $\mathcal{J}$ bet the sets of indices such
that $\mathbf{w}_{\mathcal{I}}=\mathbf{1}$ and
$\mathbf{w}_{\mathcal{J}}=\mathbf{0}$. Given an appropriate ordering, the
conditional $\mathbf{u}\mid\mathbf{w}$ is effectively the element-wise product
$[\mathbf{u}_{\mathcal{I}},\mathbf{u}_{\mathcal{J}}]^{\top}=\mathbf{u}\circ\mathbf{w}$.
This reduces the variances and covariances of some elements of $\mathbf{u}$ to
zero yielding a distribution of this form:

$p\left(\mathbf{f},\mathbf{u}\mid\mathbf{w}\right)=\mathcal{N}\left(\begin{bmatrix}\mathbf{0}\\\
\mathbf{0}\\\
\mathbf{0}\end{bmatrix},\begin{bmatrix}\mathbf{K}_{\mathbf{X}}&\mathbf{K}_{\mathbf{X},\mathbf{Z}_{\mathcal{I}}}&\mathbf{0}\\\
\mathbf{K}_{\mathbf{Z}_{\mathcal{I}},\mathbf{X}}&\mathbf{K}_{\mathbf{Z}_{\mathcal{I}}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}\end{bmatrix}\right)$ (6)

The rows and columns of $\mathbf{u}_{\mathcal{J}}$ can simply be ignored.
Regardless of the value of $\mathbf{w}$, the conditional
$\mathbf{f},\mathbf{u}_{\mathcal{I}}\mid\mathbf{w}$ is always a Gaussian
marginal, as it is a subset of Gaussian variables. The marginalized
$p(\mathbf{f},\mathbf{u})=\int
p(\mathbf{f},\mathbf{u}\mid\mathbf{w})\,p(\mathbf{w})d\mathbf{w}$ is mixture
of Gaussian densities, where the marginal over $\mathbf{f}$ is the same for
every component of the mixture.

Figure 1: The choice of inducing points does not affect the prior samples
drawn from $p(\mathbf{f})$. Left: visualizations of
$\mathbf{f}\mid\mathbf{u},\mathbf{w}$ for different samples of $\mathbf{w}$.
Right: comparison of the marginalized (w.r.t. $\mathbf{u},\mathbf{w}$) prior
over $\mathbf{f}$, against the true $p(\mathbf{f})$.

The effect on $\mathbf{f}$ is demonstrated in Figure 1, where we sample from
the (non-Gaussian) marginalized prior $p(\mathbf{u})$ in two steps: first we
consider an arbitrary random subset $\mathbf{u}_{\mathcal{I}}$, and then we
sample from $p(\mathbf{u}_{\mathcal{I}})\equiv p(\mathbf{u}\mid\mathbf{w})$.
Finally, $\mathbf{f}$ samples are drawn from
$p(\mathbf{f}\mid\mathbf{u}_{\mathcal{I}})$, which only involves the selected
inducing variables $\mathbf{u}_{\mathcal{I}}$. Following Eq. (1), the
conditional $\mathbf{f}\mid\mathbf{u}$ is normally-distributed with mean
$\mathbf{m}_{\mathbf{f}\mid\mathbf{u}_{\mathcal{I}}}=\mathbf{K}_{\mathbf{X},\mathbf{Z}_{\mathcal{I}}}\mathbf{K}_{\mathbf{Z}_{\mathcal{I}}}^{-1}\mathbf{u}_{\mathcal{I}}$
and covariance
$\mathbf{S}_{\mathbf{f}\mid\mathbf{u}_{\mathcal{I}}}=\mathbf{K}_{\mathbf{X}}-\mathbf{K}_{\mathbf{X},\mathbf{Z}_{\mathcal{I}}}\mathbf{K}_{\mathbf{Z}_{\mathcal{I}}}^{-1}\mathbf{K}_{\mathbf{Z}_{\mathcal{I}},\mathbf{X}}$.
These conditionals can be seen for different samples of
$\mathbf{u},\mathbf{w}$ in the left side of Figure 1, while in the right side
we compare the marginalized prior over $\mathbf{f}$ against the true gp prior.

Of course, although the prior remains unchanged, that is not the case for the
posterior approximation. It is well known that the choice of inducing inputs
has an effect on the variational posterior (Titsias, 2009; Burt et al., 2019).
Our choice to impose a hierarchical structure to the inducing variables
through $\mathbf{w}$ effectively changes the model compared to svgp, and we
adapt the variational scheme accordingly.

###  3.1 Lower Bound on Marginal Likelihood

By introducing $\mathbf{u},\mathbf{w}$ and using Jensen’s inequality, the
lower bound on $\log p\left(\mathbf{y}\right)$ can be obtained as follows

$\begin{split}\mathrm{E}_{q\left(\mathbf{u},\mathbf{w}\right)}\log
p\left(\mathbf{y}\left|\mathbf{u},\mathbf{w}\right.\right)-\mathrm{KL}\left(q\left(\mathbf{u},\mathbf{w}\right)\|p\left(\mathbf{u},\mathbf{w}\right)\right),\end{split}$
(7)

where we choose the variational distribution $q$ to reflect the hierarchical
structure of the prior, i.e.
$q\left(\mathbf{u},\mathbf{w}\right)=q(\mathbf{u}\mid\mathbf{w})\
p(\mathbf{w})$. This choice enforces sparsity over the approximate posterior
$q$; the variational parameters are shared among the conditionals
$q(\mathbf{u}\mid\mathbf{w})$, for which we assume:

$q(\mathbf{u}\mid\mathbf{w})=\mathcal{N}\left(\mathbf{u}\left|\mathbf{D}_{\mathbf{w}}\mathbf{m},\mathbf{D}_{\mathbf{w}}\mathbf{S}\mathbf{D}_{\mathbf{w}}\right.\right)$
(8)

By maximizing the variational bound, we aim to obtain a $q$ that performs well
under a sparsified inducing set. We continue by applying Jensen’s inequality
on $p\left(\mathbf{y}\left|\mathbf{u},\mathbf{w}\right.\right)$, obtaining:

$\log
p\left(\mathbf{y}\left|\mathbf{u},\mathbf{w}\right.\right)\geq\mathrm{E}_{p\left(\mathbf{f}|\mathbf{u},\mathbf{w}\right)}\log
p\left(\mathbf{y}\left|\mathbf{f}\right.\right)$ (9)

The bound we describe in Equations (7) and (9) is the same as in svgp, but
with a different variational distribution. In our case, the variational
distribution imposes sparsity for $\mathbf{u}$ by means of $\mathbf{w}$. We
can now substitute (9) into (7), obtaining a bound where we expand
$q(\mathbf{u},\mathbf{w})$ as
$q\left(\mathbf{u}\left|\mathbf{w}\right.\right)p\left(\mathbf{w}\right)$. By
making this assumption, we get the following evidence lower bound
$\mathcal{L}_{\textsc{elbo}}$:

$\begin{split}\sum_{n=1}^{N}\mathrm{E}_{p\left(\mathbf{w}\right)}\bigg{[}&\mathrm{E}_{q\left(\mathbf{u}\left|\mathbf{w}\right.\right)}\mathrm{E}_{p\left(f_{n}\left|\mathbf{u},\mathbf{w}\right.\right)}\log
p\left(y_{n}\left|f_{n}\right.\right)\\\
&-\frac{1}{N}\mathrm{KL}\left(q\left(\mathbf{u}\left|\mathbf{w}\right.\right)\bigg{\|}\,p\left(\mathbf{u}\left|\mathbf{w}\right.\right)\right)\bigg{]}\end{split}$
(10)

Recall that $p\left(\mathbf{w}\right)$ is implicit: although we do not make
any particular assumptions about its analytical form, we can draw samples from
it. Using MC sampling from $p\left(\mathbf{w}\right)$, we can obtain the
approximation $\tilde{\mathcal{L}}_{\textsc{elbo}}$:

$\begin{split}\sum_{n=1}^{N}\bigg{[}&\mathrm{E}_{q\left(\mathbf{u}\left|\tilde{\mathbf{w}}^{(n)}\right.\right)}\mathrm{E}_{p\left(f_{n}\left|\mathbf{u},\tilde{\mathbf{w}}^{(n)}\right.\right)}\log
p\left(y_{n}\left|f_{n}\right.\right)\\\
&-\frac{1}{N}\mathrm{KL}\left(q\left(\mathbf{u}\left|\tilde{\mathbf{w}}^{(n)}\right.\right)\bigg{\|}\,p\left(\mathbf{u}\left|\tilde{\mathbf{w}}^{(n)}\right.\right)\right)\bigg{]},\end{split}$
(11)

where $\tilde{\mathbf{w}}^{(n)}$ is sampled from $p\left(\mathbf{w}\right)$.

#### Sampling from the set of inducing points.

Recall that any sample $\tilde{\mathbf{w}}$ from $p\left(\mathbf{w}\right)$ is
a binary vector, i.e. $\mathbf{w}\in\\{0,1\\}^{M}$. In case all elements of
$\mathbf{w}$ are set to one, our approach recovers the original svgp with
computational cost of $\mathcal{O}\left(M^{3}\right)$ coming from computing
$p\left(f_{n}\left|\mathbf{u},\tilde{\mathbf{w}}=\mathbf{1}\right.\right)$ and
$\mathrm{KL}\left(q\left(\mathbf{u}\left|\mathbf{w}\right.\right)\|p\left(\mathbf{u}\left|\mathbf{w}\right.\right)\right)$
in the elbo. When any $\tilde{w}_{i}$ is set to zero, the entries of the
$i$-th row and $i$-th column of the covariance matrix in
$p\left(\mathbf{u}\left|\mathbf{w}\right.\right)$ and
$q\left(\mathbf{u}\left|\mathbf{w}\right.\right)$ are zero. This means that
the $i$-th variable becomes unnecessary, so we get rid of $i$-th row and
column in these matrices, and also eliminate the $i$-th element in mean
vectors of $q\left(\mathbf{u}\left|\mathbf{w}\right.\right)$ and
$p\left(\mathbf{u}\left|\mathbf{w}\right.\right)$. This is equivalent to
selecting a set of active inducing points in each training iteration.

###  3.2 H-nearest Inducing Inputs

Despite the fact that $p(\mathbf{w})$ is an implicit distribution, we have
been able to define and calculate a variational bound, assuming we can sample
from $p(\mathbf{w})$. We shall now describe our sampling strategy, which
relies on neighbor information of random mini-batches.

We introduce $\mathbf{Z}^{H}_{\mathbf{x}}$ as the set of $H$-nearest inducing
inputs. Intuitively, the prediction for an unseen data $\mathbf{x}$ using
$\mathbf{Z}^{H}_{\mathbf{x}}$ is a good approximation of the prediction using
all $M$ inducing points, that is $\mathbf{Z}^{M}_{\mathbf{x}}$. This can be
verified by looking at the predictive mean, which is expressed as a linear
combination of kernel functions evaluated between training points and a test
point, as in Eq. (3). The majority of the contribution is given by the
inducing points with the largest kernel values, so we can use this as a
criterion to establish whether an inducing input is “close” to an input vector
(the effect of different kernels on the definition of nearest neighbors is
explored in the supplement). With this intuition, $p\left(\mathbf{w}\right)$
becomes a deterministic function $w\left(\mathbf{x}\right)$ indicating which
inducing inputs are activated. For mini-batch based training, the value of
$\mathbf{w}$ remains random, as it depends on the elements $\mathbf{x}$ that
are selected in the random mini-batch; this materializes the sampling from the
implicit distribution $p(\mathbf{w})$. The maximization of the elbo in the
setting described is summarized in Algorithm 1 (swsgp).

Algorithm 1 Sparse within sparse gp (swsgp)

0: $\mathcal{D}$, $H$, $M$.

0: The optimum of trainable parameters $\boldsymbol{\theta}$.

1: Initialize $\boldsymbol{\theta}$ including kernel’s parameters,
$\mathbf{Z}$, $\mathbf{m}$ and $\mathbf{L}$ which can be used to construct
$\mathbf{S}$, i.e. $\mathbf{S}=\mathbf{L}\mathbf{L}^{T}$.

2: while stopping criteria is False do

3: $\mathrm{ELL}\leftarrow 0$ and $\mathrm{KL}\leftarrow 0$.

4: Sample mini-batch $\mathrm{B}$ of size $n_{\mathrm{B}}$ from $\mathcal{D}$.

5: for $\left(\mathbf{x}_{i},y_{i}\right)\in\mathrm{B}$ do

6: Find $\mathbf{Z}^{H}_{\mathbf{x}_{i}}$, i.e. the $H$-nearest $\mathbf{Z}$
to $\mathbf{x}_{i}$.

7: Compute $w\left(\mathbf{x}_{i}\right)$ using
$\mathbf{Z}^{H}_{\mathbf{x}_{i}}$ as in (12).

8: Extract $\mathbf{m}_{w\left(\mathbf{x}_{i}\right)}$ and
$\mathbf{S}_{w\left(\mathbf{x}_{i}\right)}$ from $\mathbf{m}$ and
$\mathbf{L}$.

9: Compute $q\left(f_{i}\left|w\left(\mathbf{x}_{i}\right)\right.\right)$ as
in (13).

10:
$\mathrm{ELL}\leftarrow\mathrm{ELL}+\mathrm{E}_{q\left(f_{i}\left|w\left(\mathbf{x}_{i}\right)\right.\right)}\log
p\left(y_{i}\left|f_{i}\right.\right)$.

11:
$\mathrm{KL}\leftarrow\mathrm{KL}+\mathrm{KL}\left(q\left(\mathbf{u}_{w\left(\mathbf{x}_{i}\right)}\right)\|p\left(\mathbf{u}_{w\left(\mathbf{x}_{i}\right)}\right)\right)$

12: end for

13:
$\tilde{\mathcal{L}}_{\textsc{elbo}}\leftarrow\frac{N}{n_{\mathrm{B}}}\mathrm{ELL}-\frac{1}{n_{\mathrm{B}}}\mathrm{KL}$.

14: Update $\boldsymbol{\theta}$ using the derivative of
$\tilde{\mathcal{L}}_{\textsc{elbo}}$.

15: end while

#### Predictive distribution

Contrary to what happens at training time, where mini-batches of data are
drawn randomly, at test time the inputs of interest are not random; we need to
describe the predictive distribution in terms of the deterministic function
$w\left(\mathbf{x}\right)$. In fact, if we would like to approximate the
predictive distribution at $\mathbf{x}_{n}$ using $H$-nearest inducing inputs
to $\mathbf{x}$, i.e. $\mathbf{Z}^{H}_{\mathbf{x}_{n}}$, then
$w\left(\mathbf{x}\right)=\left[w^{(1)}_{\mathbf{x}}...w^{(M)}_{\mathbf{x}}\right]^{T}$
where,

$w^{(m)}_{\mathbf{x}}=\begin{cases}1&\mbox{if
}\mathbf{z}_{m}\in\mathbf{Z}^{H}_{\mathbf{x}}\\\
0&\mbox{else}\end{cases}\text{ , with }m=1,...,M$ (12)

We extract the relevant elements using $w\left(\mathbf{x}\right)$; for the
mean, we have
$\mathbf{m}_{w\left(x_{i}\right)}=\mathbf{D}_{w\left(x_{i}\right)}\mathbf{m}$,
and for the covariance we select the appropriate rows and columns using
$\mathbf{S}_{w\left(\mathbf{x}_{i}\right)}=\mathbf{D}_{w\left(\mathbf{x}_{i}\right)}\mathbf{S}\mathbf{D}_{w\left(\mathbf{x}_{i}\right)}$.
The approximate posterior over $f_{i}$ given $w\left(\mathbf{x}_{i}\right)$,
i.e. $q\left(f_{i}\left|w\left(\mathbf{x}_{i}\right)\right.\right)$ is:

$\begin{split}\mathcal{N}\bigg{(}f_{i}\mid&\mathbf{A}_{\mathbf{x}_{i}}\mathbf{m}_{w\left(\mathbf{x}_{i}\right)},\\\
&\mathbf{K}_{\mathbf{x}_{i}}+\mathbf{A}_{\mathbf{x}_{i}}\left(\mathbf{S}_{w\left(\mathbf{x}_{i}\right)}-\mathbf{K}_{\mathbf{Z}^{H}_{\mathbf{x}_{i}}}\right)\mathbf{A}_{\mathbf{x}_{i}}^{\top}\bigg{)}\text{,}\end{split}$
(13)

where
$\mathbf{A}_{\mathbf{x}_{i}}=\mathbf{K}_{\mathbf{x}_{i},\mathbf{Z}_{\mathbf{x}_{i}}^{H}}\mathbf{K}_{\mathbf{Z}_{\mathbf{x}_{i}}^{H}}^{-1}$.

#### Limitations

What we presented allows one to compute predictive distributions for
individual test points following Eq. (13). As a result, it is not possible to
obtain a full covariance across predictions over multiple test points. While
most performance and uncertainty quantification metrics do not consider
covariances, this might be a limitation in applications where covariances are
essential. One way to work around this is to consider the union of nearest
neighbors at test time; however, this would be inconsistent with Eq. (13) and
the training procedure should be modified accordingly considering the union of
nearest neighbors for training points within mini-batches. The considerations
suggest a further limitation of swsgp, that is that the predictive
distribution becomes dependent on the number of test points that are consdered
within a batch at test time. Nevertheless, the model is built and trained to
deal with sparsity at test time, and this kind of inconsistency appears only
when testing the model in conditions which are different compared to those at
training time. Besides, the experiments show that swsgp performs extremely
well against various state-of-the-art competitors on a wide range of
experiments.

One-dimensional regression example. We visualize the posterior distribution
for a synthetic dataset generated on a one-dimensional input space. We execute
svgp and swsgp, and depict the posterior distributions of these two methods by
showing the predictive means (orange lines) and the 95% credible intervals
(shaded areas) in Figure 2. We consider identical settings for the two methods
(i.e. $128$ inducing points, kernel parameters, likelihood variance) and a
neighbor area of $16$ for swsgp; a full account of the setup can be found in
the supplement. We see that although the models are different, the predictive
distributions appear remarkably similar. A more extensive evaluation follows
in Section 4.

Figure 2: Visualization of posterior distribution of svgp and swsgp. In both
cases, we consider 128 inducing points; in terms of our scheme (swsgp) we use
16 neighbors.

###  3.3 Complexity

The computational cost of swsgp is dominated by lines 6, 8 and 9 in Algorithm
1. For each data-point $\left(\mathbf{x}_{i},y_{i}\right)$ in mini-batch
$\mathrm{B}$, we need to find the $H$ nearest inducing neighbors
$\mathbf{Z}^{H}_{\mathbf{x}_{i}}$ for $n_{\mathrm{B}}$ points in line 6, where
$n_{\mathrm{B}}=|\mathrm{B}|$; this contributes to the worst-case complexity
by $\mathcal{O}\left(n_{\mathrm{B}}MH\right)$.

In line 8, we extract relevant parameters from $\mathbf{m}$ and $\mathbf{L}$.
We focus on the cost of extracting $\mathbf{S}_{w\left(\mathbf{x}_{i}\right)}$
from $\mathbf{L}$. Similar to svgp (Section 2.1), we consider
$\mathbf{S}=\mathbf{L}\mathbf{L}^{T}$, where $\mathbf{L}$ is lower triangular.
We extract
$\mathbf{L}_{w\left(\mathbf{x}_{i}\right)}=\mathbf{D}_{w\left(\mathbf{x}_{i}\right)}\mathbf{L}$
which contains the rows of $\mathbf{L}$. Then, we compute
$\mathbf{S}_{w\left(\mathbf{x}_{i}\right)}$ by
$\mathbf{S}_{w\left(\mathbf{x}_{i}\right)}=\mathbf{L}_{w\left(\mathbf{x}_{i}\right)}\mathbf{L}_{w\left(\mathbf{x}_{i}\right)}^{T}$.
The computational complexity of selecting the variational parameters is
$\mathcal{O}\left(n_{\mathrm{B}}MH^{2}\right)$.

Finally, the computation of approximating the predictive distribution in line
9 requires $\mathcal{O}\left(n_{\mathrm{B}}H^{3}\right)$. The overall
complexity for swsgp in the general case is
$\mathcal{O}\left(n_{\mathrm{B}}MH+n_{\mathrm{B}}MH^{2}+n_{\mathrm{B}}H^{3}\right)$,
which is a significant improvement over the $\mathcal{O}\left(M^{3}\right)$
complexity of standard svgp, assuming that $n_{\mathrm{B}},H\ll M$. If we
choose $\mathbf{S}$ to be diagonal, the total complexity reduces to
$\mathcal{O}\left(n_{\mathrm{B}}MH+n_{\mathrm{B}}H^{3}\right)$; if we
additionally consider $\mathbf{Z}$ to be fixed, the computational cost is
$O(n_{\mathrm{B}}H^{3})$. In the experiments of Section 4 we explore all these
settings.

##  4 Experiments

We conduct experiments to evaluate swsgp on a variety of experimental
conditions. Our approach is denoted by swsgp-m-h, where $M$ inducing points
are used and $H$ determines how many neighbors are selected. The locations of
the inducing points are optimized, unless stated otherwise. We introduce
svgp-m and svgp-m-h as competitors; svgp-m uses $M$ inducing points. svgp-m-h,
instead, refers to svgp using $M$ inducing points at training time and
$H$-nearest inducing inputs at test time. We also use svgp-16m using $16$
inducing points as a reference. The comparison is carried out on some UCI data
sets for regression and classification, i.e., powerplant, kin, protein, eeg,
credit and spam. We also consider larger scale data sets, such as, wave, query
and the airline data or images classification on mnist. The task for the
airline data set is the classification of whether flights are to subject to a
delay, and we follow the same setup as in Hensman et al. (2013) and Wilson et
al. (2016). Regarding mnist, we consider the binary classification setup of
separating odd and even digits (pixels normalized between $0$ and $1$), and we
use a Bernoulli likelihood with probit inverse link function. We use the
Matérn-$\nicefrac{{5}}{{2}}$ kernel in all cases except for the airline
dataset, where the sum of a Matérn-$\nicefrac{{3}}{{2}}$ and a linear kernel
is used, similar to Hensman et al. (2015). All models are trained using the
Adam optimizer (Kingma & Ba, 2015) with a learning rate of $0.001$ and a mini-
batch size of $64$. The likelihood for regression and binary classification
are set to Gaussian and probit function, respectively. In regression tasks, we
report the test root mean squared error (rmse) and the test mean negative log-
likelihood (mnll), whereas we report the test error rate (err) and mnll in
classification tasks. The results are averaged over $10$ folds. All models in
section 4.1 are trained over $300,000$ iterations.

###  4.1 Impact of M and H

We begin our evaluation by investigating the behavior of swsgp-m-h with
respect to $M$ and $H$. Figure 3 shows that swsgp-m-h consistently outperforms
svgp-m-h. This suggests that including neighbor information at prediction
time, combined with the use of a larger set of inducing points alone is not
enough to obtain competitive performance, and that only thanks to the
sparsity-inducing prior over latent variables, this yields improvements.
Crucially, the performance metrics obtained by swsgp are comparable with those
obtained by svgp-m, while at each iteration only a small subset of $H$ out of
$M$ inducing points are updated, carrying a significant complexity reduction.

Figure 3:  Evaluation of swsgp on uci data sets with various configuration for
$M$ and $H$. The title of each sub-figure follows the format of [name of data
set]-[data dimensions]. The black up-triangles are for svgp with $M$ inducing
points. The cyan down-triangles are for svgp with $16$ inducing points. The
red circles are for svgp training with $M$ inducing points and the prediction
at an unseen data $\mathbf{x}$ are made by $\mathbf{Z}_{\mathbf{x}}^{H}$. The
green squares, stars and plus are for swsgp with $H$ of $4$, $8$ and $16$
respectively. In these experiments, $M$ varies from $32$ to $128$, as shown on
horizontal axes. The standard deviation of the error metrics over 10 folds is
represented by vertical bars; they are very small for most configurations.

###  4.2 Running Time

Table 1: Comparison of running time between svgp and swsgp. The running time
for a training iteration is denoted by $t_{1}$. The testing time for an unseen
example is denoted by $t_{2}$. Times are in milliseconds. The corresponding
evaluated metrics ,i.e. err and mnll are also shown in terms of training time.
These figures are averaged over 5 folds. The models for eeg and mnist are
trained in 90 minutes and 6 hours respectively.

(a) Running times for svgp-256m and swsgp-256m-4h on eeg. Configuration |  $t_{1}$(ms) |  $t_{2}$(ms) |  err |  mnll
---|---|---|---|---
svgp-256m |  $\mathbf{21.42}$ |  1.43 |  $\mathbf{0.17}$ |  $\mathbf{0.34}$
swsgp-256m-4h |  26.18 |  $\mathbf{0.56}$ |  0.18 |  0.36

(b) Running times for svgp-1024m and swsgp-1024m-4h on mnist. Configuration |  $t_{1}$(ms) |  $t_{2}$(ms) |  err |  mnll
---|---|---|---|---
svgp-1024m |  516 |  21.6 |  0.02 |  0.066
swsgp-1024m-4h |  $\mathbf{233}$ |  $\mathbf{1.77}$ |  $\mathbf{0.016}$ |  $\mathbf{0.05}$

(c) err and mnll over training time (in minutes) on eeg.

(d) err and mnll over training time (in minutes) on mnist.

In Table 1, we report the training and testing times of swsgp and svgp In
svgp, we set $M=256$ for eeg, and $1024$ for mnist, i.e. svgp-256m and
svgp-1024m. In our approach, we use the same $M$ and we set $H$ to $4$, i.e.
swsgp-256m-4h and swsgp-1024m-4h. In Table 1, we stress that $t_{1}$ and
$t_{2}$ in swsgp take into account the computation of finding neighbors
inducing inputs for each data-point. In svgp, we assume that
$\mathbf{K}_{\mathbf{Z}}^{-1}$ is pre-computed and saved after the training
phase. Therefore, the computational cost to evaluate the predictive
distribution on a single test point is $\mathcal{O}\left(M^{2}\right)$. The
time $t_{2}$ in svgp refers to the execution time of carrying out predictions
with the complexity of $\mathcal{O}\left(M^{2}\right)$.

The results in Tab 1 show a consistent improvement at test time compared to
svgp across all values of $H$ and $M$. At training time, the results show a
trend dependent on the number $M$ of inducing points. Not surprisingly, swsgp
offers limited improvements when $M$ is small. Considering eeg in which $M$ is
set to $256$, svgp is faster than swsgp in terms of training time. This is
because the inversion of a $256\times 256$ matrix requires less time than
finding the neighbors and inverting several $4\times 4$ matrices. However, Tab
1 shows dramatic speedups compared to svgp when the number of inducing points
$M$ is large. When $M=1024$ on mnist, swsgp-1024m-4h is faster than svgp-1024m
at training time. This is due to the inversion of the $1024\times 1024$ kernel
matrix being a burden for svgp, whereas swsgp deals with much cheaper
computations. In addition, we show the corresponding err and mnll of each
model when we train svgp and swsgp on eeg and mnist. On the eeg data set, our
method is comparable with svgp. On mnist, swsgp reaches high accuracies
significantly faster compared to svgp.

###  4.3 Large-scale Sparse GP Modeling with a Huge Number of Inducing Points

Here we show that swsgp allows one to use sparse gps with a massive number of
inducing points without incurring a prohibitive computational cost. We employ
several large-scale data sets, i.e. protein, eeg, wave, query and airline with
5 millions training samples. We test swsgp with a large number of inducing
points, ranging from $3,000$ to $100,000$. In this case, we keep the inducing
locations fixed for swsgp. We have attempted to run svgp with such large
values of $M$ without success (out of memory in a system with 32GB of RAM).
Therefore, as a baseline we execute svgp-128m and svgp-512m and report the
results of svgp with the configuration in Hensman et al. (2015).

In swsgp, we impose a diagonal matrix $\mathbf{S}$ in the variational
distribution $q\left(\mathbf{u}\mid\mathbf{w}\right)$, and we fix the position
of the inducing inputs during training. By fixing the inducing inputs, we can
operate with pre-computed information about which inducing inputs are
neighbors of training inputs. Thanks to these settings, swsgp’s training phase
requires $\mathcal{O}\left(n_{\mathrm{B}}H^{3}\right)$ operations only, where
$n_{\mathrm{B}}$ is the mini-batch size. Due to the appropriate choice of $H$
and $n_{\mathrm{B}}$, and the computational cost being independent of $M$,
unlike svgp, we can successfully run swsgp with an unprecedented number of
inducing point, e.g. $M=100,000$.

Figure 4: swsgp with large number of inducing points. The figure shows the
progression of rmse (err) and mnll over time. Horizontal axes indicate the
running times in minute. The title of each sub-figure follows the format of
[name of data set]-[training size]-[testing size]-[data dimensions].

In airline experiments, by setting $H$ and the mini-batch size $n$ to $100$
and $64$ respectively, in about $24$ hours of training we could run
swsgp-100,000m-100h for one million iterations. The err and mnll of
swsgp-100,000m-100h evaluated on the test set are $21\%$ and $0.48$,
respectively, while the err and mnll of svgp-200m published in Hensman et al.
(2015) are about $34\%$ and $0.61$, respectively. To the best of our
knowledge, swsgp is the first to enable sparse gps with such a large set of
inducing points without imposing a grid structure on the inducing inputs. In
addition, in Figure 4 we also show that swsgp using $3,000$ or $5,000$ or
$10,000$ inducing points outperforms svgp using $128$ and $512$ inducing
points on other data sets.

We conclude by reporting comparisons with other gp-based models. In
particular, we compare against the Stochastic Variational Deep Kernel Learning
(svdkl) (Wilson et al., 2016) and the Deep gp approximated with random
features (dgp-rbf) (Cutajar et al., 2017). In the former, kiss-gp is trained
on top of a deep neural network which is optimized during training, and in the
latter the layers of a deep gp are approximated as parametric models using
random feature expansions. Both competitors feature mini-batch based learning,
so this represents a challenging test for swsgp. The results in Table 2 show
that swsgp is comparable with these competitors on various data sets. We
believe that this is a remarkable result obtained by our shallow swsgp,
supporting the conclusions of previous works showing that advances in kernel
methods can result in performance which are competitive with deep learning
approaches (see, e.g., Rudi et al. (2017)).

Table 2: Comparison of swsgp, kiss-gp (Wilson & Nickisch, 2015), svdkl (Wilson
et al., 2016) and dgp-rbf (Cutajar et al., 2017). The results swsgp and kiss-
gp are averaged over 5 folds. In order to deal with the difficulties of kiss-
gp to handle large-dimensional input spaces, we followed (Wilson et al.,
2015), and we linearly projected the inputs to a two-dimensional space using a
linear transformation that is learned at training time, and used a grid of
size 100.

(a) airline Method |  err |  mnll
---|---|---
swsgp-100km-100h |  $\mathbf{0.210\pm 0.012}$ |  $0.48\pm 0.015$
svdkl |  $0.22$ |  $\mathbf{0.46}$
dgp-rbf |  $\mathbf{0.21}$ |  $\mathbf{0.46}$

(b) powerplant Method |  rmse |  mnll
---|---|---
swsgp-64m-4h |  $\mathbf{4.095\pm 0.145}$ |  $\mathbf{2.371\pm 0.038}$
kiss-gp |  $4.459\pm 0.355$ |  $3.074\pm 0.037$

(c) protein Method |  rmse |  mnll
---|---|---
swsgp-64m-4h |  $\mathbf{0.773\pm 0.002}$ |  $\mathbf{0.702\pm 0.003}$
kiss-gp |  $0.816\pm 0.001$ |  $0.904\pm 0.017$

###  4.4 Comparison to Local GPs

We finally demonstrate that swsgp behaves differently from other approaches
that use local approximations of gps. We consider two well-established
approaches of local gps proposed by Kim et al. (2005) and Urtasun & Darrell
(2008). Following Liu et al. (2020), we refer to these methods as _Inductive_
gps and _Transductive_ gps, respectively. We run all methods on two regression
data sets: powerplant and kin. We set the number of local experts to $64$, and
we use the same number of inducing points for swsgp (with $H$ either 4 or 8).
As the size of powerplant and kin are approximately $7,000$, we set the number
of training points governed by a local expert to $100$. For the local gp
approaches, we choose $64$ locations in the input space using the $K$-means
algorithm, and for each location we choose $100$ neighboring points; we then
train the corresponding local gp expert. For the testing phase, inductive gps
simply rely on the nearest local experts to an unseen point $\mathbf{x}_{*}$.
For transductive gps, we use $100$ neighbors of $\mathbf{x}_{*}$ and the
nearest local expert to make predictions. In Table 3, we summarize rmse and
mnll for all methods; swsgp clearly outperforms the local gp approaches in
terms of mnll.

Table 3: Comparison with Local gp approximations. Method |  powerplant |  kin
---|---|---
|  rmse $\mid$ mnll |  rmse $\mid$ mnll
swsgp-64-4 |  $4.27\mid 2.41$ |  $0.11\mid-1.27$
swsgp-64-8 |  $\textbf{4.24}\mid\textbf{2.40}$ |  $0.10\mid\textbf{-1.38}$
Inductive GPs |  $9.93\mid 38.38$ |  $0.13\mid-0.40$
Transductive GPs |  $6.17\mid 18.78$ |  $\textbf{0.09}\mid-0.65$

###  4.5 Joint Predictive Covariances

Figure 5:  swsgp-u on powerplant. The horizontal axis shows the mini-batch
size of the testing set. We evaluate swsgp-u with various numbers of nearest
inducing points from $4$ to $64$, i.e. $4H$, $8H$, $16H$, $32H$, and $64H$.
The green markers show swsgp-u. For example, in order to plot the green
squares $4H$, we firstly find the 4-nearest inducing points of each sample
within a testing batch. Then, we use the union of these nearest inducing
points to make a joint prediction for all samples in the testing batch.

In this section, we propose an extension of swsgp where the training procedure
is modified by considering the union of nearest neighbors for training points
within mini-batches. We refer to this as swsgp-u, where u stands for “union”.
In the experiments, we execute swsgp-u-512m-128m where the total number of
inducing points is $512$, and the number of active inducing points in each
training iteration is fixed to 128. To illustrate the effectiveness of the
modification, we compare swsgp-u-512m-128m against svgp-512m and full gps. Due
to the computational complexity of full gps, we only use $2,000$ training
samples and fix the computational budget to six hours. In the evaluation
phase, we randomly group the test set based on the specified mini-batch size.
Next, the joint predictive distributions of each testing batch are obtained,
which we use to compute rmse and mnll. The results in Figure 5 are averaged by
repeating the process $10$ times. The results indicate that swsgp-u works well
when using an adequate number of nearest inducing points. In the figure, we
also assess the quality of the covariance approximation of swsgp-u and svgp by
reporting the mean square error with respect to the covariance of full gps.
These results show that swsgp-u provides a good approximation of joint
predictive covariances.

##  5 Conclusions

Sparse approaches that rely on inducing points have met with success in
reducing the complexity of gp regression and classification. However, these
methods are limited by the number of inducing inputs that is required to
obtain an accurate approximation of the true gp model. A large number of
inducing inputs is often necessary in cases of very large datasets, which
marks the limits of practical applications for most gp-based approaches.

In this work, we further improve the computational gains of sparse gps. We
proposed swsgp, a novel methodology that imposes a hierarchical and sparsity-
inducing effect on the prior over the inducing variables. This has been
realized as a conditional gp given a random subset of the inducing points,
which is defined as the nearest neighbors of random mini-batches of data. We
have developed an appropriate variational bound which can be estimated in an
unbiased way by means of mini-batches. We have performed an extensive
experimental campaign that demonstrated the superior scalability properties of
swsgp compared to the state-of-the-art.

#### Acknowledgments

MF gratefully acknowledges support from the AXA Research Fund and the Agence
Nationale de la Recherche (grant ANR-18-CE46-0002 and ANR-19-P3IA-0002).

## References

  * Burt et al. (2019) Burt, D., Rasmussen, C. E., and Van Der Wilk, M.  Rates of convergence for sparse variational Gaussian process regression.  In _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pp. 862–871. PMLR, 2019.
  * Csató & Opper (2002) Csató, L. and Opper, M.  Sparse on-line gaussian processes.  _Neural Computation_ , 14(3):641–668, 2002.  ISSN 0899-7667.
  * Cutajar et al. (2017) Cutajar, K., Bonilla, E. V., Michiardi, P., and Filippone, M.  Random feature expansions for deep Gaussian processes.  In _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pp. 884–893. PMLR, 2017.
  * Datta et al. (2016) Datta, A., Banerjee, S., Finley, A., and Gelfand, A.  On nearest-neighbor gaussian process models for massive spatial data: Nearest-neighbor gaussian process models.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 8, 08 2016.
  * Hensman et al. (2013) Hensman, J., Fusi, N., and Lawrence, N. D.  Gaussian processes for big data.  In _Proceedings of the 29th Conference on Uncertainty in Artificial Intelligence_ , pp. 282–290. AUAI Press, 2013.
  * Hensman et al. (2015) Hensman, J., Matthews, A., and Ghahramani, Z.  Scalable Variational Gaussian Process Classification.  In _Proceedings of the 18th International Conference on Artificial Intelligence and Statistics_ , volume 38 of _Proceedings of Machine Learning Research_ , pp. 351–360. PMLR, 2015.
  * Kim et al. (2005) Kim, H.-M., Mallick, B., and Holmes, C.  Analyzing nonstationary spatial data using piecewise gaussian processes.  _Journal of the American Statistical Association_ , 100:653–668, 02 2005.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  In _3rd International Conference on Learning Representations_ , 2015\.
  * Lawrence et al. (2002) Lawrence, N. D., Seeger, M., and Herbrich, R.  Fast Sparse Gaussian Process Methods: The Informative Vector Machine.  In _Advances in Neural Information Processing Systems 15_ , pp. 625–632. MIT Press, 2002.
  * Liu et al. (2020) Liu, H., Ong, Y., Shen, X., and Cai, J.  When gaussian process meets big data: A review of scalable gps.  _IEEE Transactions on Neural Networks and Learning Systems_ , 31:4405–4423, 2020.
  * Liu & Liu (2019) Liu, L. and Liu, L.  Amortized variational inference with graph convolutional networks for gaussian processes.  In _Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics_ , volume 89 of _Proceedings of Machine Learning Research_ , pp. 2291–2300. PMLR, 2019.
  * Louizos et al. (2017) Louizos, C., Ullrich, K., and Welling, M.  Bayesian compression for deep learning.  In _Advances in Neural Information Processing Systems 30_ , pp. 3288–3298. Curran Associates, Inc., 2017.
  * Molchanov et al. (2017) Molchanov, D., Ashukha, A., and Vetrov, D.  Variational dropout sparsifies deep neural networks.  In _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pp. 2498–2507. PMLR, 2017.
  * Naish-Guzman & Holden (2007) Naish-Guzman, A. and Holden, S.  The generalized FITC approximation.  In _Advances in Neural Information Processing Systems 20_ , pp. 1057–1064. Curran Associates Inc., 2007.  ISBN 978-1-60560-352-0.
  * Panos et al. (2018) Panos, A., Dellaportas, P., and Titsias, M. K.  Fully scalable gaussian processes using subspace inducing inputs.  _CoRR_ , abs/1807.02537, 2018.
  * Park & Apley (2017) Park, C. and Apley, D.  Patchwork kriging for large-scale gaussian process regression.  _Journal of Machine Learning Research_ , 19, 01 2017.
  * Park & Huang (2016) Park, C. and Huang, J. Z.  Efficient computation of gaussian process regression for large spatial data sets by patching local gaussian processes.  _J. Mach. Learn. Res._ , 17(1):6071–6099, January 2016.  ISSN 1532-4435.
  * Pleiss et al. (2018) Pleiss, G., Gardner, J., Weinberger, K., and Wilson, A. G.  Constant-time predictive distributions for Gaussian processes.  In _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 4114–4123. PMLR, 2018.
  * Pourhabib et al. (2014) Pourhabib, A., Liang, F., and Ding, Y.  Bayesian site selection for fast Gaussian process regression.  _Institute of Industrial Engineers Transactions_ , 46(5):543–555, 2014.
  * Quiñonero Candela & Rasmussen (2005) Quiñonero Candela, J. and Rasmussen, C. E.  A unifying view of sparse approximate Gaussian process regression.  _Journal of Machine Learning Research_ , 6:1939–1959, 2005\.  ISSN 1532-4435.
  * Rasmussen & Williams (2006) Rasmussen, C. E. and Williams, C.  _Gaussian Processes for Machine Learning_.  MIT Press, 2006.
  * Rudi et al. (2017) Rudi, A., Carratino, L., and Rosasco, L.  FALKON: An Optimal Large Scale Kernel Method.  In _Advances in Neural Information Processing Systems 30_ , pp. 3888–3898. Curran Associates, Inc., 2017.
  * Seeger et al. (2003) Seeger, M., Williams, C. K. I., and Lawrence, N. D.  Fast forward selection to speed up sparse Gaussian process regression.  In _Artificial Intelligence and Statistics 9_ , 2003.
  * Snelson & Ghahramani (2005) Snelson, E. and Ghahramani, Z.  Sparse Gaussian Processes using Pseudo-inputs.  In _Advances in Neural Information Processing Systems 18_ , pp. 1257–1264. MIT Press, 2005.
  * Titsias (2009) Titsias, M. K.  Variational Learning of Inducing Variables in Sparse Gaussian Processes.  In _Proceedings of the 12th International Conference on Artificial Intelligence and Statistics_ , volume 5 of _Proceedings of Machine Learning Research_ , pp. 567–574. PMLR, 2009.
  * Urtasun & Darrell (2008) Urtasun, R. and Darrell, T.  Sparse probabilistic regression for activity-independent human pose inference.  In _2008 IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 1–8, 2008.
  * Wilson & Nickisch (2015) Wilson, A. and Nickisch, H.  Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP).  In _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pp. 1775–1784. PMLR, 2015.
  * Wilson et al. (2015) Wilson, A., Dann, C., and Nickisch, H.  Thoughts on massively scalable gaussian processes.  _ArXiv_ , abs/1511.01870, 2015.
  * Wilson et al. (2016) Wilson, A. G., Hu, Z., Salakhutdinov, R. R., and Xing, E. P.  Stochastic Variational Deep Kernel Learning.  In Lee, D. D., Sugiyama, M., Luxburg, U. V., Guyon, I., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 29_ , pp. 2586–2594. Curran Associates, Inc., 2016.