\subsection{Scalable Variational Gaussian Processes}
\label{subsection: svgp}
Consider a supervised learning problem with inputs $\Xmat = (\xvect_1, \ldots, \xvect_N)^{\top}$ associated with labels $\yvect = ( y_1, \ldots, y_N )^{\top}$.
Given a set of latent variables $\fvect = (f_1, \ldots, f_N)^{\top}$, \name{gp} models assume that labels are stochastic realizations based on $\fvect$ and a likelihood function $p(\yvect \mid \fvect)$.
In \textsc{svgp}s\xspace, the set of inducing points is characterized by inducing inputs $\mathbf{Z} = (\zvect_1, \ldots, \zvect_{M})^{\top}$ and inducing variables $\mathbf{u} = (u_1, \ldots, u_{M})^{\top}$.
Regarding $\fvect$ and $\mathbf{u}$, we have the following joint prior:
\begin{equation}
	p(\fvect, \mathbf{u}) = \norm\left(0,
	\begin{bmatrix}
	\Kmat_\Xmat 	& \Kmat_{\Xmat, \mathbf{Z}} \\
	\Kmat_{\mathbf{Z}, \Xmat} 	& \Kmat_\mathbf{Z}
	\end{bmatrix}
	\right) \text{,}
\label{eq:fu_joint}
\end{equation}
where $\Kmat_\Xmat$, $\Kmat_\mathbf{Z}$ and $\Kmat_{\Xmat, \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)=\norm\left(\mathbf{u} \mid \mathbf{m}, \mathbf{S}\right)$, while keeping the exact conditional $p(\fvect \mid \mathbf{u})$ intact, that is $q(\fvect, \mathbf{u}) = p(\fvect \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(\yvect \mid \Xmat) = \int p(\yvect \mid \fvect) p(\fvect \mid \Xmat) d\fvect$.
The lower bound on $\log p\left(\yvect \mid \Xmat \right)$ can be obtained by considering the form of $q(\fvect, \mathbf{u})$ above and by applying Jensen's inequality:
\begin{equation}
\label{eq: lower bound on marginal likelihood 1}
\begin{split}
\E_{q\left(\fvect\right)}\left[\log p\left(\yvect \mid \fvect\right)\right] - \mathrm{KL}\left(q\left(\mathbf{u}\right) \| \, p\left(\mathbf{u}\right)\right) \text{.}
\end{split}
\end{equation}
The approximate posterior $q\left(\fvect\right)$ can be computed by integrating out $\mathbf{u}$:
$
q\left(\fvect\right) = \int q\left(\mathbf{u}\right) p\left(\fvect\left|\mathbf{u}\right.\right) d\mathbf{u}.
$
Thanks to the Gaussian form of $q\left(\mathbf{u}\right)$, $q\left(\fvect\right)$ can be computed analytically:
\begin{equation}
\label{eq: q(f)}
q\left(\fvect\right) = \norm\left(\fvect\ |\  \mathbf{A}\mathbf{m}, \ \Kmat_\Xmat + \mathbf{A}\left(\mathbf{S} - \Kmat_\mathbf{Z}\right)\mathbf{A}\right) \text{,}
\end{equation}
where $\mathbf{A} = \Kmat_{\Xmat, \mathbf{Z}} \Kmat_{\mathbf{Z}}^{-1}$.
When the likelihood factorizes over training points,

the lower bound can be re-written as:
\begin{equation}
\textstyle\sum_{i=1}^N\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{.}
\end{equation}
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}\Lmat^T$, with $\mathbf{L}$ lower triangular.

\subsection{Comparison to Local GPs}
We finally demonstrate that \name{swsgp} behaves differently from other approaches that use local approximations of \textsc{gp}s\xspace.
We consider two well-established approaches of local \textsc{gp}s\xspace proposed by \citet{Kim05} and \citet{Urtasun08}.
Following \citet{Liu2020WhenGP}, we refer to these methods as \emph{Inductive} \textsc{gp}s\xspace and \emph{Transductive} \textsc{gp}s\xspace, respectively.
We run all methods on two regression data sets: \name{powerplant} and \name{kin}.
We set the number of local experts to $64$, and we use the same number of inducing points for \name{swsgp} (with $H$ either 4 or 8).
As the size of \name{powerplant} and \name{kin} are approximately $7,000$, we set the number of training points governed by a local expert to $100$.
For the local \name{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 \name{gp} expert.
For the testing phase, inductive \textsc{gp}s\xspace simply rely on the nearest local experts to an unseen point $\xvect_*$.
For transductive \textsc{gp}s\xspace, we use $100$ neighbors of $\xvect_*$ and the nearest local expert to make predictions.
In Table \ref{tab: Comparison to Local approximation of gps}, we summarize \name{rmse} and \name{mnll} for all methods; \name{swsgp} clearly outperforms the local \name{gp} approaches in terms of \name{mnll}.

\begin{table}[!htb]
	\centering
    \caption{Comparison with Local \name{gp} approximations.}
    \vspace{8pt}
    \label{tab: Comparison to Local approximation of gps}
	\begin{tabular}{C{2.8cm}C{2.0cm}C{2.0cm}}
	\toprule
	Method & \name{powerplant} & \name{kin}\\
	       & \name{rmse} $\mid$ \name{mnll} & \name{rmse} $\mid$ \name{mnll}\\
	\midrule
	\name{swsgp}-64-4 & $4.27 \mid 2.41$ & $0.11 \mid -1.27$\\
	\midrule
	\name{swsgp}-64-8 & $\textbf{4.24} \mid \textbf{2.40}$ & $0.10 \mid \textbf{-1.38}$\\
	\midrule
	Inductive GPs & $9.93 \mid 38.38$ & $0.13 \mid -0.40$\\
	\midrule
	Transductive GPs & $6.17 \mid 18.78$ & $\textbf{0.09} \mid -0.65$\\
	\bottomrule
	\end{tabular}
\end{table}

\subsection{Comparing with Structured Kernel Interpolation (\name{ski})}
We conclude by reporting comparisons with other \name{gp}-based models.
In particular, we compare against the Stochastic Variational Deep Kernel Learning (\name{svdkl}) \citep{Wilson16} and the Deep \name{gp} approximated with random features (\textsc{dgp}-\textsc{rbf}\xspace) \citep{Cutajar17}.
In the former, \name{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 \name{gp} are approximated as parametric models using random feature expansions.
Both competitors feature mini-batch based learning, so this represents a challenging test for \name{swsgp}.
The results in Table \ref{tab: comparison of swsgp and ski} show that \name{swsgp} is comparable with these competitors on various data sets.
We believe that this is a remarkable result obtained by our shallow \name{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., \citet{Rudi17}).
\begin{table}[!htb]
\caption{Comparison of \name{swsgp}, \name{kiss-gp} \citep{Wilson15}, \name{svdkl}\citep{Wilson16} and \textsc{dgp}-\textsc{rbf}\xspace \citep{Cutajar17}.
	The results \name{swsgp} and \name{kiss-gp} are averaged over 5 folds.
In order to deal with the difficulties of \name{kiss-gp} to handle large-dimensional input spaces, we followed \citep{Wilson15a}, 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.
}
\label{tab: comparison of swsgp and ski}
\begin{subtable}{.98\linewidth}
	\centering
	\caption{\name{airline}}
	\label{tab: comparing ski airline}
	\begin{tabular}{C{2.9cm}C{2.0cm}C{2.0cm}}
		\toprule
		Method & \name{err} & \name{mnll}\\
		\midrule
		\name{swsgp}-{\footnotesize{100k}}\name{m}-{\footnotesize{100}}\name{h} & $\mathbf{0.210 \pm 0.012}$ & $0.48 \pm 0.015$\\
		\name{svdkl} & $0.22$ & $\mathbf{0.46}$ \\
		\textsc{dgp}-\textsc{rbf}\xspace & $\mathbf{0.21}$ & $\mathbf{0.46}$ \\
		\bottomrule
	\end{tabular}
\end{subtable}
\begin{subtable}{0.98\linewidth}
	\vspace{10pt}
	\centering
	\caption{\name{powerplant}}
	\label{tab: comparing ski powerplant}
	\begin{tabular}{C{2.3cm}C{2.1cm}C{2.1cm}}
		\toprule
		Method & \name{rmse} & \name{mnll}\\
		\midrule
		\name{swsgp}-{\footnotesize{64}}\name{m}-{\footnotesize{4}}\name{h} & $\mathbf{4.095 \pm 0.145}$ & $\mathbf{2.371 \pm 0.038}$ \\
		\name{kiss-gp} & $4.459 \pm 0.355$ & $3.074 \pm 0.037$\\
		\bottomrule
	\end{tabular}
\end{subtable}
\begin{subtable}{0.98\linewidth}
	\vspace{10pt}
	\centering
	\caption{\name{protein}}
	\label{tab: comparing ski protein}
	\begin{tabular}{C{2.3cm}C{2.1cm}C{2.1cm}}
		\toprule
		Method & \name{rmse} & \name{mnll}\\
		\midrule
		\name{swsgp}-{\footnotesize{64}}\name{m}-{\footnotesize{4}}\name{h} & $\mathbf{0.773 \pm 0.002}$ & $\mathbf{0.702 \pm 0.003}$\\
		\name{kiss-gp} & $0.816 \pm 0.001$ & $0.904 \pm 0.017$\\
		\bottomrule
	\end{tabular}
\end{subtable}
\end{table}

\section{Conclusions}
Sparse approaches that rely on inducing points
have met with success in reducing the complexity of \name{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 \name{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 \name{gp}-based approaches.

In this work, we further improve the computational gains of sparse \textsc{gp}s\xspace.
We proposed \name{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 \name{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 \name{swsgp} compared to the state-of-the-art.
\section{Experiments}
\label{sec:experiments}
We conduct experiments to evaluate \name{swsgp} on a variety of experimental conditions.
Our approach is denoted by \name{swsgp}-\name{m}-\name{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 \textsc{svgp}\xspace-\name{m} and \textsc{svgp}\xspace-\name{m}-\name{h} as competitors; \textsc{svgp}\xspace-\name{m} uses $M$ inducing points.
\textsc{svgp}\xspace-\name{m}-\name{h}, instead, refers to \textsc{svgp}\xspace using $M$ inducing points at training time and $H$-nearest inducing inputs at test time.
We also use \textsc{svgp}\xspace-{\footnotesize{16}\name{m}} using $16$ inducing points as a reference.
The comparison is carried out on some UCI data sets for regression and classification, i.e., \name{powerplant}, \name{kin}, \name{protein}, \name{eeg}, \name{credit} and \name{spam}.
We also consider larger scale data sets, such as, \name{wave}, \name{query} and the \name{airline} data or images classification on \mnist.
The task for the \name{airline} data set is the classification of whether flights are to subject to a delay, and we follow the same setup as in \citet{Hensman13} and \citet{Wilson16}.
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\'{e}rn-$\nicefrac{5}{2}$ kernel in all cases except for the \name{airline} dataset, where the sum of a Mat\'{e}rn-$\nicefrac{3}{2}$ and a linear kernel is used, similar to \citet{Hensman15b}.
All models are trained using the Adam optimizer \citep{Kingma14b} 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 (\name{rmse}) and the test mean negative log-likelihood (\name{mnll}), whereas we report the test error rate (\name{err}) and \name{mnll} in classification tasks.
The results are averaged over $10$ folds.
All models in section \ref{sec: listM-listH} are trained over $300,000$ iterations.
\input{increasing_M_H}
\input{runtime}
\input{large_M}
\input{comparing_ski}
\input{comparing_localgp}
\input{joint_predictive_covariances}

\subsection{Impact of M and H}
\label{sec: listM-listH}
We begin our evaluation by investigating the behavior of \name{swsgp}-\name{m}-\name{h} with respect to $M$ and $H$.
Figure \ref{fig: listM-listH} shows that \name{swsgp}-\name{m}-\name{h} consistently outperforms \textsc{svgp}\xspace-\name{m}-\name{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 \name{swsgp} are comparable with those obtained by \textsc{svgp}\xspace-\name{m}, while at each iteration only a small subset of $H$ out of $M$ inducing points are updated, carrying a significant complexity reduction.
\begin{figure}[!htb]
	\centering
	\begin{subfigure}{.9\linewidth}
		\centering

		\includegraphics[width=.98\linewidth]{04_powerplant_kin8nm_protein_listM_listH.pdf}

		\label{fig: listM-listH-regression}
	\end{subfigure}
	\begin{subfigure}{.9\linewidth}
		\centering

		\includegraphics[width=.98\linewidth]{05_eeg_credit_spam_listM_listH.pdf}

		\label{fig: listM-listH-binary-classification}
	\end{subfigure}
	\caption{
		Evaluation of \name{swsgp} on \name{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 \textsc{svgp}\xspace with $M$ inducing points.
		The cyan down-triangles are for \textsc{svgp}\xspace with $16$ inducing points.
		The red circles are for \textsc{svgp}\xspace training with $M$ inducing points and the prediction at an unseen data $\xvect$ are made by $\mathbf{Z}_\xvect^H$.
		The green squares, stars and plus are for \name{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.

	}
	\label{fig: listM-listH}
\end{figure}

\section{Introduction}
Gaussian Processes (\textsc{gp}s\xspace) \citep{Rasmussen06} offer a powerful framework to perform inference over functions; being Bayesian, \textsc{gp}s\xspace provide rigorous uncertainty quantification and prevent overfitting.
However, the applicability of \textsc{gp}s\xspace on large datasets is hindered by their computational complexity of $\bigO\left(N^3\right)$, where $N$ is the training size.
This issue has fuelled a considerable amount of research towards scalable \name{gp} methodologies that operate on a set of \emph{inducing variables} \citep{Candela05}.
In the literature, there is a plethora of approaches that offer different treatments of the inducing variables \citep{Lawrence02, Seeger03fastforward, Snelson05, Guzman07, Titsias09, Hensman13, Wilson15, Hensman15b}.
Some of the more recent approaches, such as Scalable Variational Gaussian Processes (\textsc{svgp}s\xspace) \citep{Hensman15b}, allow for the application of \textsc{gp}s\xspace to problems with millions of data-points.
In most applications of scalable \textsc{gp}s\xspace, these are approximated using $M$ inducing points, which results in a complexity of $\bigO\left(M^3\right)$.
It has been shown recently by \citet{Burt19a} that it is possible to obtain an arbitrarily good approximation for a certain class of \name{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 \textsc{gp}s\xspace 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 \textsc{svgp}\xspace 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 \emph{Sparse within a Sparse \name{gp}} (\name{swsgp}).
\name{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 \name{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 \citep{Molchanov17, Louizos17}.
To the best of our knowledge, this work is the first to explore these ideas as means to sparsify the inducing set in sparse \textsc{gp}s\xspace .

\subsection{Joint Predictive Covariances}
\begin{figure}[!htb]
	\centering

	\includegraphics[width=.98\linewidth]{12_swsgpbatch_powerplant.pdf}

	\caption
	{
		\name{swsgp-u} on \name{powerplant}.
		The horizontal axis shows the mini-batch size of the testing set.
		We evaluate \name{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 \name{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.
	}
	\label{fig: swsgpu on powerplant}
\end{figure}
In this section, we propose an extension of \name{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 $\name{swsgp-u}$, where \name{u} stands for ``union''.
In the experiments, we execute \name{swsgp-u}-{\footnotesize{512}}\name{m}-{\footnotesize{128}}m 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 \name{swsgp-u}-{\footnotesize{512}}\name{m}-{\footnotesize{128}}m against \textsc{svgp}\xspace-{\footnotesize{512}}\name{m} and full \textsc{gp}s\xspace.
Due to the computational complexity of full \textsc{gp}s\xspace, 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 \name{rmse} and \name{mnll}.
The results in Figure~\ref{fig: swsgpu on powerplant} are averaged by repeating the process $10$ times.
The results indicate that \name{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  \name{swsgp-u} and \textsc{svgp}\xspace by reporting the mean square error with respect to the covariance of full \textsc{gp}s\xspace.
These results show that \name{swsgp-u} provides a good approximation of joint predictive covariances.

\subsection{Large-scale Sparse GP Modeling with a Huge Number of Inducing Points}
Here we show that \name{swsgp} allows one to use sparse \textsc{gp}s\xspace with a massive number of inducing points without incurring a prohibitive computational cost.
We employ several large-scale data sets, i.e. \name{protein}, \name{eeg}, \name{wave}, \name{query} and \name{airline} with 5 millions training samples.
We test \name{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 \name{swsgp}.
We have attempted to run \textsc{svgp}\xspace with such large values of $M$ without success (out of memory in a system with 32GB of RAM).
Therefore, as a baseline we execute \textsc{svgp}\xspace-{\footnotesize{128}}\name{m} and \textsc{svgp}\xspace-{\footnotesize{512}}\name{m} and report the results of \textsc{svgp}\xspace with the configuration in \citet{Hensman15b}.

In \name{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, \name{swsgp}'s training phase requires $\bigO\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 \textsc{svgp}\xspace, we can successfully run \name{swsgp} with an unprecedented number of inducing point, e.g. $M = 100,000$.
\iffalse
\begin{table}[!htb]
	\caption{\name{swsgp} with large number of inducing points.}
	\label{tab: swsgp large M}
	\begin{subtable}{.98\linewidth}
		\centering
		\label{tab: protein}
		\begin{tabular}{C{2.8cm}C{2cm}C{2.2cm}}
			\toprule
			Method & \name{rmse} & \name{mnll}\\
			\midrule
			\textsc{svgp}\xspace-128M & $0.739 \pm 0.008$ & $0.656 \pm 0.012$\\
			\textsc{svgp}\xspace-512M & $0.712 \pm 0.010$ & $0.621 \pm 0.013$\\
			\name{swsgp}-3kM-50H & $0.641 \pm 0.002$ & $0.507 \pm 0.003$\\
			\bottomrule
		\end{tabular}
		\caption{\name{protein}-41157-4573-9}
	\end{subtable}
    \begin{subtable}{.98\linewidth}
    	\centering
    	\label{tab: eeg}
    	\begin{tabular}{C{2.8cm}C{2cm}C{2.2cm}}
    		\toprule
    		Method & \name{err} & \name{mnll}\\
    		\midrule
    		\textsc{svgp}\xspace-128M & $0.167 \pm 0.011$ & $0.360 \pm 0.014$\\
    		\textsc{svgp}\xspace-512M & $0.121 \pm 0.002$ & $0.274 \pm 0.004$\\
    		\name{swsgp}-5kM-50H & $0.034 \pm 0.001$ & $0.094 \pm 0.004$\\
    		\bottomrule
    	\end{tabular}
    	\caption{\name{eeg}-70000-4900-14}
    \end{subtable}
	\begin{subtable}{0.98\linewidth}
		\centering
		\label{tab: wave}
		\begin{tabular}{C{3cm}C{2cm}C{2cm}}
			\toprule
			Method & \name{rmse} & \name{mnll}\\
			\midrule
			\textsc{svgp}\xspace-128M & $0.756 \pm 0.004$ & $0.677 \pm 0.005$\\
			\textsc{svgp}\xspace-512M & $0.744 \pm 0.003$ & $0.661 \pm 0.004$\\
			\name{swsgp}-10kM-50H & $0.709 \pm 0.006$ & $0.613 \pm 0.009$\\
			\bottomrule
		\end{tabular}
		\caption{\name{wave}-250000-37999-32}
	\end{subtable}
	\begin{subtable}{.98\linewidth}
		\centering
		\label{tab: query}
		\begin{tabular}{C{2.8cm}C{2cm}C{2.2cm}}
			\toprule
			Method & \name{rmse} & \name{mnll}\\
			\midrule
			\textsc{svgp}\xspace-128M & $0.064 \pm 0.002$ & $-1.771 \pm 0.019$\\
			\textsc{svgp}\xspace-512M & $0.064 \pm 0.001$ & $-1.773 \pm 0.015$\\
			\name{swsgp}-10kM-50H & $0.046 \pm 0.003$ & $-2.124 \pm 0.066$\\
			\bottomrule
		\end{tabular}
		\caption{\name{query}-190000-9843-4}
	\end{subtable}
    \begin{subtable}{.98\linewidth}
    	\centering
    	\label{tab: query}
    	\begin{tabular}{C{5.0cm}C{1.2cm}C{1.2cm}}
    		\toprule
    		Method & \name{err}& \name{mnll}\\
    		\midrule
    		\textsc{svgp}\xspace-200M\cite{Hensman15b} & $0.34$ & $0.61$\\
    		\name{swsgp}-100kM-100H & $0.21$ & $0.48$\\
    		\bottomrule
    	\end{tabular}
    	\caption{\name{airline}-5000000-1000000-8}
    \end{subtable}
\end{table}
\fi

\begin{figure}[!htb]
	\centering
	\begin{subfigure}{.9\linewidth}
		\centering
		\includegraphics[width=.98\linewidth]{08_protein_runtime.pdf}

		\label{fig: protein runtime}
	\end{subfigure}
    \begin{subfigure}{.9\linewidth}
    	\centering
    	\includegraphics[width=.98\linewidth]{09_eeg_runtime.pdf}

    	\label{fig: eeg runtime}
    \end{subfigure}
	\begin{subfigure}{.9\linewidth}
		\centering

		\includegraphics[width=.98\linewidth]{10_query_runtime.pdf}

		\label{fig: query runtime}
	\end{subfigure}
    \begin{subfigure}{.9\linewidth}
    	\centering

    	\includegraphics[width=.98\linewidth]{11_wave_runtime.pdf}

    	\label{fig: wave runtime}
    \end{subfigure}
    \caption{\name{swsgp} with large number of inducing points. The figure shows the progression of \name{rmse} (\name{err}) and \name{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].}
	\label{fig: large-M}
\end{figure}
In \name{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 \name{swsgp}-{\footnotesize{100,000}}\name{m}-{\footnotesize{100}}\name{h} for one million iterations.
The \name{err} and \name{mnll} of \name{swsgp}-{\footnotesize{100,000}}\name{m}-{\footnotesize{100}}\name{h} evaluated on the test set are $21\%$ and $0.48$, respectively, while the \name{err} and \name{mnll} of \textsc{svgp}\xspace-{\footnotesize{200}}\name{m} published in \citet{Hensman15b} are about $34\%$ and $0.61$, respectively.
To the best of our knowledge, \name{swsgp} is the first to enable sparse \textsc{gp}s\xspace with such a large set of inducing points without imposing a grid structure on the inducing inputs.
In addition, in Figure~\ref{fig: large-M} we also show that \name{swsgp} using $3,000$ or $5,000$ or $10,000$ inducing points outperforms \textsc{svgp}\xspace using $128$ and $512$ inducing points on other data sets.

\section{Sparse Within Sparse GP
}
We present a novel formulation of sparse \textsc{gp}s\xspace, 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:
\begin{equation}
\label{eq: hierarchical prior}
p\left(\mathbf{u} \mid \mathbf{w} \right) = \norm\left(\zerovect, \mathbf{D}_\mathbf{w}\Kmat_{\mathbf{Z}}\mathbf{D}_\mathbf{w} \right),
\end{equation}
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.

\paragraph{Remarks on the prior over $\fvect$}
Our strategy simply assumes a certain structure on the auxiliary variables, but it \textit{has no effect} on the prior over $\fvect$; the latter remains unchanged.
Let $\mathcal{I}$ and $\mathcal{J}$ bet the sets of indices such that $\mathbf{w}_{\mathcal{I}} = \onevect$ and $\mathbf{w}_{\mathcal{J}} = \zerovect$.
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:
\begin{equation}
p\left(\fvect, \mathbf{u}\mid\mathbf{w}\right) =
\norm\left(
\begin{bmatrix}
\zerovect\\
\zerovect\\
\zerovect
\end{bmatrix}
,
\begin{bmatrix}
\Kmat_{\Xmat} & \Kmat_{\Xmat, \mathbf{Z}_{\mathcal{I}}} & \zerovect\\
\Kmat_{\mathbf{Z}_{\mathcal{I}}, \Xmat} & \Kmat_{\mathbf{Z}_{\mathcal{I}}} & \zerovect\\
\zerovect & \zerovect & \zerovect
\end{bmatrix}
\right)
\end{equation}
The rows and columns of $\mathbf{u}_{\mathcal{J}}$ can simply be ignored.
Regardless of the value of $\mathbf{w}$, the conditional $\fvect,\mathbf{u}_{\mathcal{I}}\mid\mathbf{w}$ is always a Gaussian marginal, as it is a subset of Gaussian variables.
The marginalized $p(\fvect,\mathbf{u}) = \int p(\fvect,\mathbf{u}\mid\mathbf{w})\, p(\mathbf{w}) d\mathbf{w}$ is mixture of Gaussian densities, where the marginal over $\fvect$ is the same for every component of the mixture.

\begin{figure}
\includegraphics[width=\linewidth]{01_prior_f_samples.pdf}
\caption{The choice of inducing points does not affect the prior samples drawn from $p(\fvect)$. Left: visualizations of $\fvect\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 $\fvect$, against the true $p(\fvect)$.}
\label{fig:prior_f}
\end{figure}

The effect on $\fvect$ is demonstrated in Figure \ref{fig:prior_f}, 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, $\fvect$ samples are drawn from $p(\fvect\mid\mathbf{u}_{\mathcal{I}})$, which only involves the selected inducing variables $\mathbf{u}_{\mathcal{I}}$.
Following Eq.~\eqref{eq:fu_joint}, the conditional $\fvect\mid\mathbf{u}$ is normally-distributed with mean
$\mathbf{m}_{\fvect\mid\mathbf{u}_{\mathcal{I}}} = \Kmat_{\Xmat, \mathbf{Z}_{\mathcal{I}}} \Kmat_{\mathbf{Z}_{\mathcal{I}}}^{-1} \mathbf{u}_{\mathcal{I}}$
and covariance
$\mathbf{S}_{\fvect\mid\mathbf{u}_{\mathcal{I}}} = \Kmat_\Xmat - \Kmat_{\Xmat, \mathbf{Z}_{\mathcal{I}}} \Kmat_{\mathbf{Z}_{\mathcal{I}}}^{-1} \Kmat_{\mathbf{Z}_{\mathcal{I}}, \Xmat}$.
These conditionals can be seen for different samples of $\mathbf{u}, \mathbf{w}$ in the left side of Figure \ref{fig:prior_f}, while in the right side we compare the marginalized prior over $\fvect$ against the true \name{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 \citep{Titsias09,Burt19a}.
Our choice to impose a hierarchical structure to the inducing variables through $\mathbf{w}$ effectively changes the model compared to \textsc{svgp}\xspace, and we adapt the variational scheme accordingly.

\subsection{Lower Bound on Marginal Likelihood}
By introducing $\mathbf{u}, \mathbf{w}$ and using Jensen's inequality, the lower bound on $\log p\left(\yvect\right)$ can be obtained as follows

\begin{equation}
\label{eq:Jensen1}
\begin{split}
\E_{q\left(\mathbf{u},\mathbf{w}\right)} \log p\left(\yvect\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}
\end{equation}
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:
\begin{equation}
q(\mathbf{u} \mid \mathbf{w}) = \norm\left(\mathbf{u}\left|\mathbf{D}_\mathbf{w} \mathbf{m}, \mathbf{D}_\mathbf{w} \mathbf{S} \mathbf{D}_\mathbf{w}\right.\right)
\label{eq: variational distribution q(u|w)}
\end{equation}
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(\yvect\left|\mathbf{u},\mathbf{w}\right.\right)$, obtaining:
\begin{equation}
\label{eq:Jensen2}
\log p\left(\yvect\left|\mathbf{u},\mathbf{w}\right.\right) \geq \E_{p\left(\fvect|\mathbf{u},\mathbf{w}\right)}\log p\left(\yvect\left|\fvect\right.\right)
\end{equation}
The bound we describe in Equations \eqref{eq:Jensen1} and \eqref{eq:Jensen2} is the same as in \textsc{svgp}\xspace, 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 \eqref{eq:Jensen2} into \eqref{eq:Jensen1}, 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}_\name{elbo}$:
\begin{equation}
\label{eq:ll_elbo}
	\begin{split}
		\sum_{n=1}^\nobs \E_{p\left(\mathbf{w}\right)} \bigg[&\E_{q\left(\mathbf{u}\left|\mathbf{w}\right.\right)} \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}
\end{equation}
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}}_\name{elbo}$:
\begin{equation}
	\begin{split}
		\sum_{n=1}^N \bigg[&\E_{q\left(\mathbf{u}\left|\tilde{\mathbf{w}}^{(n)}\right.\right)} \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}
\end{equation}
where $\tilde{\mathbf{w}}^{(n)}$ is sampled from $p\left(\mathbf{w}\right)$.
\paragraph{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 \textsc{svgp}\xspace with computational cost
of $\bigO\left(M^3\right)$ coming from computing $p\left(f_n\left|\mathbf{u},\tilde{\mathbf{w}}=\onevect\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 \name{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.

\subsection{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_\xvect$ as the set of $H$-nearest inducing inputs.
Intuitively, the prediction for an unseen data $\xvect$ using $\mathbf{Z}^H_\xvect$ is a good approximation of the prediction using all $M$ inducing points, that is $\mathbf{Z}^M_\xvect$.
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.~\eqref{eq: q(f)}.
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(\xvect\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 $\xvect$ that are selected in the random mini-batch; this materializes the sampling from the implicit distribution $p(\mathbf{w})$.
The maximization of the \name{elbo} in the setting described is summarized in Algorithm~\ref{alg: MAIP} (\name{swsgp}).
\begin{algorithm}
	\caption{Sparse within sparse \name{gp}  (\name{swsgp})}
	\label{alg: MAIP}
	\begin{algorithmic}[1]
		\REQUIRE $\mathcal{D}$, $H$, $M$.
		\ENSURE The optimum of trainable parameters $\thetavect$.
		\STATE \alglinelabel{alg: initialize parameters} Initialize $\thetavect$ 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}\Lmat^T$.
		\WHILE{stopping criteria is False}\label{alg: while}
			\STATE \alglinelabel{alg: initialize ell and kl} $\mathrm{ELL} \gets 0$ and $\mathrm{KL} \gets 0$.
			\STATE \alglinelabel{alg: extract minibatch-size} Sample mini-batch $\mathrm{B}$ of size $n_{\mathrm{B}}$ from $\mathcal{D}$.
			\FOR{$\left(\xvect_i, y_i\right) \in \mathrm{B}$} \alglinelabel{alg: for}
				\STATE \alglinelabel{alg: find neighbor inducing points} Find $\mathbf{Z}^H_{\xvect_i}$, i.e. the $H$-nearest $\mathbf{Z}$ to $\xvect_i$.
				\STATE \alglinelabel{alg: approximate posterior} Compute $w\left(\xvect_i\right)$ using $\mathbf{Z}^H_{\xvect_i}$ as in \eqref{eq: binary form of w(xn)}.
				\STATE \alglinelabel{alg: extract relevant parameters} Extract $\mathbf{m}_{w\left(\xvect_i\right)}$ and $\mathbf{S}_{w\left(\xvect_i\right)}$ from $\mathbf{m}$ and $\mathbf{L}$.
				\STATE \alglinelabel{alg: compute approximate predictive distribution} Compute $q\left(f_i\left|w\left(\xvect_i\right)\right.\right)$ as in \eqref{eq: compute posterior using H-nearest IP}.
				\STATE \alglinelabel{alg: update ell} $\mathrm{ELL} \gets \mathrm{ELL} + \E_{q\left(f_i\left|w\left(\xvect_i\right)\right.\right)}\log p\left(y_i\left|f_i\right.\right)$.
				\STATE \alglinelabel{alg: update kl} $\mathrm{KL} \gets \mathrm{KL} + \mathrm{KL}\left(q\left(\mathbf{u}_{w\left(\xvect_i\right)}\right) \| p\left(\mathbf{u}_{w\left(\xvect_i\right)}\right)\right)$
			\ENDFOR
			\STATE \alglinelabel{alg: update ll} $\tilde{\mathcal{L}}_\name{elbo} \gets \frac{N}{n_{\mathrm{B}}}\mathrm{ELL} - \frac{1}{n_{\mathrm{B}}}\mathrm{KL}$.
			\STATE \alglinelabel{alg: update parameters} Update $\thetavect$ using the derivative of $\tilde{\mathcal{L}}_\name{elbo}$.
		\ENDWHILE
	\end{algorithmic}
\end{algorithm}

\paragraph{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(\xvect\right)$.
In fact, if we would like to approximate the predictive distribution at $\xvect_n$ using $H$-nearest inducing inputs to $\xvect$, i.e. $\mathbf{Z}^H_{\xvect_n}$, then $w\left(\xvect\right) = \left[w^{(1)}_\xvect ... w^{(M)}_\xvect\right]^T$ where,
\begin{equation}
\label{eq: binary form of w(xn)}
w^{(m)}_\xvect =
\begin{cases}
1 & \mbox{if } \zvect_m \in \mathbf{Z}^H_{\xvect}\\
0 & \mbox{else}
\end{cases}
\text{ , with } m = 1, ..., M
\end{equation}
We extract the relevant elements using $w\left(\xvect\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(\xvect_i\right)} = \mathbf{D}_{w\left(\xvect_i\right)} \mathbf{S} \mathbf{D}_{w\left(\xvect_i\right)}$.
The approximate posterior over $f_i$ given $w\left(\xvect_i\right)$, i.e. $q\left(f_i\left|w\left(\xvect_i\right)\right.\right)$ is:
\begin{equation}
	\label{eq: compute posterior using H-nearest IP}
	\begin{split}
		\norm\bigg(f_i \mid &\mathbf{A}_{\xvect_i} \mathbf{m}_{w\left(\xvect_i\right)}, \\
		&\Kmat_{\xvect_i} + \mathbf{A}_{\xvect_i}\left(\mathbf{S}_{w\left(\xvect_i\right)} - \Kmat_{\mathbf{Z}^H_{\xvect_i}}\right)\mathbf{A}_{\xvect_i}^\top\bigg) \text{,}
	\end{split}
\end{equation}
where $\mathbf{A}_{\xvect_i} = \Kmat_{\xvect_i, \mathbf{Z}_{\xvect_i}^H}\Kmat_{\mathbf{Z}_{\xvect_i}^H}^{-1}$.

\paragraph{Limitations}
What we presented allows one to compute predictive distributions for individual test points following Eq.~\eqref{eq: compute posterior using H-nearest IP}.
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.~\eqref{eq: compute posterior using H-nearest IP} 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 \name{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 \name{swsgp} performs extremely well against various state-of-the-art competitors on a wide range of experiments.

\begin{comment}
{\color{red} comments on the effect on the model $\downarrow$ \\
In the definition of $\mathcal{L}_\name{elbo}$ in Eq.~\eqref{eq:ll_elbo} the quality of the model is reflected in the expectation likelihood.
We denote $q(\fvect\mid\mathbf{w})$ to be the approximate posterior of $\fvect$ fiven a subset of inducing variables indexed by $\mathbf{w}$.
Then we let $q(\fvect)$ be the approximate posterior for the entire set $\mathbf{u}$.
Ideally, we should have $\int q(\fvect\mid\mathbf{w}) p(\mathbf{w}) d\mathbf{w} \approx q(\fvect)$.
}
\end{comment}

\textbf{One-dimensional regression example.}
We visualize the posterior distribution for a synthetic dataset generated on a one-dimensional input space.
We execute \textsc{svgp}\xspace and \name{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 \ref{fig: posterior svgp, swsgp}.
We consider identical settings for the two methods (i.e.\ $128$ inducing points, kernel parameters, likelihood variance) and a neighbor area of $16$ for \name{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 \ref{sec:experiments}.

\begin{figure}[t]
	\begin{subfigure}{.98\linewidth}
		\centering

		\includegraphics[width=.98\linewidth]{02_svgp_128_synthetic.pdf}

		\label{fig: swsgp_irradiance_64_8 trainable Z}
	\end{subfigure}
	\\
	\begin{subfigure}{.98\linewidth}
		\centering

		\includegraphics[width=.98\linewidth]{03_swsgp_synthetic_128_16_fixedZ.pdf}

		\label{fig: swsgp_irradiance_64_8 fixed Z}
	\end{subfigure}
	\caption{Visualization of posterior distribution of \textsc{svgp}\xspace and \name{swsgp}. In both cases, we consider 128 inducing points; in terms of our scheme (\name{swsgp}) we use 16 neighbors.}
	\label{fig: posterior svgp, swsgp}
\end{figure}

\subsection{Complexity}
The computational cost of \name{swsgp} is dominated by lines \ref{alg: find neighbor inducing points}, \ref{alg: extract relevant parameters} and \ref{alg: compute approximate predictive distribution} in Algorithm \ref{alg: MAIP}.
For each data-point $\left(\xvect_i, y_i\right)$ in mini-batch $\mathrm{B}$, we need to find the $H$ nearest inducing neighbors $\mathbf{Z}^H_{\xvect_i}$ for $n_{\mathrm{B}}$ points in line \ref{alg: find neighbor inducing points}, where $n_{\mathrm{B}} = |\mathrm{B}|$; this contributes to the worst-case complexity by $\bigO\left(n_{\mathrm{B}} M H\right)$.

In line~\ref{alg: extract relevant parameters}, we extract relevant parameters from $\mathbf{m}$ and $\mathbf{L}$.
We focus on the cost of extracting $\mathbf{S}_{w\left(\xvect_i\right)}$ from $\mathbf{L}$.
Similar to \textsc{svgp}\xspace (Section \ref{subsection: svgp}), we consider $\mathbf{S}=\mathbf{L}\Lmat^T$, where $\mathbf{L}$ is lower triangular.
We extract $\mathbf{L}_{w\left(\xvect_i\right)} = \mathbf{D}_{w\left(\xvect_i\right)} \mathbf{L}$ which contains the rows of $\mathbf{L}$. Then, we compute $\mathbf{S}_{w\left(\xvect_i\right)}$ by $\mathbf{S}_{w\left(\xvect_i\right)} = \mathbf{L}_{w\left(\xvect_i\right)} \mathbf{L}_{w\left(\xvect_i\right)}^T$.
The computational complexity of selecting the variational parameters is $\bigO\left(n_{\mathrm{B}} M H^2\right)$.

Finally, the computation of approximating the predictive distribution in line \ref{alg: compute approximate predictive distribution} requires $\bigO\left(n_{\mathrm{B}} H^3\right)$.
The overall complexity for \name{swsgp} in the general case is $\bigO\left(n_{\mathrm{B}} M H + n_{\mathrm{B}} M H^2 + n_{\mathrm{B}} H^3\right)$, which is a significant improvement over the $\bigO\left(M^3\right)$ complexity of standard \textsc{svgp}\xspace, assuming that $n_{\mathrm{B}}, H \ll M$.
If we choose $\mathbf{S}$ to be diagonal, the total complexity reduces to $\bigO\left(n_{\mathrm{B}} M H + 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 \ref{sec:experiments} we explore all these settings.

\subsection{Running Time}
\iffalse
\begin{table}[!htb]
	\centering
	\caption{Comparison of running time between \textsc{svgp}\xspace and \name{swsgp}.
		In the table, each cell follows the format of [training time]$\mid$[testing time] (times are in milliseconds).
		In the figure, we show the progression of \name{err} (\name{rmse} for regression case) and \name{mnll} over training time.
		The black lines refer \textsc{svgp}\xspace, and the green lines indicate \name{swsgp}.}
	\begin{subtable}[tb]{0.98\linewidth}
		\centering
		\begin{tabular}{C{2.2cm}C{2.0cm}C{2.0cm}}
			\toprule
			Configuration & \name{powerplant} & \name{eeg}\\
			\midrule
			\textsc{svgp}\xspace-256 & $\textbf{22.83} \mid 2.89$ & $\textbf{21.42} \mid 1.43$\\
			\midrule
			\name{swsgp}-256-4 & $25.51 \mid \textbf{0.51}$ & $26.18 \mid \textbf{0.56}$\\
			\bottomrule
		\end{tabular}
	    \caption{Running times for \textsc{svgp}\xspace-256 and \name{swsgp}-256-4.}

	\end{subtable}

    \begin{subtable}[tb]{0.98\linewidth}
    	\centering
    	\begin{tabular}{C{2.2cm}C{2.0cm}C{2.0cm}}
    		\toprule
    		Configuration & \mnist & \name{airline}\\
    		\midrule
    		\textsc{svgp}\xspace-1024& $516 \mid 21.6$ & $465 \mid 45.8$\\
    		\name{swsgp}-1024-4 & $\textbf{233} \mid \textbf{1.77}$ & $\textbf{157} \mid \textbf{0.78}$\\
    		\bottomrule
    	\end{tabular}
        \caption{Running times for \textsc{svgp}\xspace-1024 and \name{swsgp}-1024-4.}

    \end{subtable}

    \begin{subtable}[tb]{0.98\linewidth}
    	\centering
    	\includegraphics[width=.98\linewidth]{figure/metrics_runtime_mnist.pdf}
    	\caption{\name{err} and \name{mnll} over training time (in hours).}
    \end{subtable}
\label{tab: running time}
\end{table}
\fi
\begin{table}[!htb]
	\centering
    \caption{Comparison of running time between \textsc{svgp}\xspace and \name{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. \name{err} and \name{mnll} are also shown in terms of training time.
    	These figures are averaged over 5 folds.
    	The models for \name{eeg} and \mnist are trained in 90 minutes and 6 hours respectively.
    	\vspace{5pt}
    }
	\label{tab: running time}
	\begin{subtable}[tb]{0.98\linewidth}
		\centering
		\caption{Running times for \textsc{svgp}\xspace-{\scriptsize{256}}\name{m} and \name{swsgp}-{\scriptsize{256}}\name{m}-{\scriptsize{4}}\name{h} on \name{eeg}.}
		\begin{tabular}{C{2.4cm}C{0.9cm}C{0.9cm}C{0.9cm}C{0.9cm}}
			\toprule
			Configuration & $t_1$(ms) & $t_2$(ms) & \name{err} & \name{mnll}\\
			\midrule
			\textsc{svgp}\xspace-{\footnotesize{256}}\name{m} & $\mathbf{21.42}$ & 1.43 & $\mathbf{0.17}$ & $\mathbf{0.34}$\\
			\name{swsgp}-{\footnotesize{256}}\name{m}-{\footnotesize{4}}\name{h} & 26.18 & $\mathbf{0.56}$ & 0.18 & 0.36\\
			\bottomrule
		\end{tabular}

	\end{subtable}

	\begin{subtable}[tb]{0.98\linewidth}
		\centering
		\vspace{5pt}
		\caption{Running times for \textsc{svgp}\xspace-{\scriptsize{1024}}\name{m} and \name{swsgp}-{\scriptsize{1024}}\name{m}-{\scriptsize{4}}\name{h} on \mnist.}
		\begin{tabular}{C{2.6cm}C{0.9cm}C{0.9cm}C{0.9cm}C{0.9cm}}
			\toprule
			Configuration & $t_1$(ms) & $t_2$(ms) & \name{err} & \name{mnll}\\
			\midrule
			\textsc{svgp}\xspace-{\footnotesize{1024}}\name{m} & 516 & 21.6 & 0.02 & 0.066\\
			\name{swsgp}-{\footnotesize{1024}}\name{m}-{\footnotesize{4}}\name{h} & $\mathbf{233}$ & $\mathbf{1.77}$ & $\mathbf{0.016}$ & $\mathbf{0.05}$\\
			\bottomrule
		\end{tabular}
	\end{subtable}

	\newline
	\begin{subtable}[tb]{0.9\linewidth}
		\centering
		\vspace{10pt}
		\includegraphics[width=.98\linewidth]{06_eeg_256M_runtime.pdf}

		\caption{\name{err} and \name{mnll} over training time (in minutes) on \name{eeg}.}
	\end{subtable}
    \begin{subtable}[tb]{0.9\linewidth}
    	\centering
    	\includegraphics[width=.98\linewidth]{07_mnist_runtime.pdf}

    	\caption{\name{err} and \name{mnll} over training time (in minutes) on \mnist.}
    \end{subtable}
\end{table}

In Table \ref{tab: running time}, we report the training and testing times of \name{swsgp} and \textsc{svgp}\xspace
In \textsc{svgp}\xspace, we set $M=256$ for \name{eeg}, and $1024$ for \mnist, i.e. \textsc{svgp}\xspace-{\footnotesize{256}}\name{m} and \textsc{svgp}\xspace-{\footnotesize{1024}}\name{m}.
In our approach, we use the same $M$ and we set $H$ to $4$, i.e. \name{swsgp}-{\footnotesize{256}}\name{m}-{\footnotesize{4}}\name{h} and \name{swsgp}-{\footnotesize{1024}}\name{m}-{\footnotesize{4}}\name{h}.
In Table \ref{tab: running time}, we stress that $t_1$ and $t_2$ in \name{swsgp} take into account the computation of finding neighbors inducing inputs for each data-point.
In \textsc{svgp}\xspace, we assume that $\Kmat_{\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 $\bigO\left(M^2\right)$.
The time $t_2$ in \textsc{svgp}\xspace refers to the execution time of carrying out predictions with the complexity of $\bigO\left(M^2\right)$.

The results in  Tab \ref{tab: running time} show a consistent improvement at test time compared to \textsc{svgp}\xspace 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, \name{swsgp} offers limited improvements when $M$ is small.
Considering \name{eeg} in which $M$ is set to $256$, \textsc{svgp}\xspace is faster than \name{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 \ref{tab: running time} shows dramatic speedups compared to \textsc{svgp}\xspace when the number of inducing points $M$ is large.
When $M=1024$ on \mnist, \name{swsgp}-{\footnotesize{1024}}\name{m}-{\footnotesize{4}}\name{h} is faster than \textsc{svgp}\xspace-{\footnotesize{1024}}\name{m} at training time.
This is due to the inversion of the $1024 \times 1024$ kernel matrix being a burden for \textsc{svgp}\xspace, whereas \name{swsgp} deals with much cheaper computations.
In addition, we show the corresponding \name{err} and \name{mnll} of each model when we train \textsc{svgp}\xspace and \name{swsgp} on \name{eeg} and \mnist.
On the \name{eeg} data set, our method is comparable with \textsc{svgp}\xspace.
On \mnist,
\name{swsgp} reaches high accuracies significantly faster compared to \textsc{svgp}\xspace.

\section{Related Work and Background}
\input{related_work}
\input{background}

\input{methods}
\input{experiments}
\input{conclusions}

\paragraph{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).

\newpage