Convergence of Sparse Variational Inference in Gaussian Processes Regression

Introduction

Gaussian process (GP) priors are commonly used in Bayesian modelling due to their mathematical convenience and empirical success. The resulting models give flexible mean predictions, as well as useful estimates of uncertainty. GP priors are often used with a Gaussian likelihood for regression tasks, as the Bayesian posterior can be computed in closed form in this case. Additionally, in many instances, the kernel is differentiable with respect to hyperparameters, in which case hyperparameters can be efficiently learned using gradient-based optimization by maximizing the marginal likelihood, which can be computed analytically (also known as empirical Bayes, or type-II maximum likelihood). However, standard implementations of exact inference in Gaussian process regression models require storing and inverting a kernel matrix, imposing an $\mathcal{O}(N^{2})$ memory cost and an $\mathcal{O}(N^{3})$ computational cost, where $N$ is the number of training examples. These computational constraints have pushed researchers to adopt approximate methods in order to allow Gaussian process models to scale to large data sets.

Sparse methods (e.g. Seeger et al., 2003; Snelson and Ghahramani, 2006; Titsias, 2009b) rely on a set of inducing variables to represent the posterior distribution. While these methods have been widely adopted in research and application areas, there is a limited theoretical understanding of the effects of these approximations on the quality of posterior predictions, as well as what biases are introduced into hyperparameter selection when using approximations to the marginal likelihood. In this work, we aim to characterize the accuracy of sparse approximations. If all of the key properties of the exact model, i.e. the predictive mean and uncertainties and the marginal likelihood, are maintained by very sparse models, then a great deal of computation can be saved through these approximations.

We focus on the case of sparse inference in the variational framework of Titsias (2009b). We analyze the relationship between the level of sparsity used in performing inference, which dictates the computational cost, and the quality of the approximate posterior distribution. In particular, we analyze how many inducing variables should be used in order for the KL-divergence between the approximate posterior and the Bayesian posterior to be small. This offers theoretical insight into the trade-off between computation and quality of inference within the variational framework. From a practical perspective, our work suggests new methods for choosing which inducing variables to use to construct the approximation and provides theoretically grounded insight into the types of problems to which the sparse variational approach is particularly well-suited.

We derive bounds on the quality of variational inference in Gaussian process models. When our bounds are applied in the case of the squared exponential (SE) kernel and Gaussian or compactly supported inputs, we prove that the variational approximation can be made arbitrarily close to the true posterior with arbitrarily high probability using $\mathcal{O}((\log N)^{D})$ inducing variables, where $D$ is the dimensionality of the training inputs, leading to an overall computational cost of $\mathcal{O}(N(\log N)^{2D}(\log\log N)^{2})$. Note that we consider $D$ fixed throughout, implying a scaling in $N$ that is nearly linear, i.e. $\mathcal{O}(N^{1+\epsilon})$, $\forall\epsilon>0$.

Our bounds measure the discrepancy to the true posterior using the KL-divergence between the approximate and exact posteriors. We also show that this implies convergence of the point-wise predictive means and variances.

We show that theoretical guarantees on the quality of matrix approximation for existing methods for selecting regressors in sparse kernel ridge regression can be directly translated into guarantees on variational sparse GP regression. We demonstrate this for ridge leverage scores.

We derive lower bounds on the number of inducing variables needed to ensure that the KL-divergence remains small. For the SE kernel and Gaussian covariate distribution, these lower bounds have the same dependence on the size of the data set as the upper bounds.

Based on the theoretical results, we provide recommendations on how to select inducing variables in practice, and demonstrate empirical improvements.

This paper is an extension of the work Burt et al. (2019) presented at ICML 2019.

2 Overview of this Paper

In Section 2 we introduce notation and review the Gaussian process regression model, as well as sparse variational inference for Gaussian process models. In Section 3, we discuss practical considerations regarding assessing the quality of sparse variational inference using upper bounds on the log marginal likelihood that can be computed after observing a data set. In Section 4, we prove our main results, which bound the quality of the sparse approximate posterior, as measured by the KL-divergence. In order to do this, we consider methods for selecting inducing inputs inspired by methods used to obtain theoretical guarantees on sparse kernel ridge regression. Section 5 considers specific, commonly studied kernels and covariate distributions and investigates the implications of our results in these instances. We provide concrete computational complexities for finding arbitrarily accurate approximations to GPs. In Section 6, we consider the inverse problem, and show that in certain instances the KL-divergence will be large unless the number of inducing variables increases sufficiently quickly as a function of the size of the data set. Section 7 discusses practical insights and limitations of the theory as applied to real-world problems.

Background and Notation

In this section, we review exact inference in Gaussian process models, as well as sparse methods for approximate inference in these models. We particularly focus on the formulation of sparse methods based on variational inference (Titsias, 2009b). Throughout the paper, we use boldface letters to denote random variables, and the same letter in non-bold to denote a realization of this random variable. We follow the standard shorthand notation adopted in many Bayesian machine learning papers and denote probability densities by lower case letters $p$ and $q$, with the distribution to which they are associated inferred by the name of the argument; e.g. $p(X,y)$ is the density of a joint distribution over random variables $\mathbf{X}$ and $\mathbf{y}$ evaluated at $\mathbf{X}=X$ and $\mathbf{y}=y$.

2 Gaussian Process Regression

We specify our Bayesian model through a prior and likelihood. We place a GP prior, which for notational convenience we assume has zero mean function i.e. $\mu\equiv 0$, on the function $\mathbf{f}$ so that

To allow for deviations from $\mathbf{f}$ in the observations, we model the data $y$ as a noisy observation of this process through the likelihood

where the noise variance $\sigma^{2}$, is a model hyperparameter and I is an $N\times N$ identity matrix.

Both $p(f_{X^{\star}}\,|\,f_{X})$ and $p(f_{X}\,|\,\mathcal{D})$ are Gaussian densities and the marginal distribution of a Gaussian is Gaussian, so $p(f_{X^{\star}}\,|\,\mathcal{D})$ is also a Gaussian density. The posterior predictive distribution over the inputs $\mathbf{f}_{X^{\star}}$ has mean vector and covariance matrix

where $\textup{K}_{\star\text{f}}$ is $T\times N$ matrix with $[\textup{K}_{\star\text{f}}]_{t,n}=k(x^{\star}_{t},x_{n})$ and $\textup{K}_{\star\star}$ is a $T\times T$ matrix with $[\textup{K}_{\star\star}]_{t,t^{\prime}}=k(x^{\star}_{t},x^{\star}_{t^{\prime}})$.

The marginal likelihood is of interest in Bayesian models for selecting the properties of the model, which are determined by hyperparameters. Point estimates of model hyperparameters are commonly obtained by maximizing the marginal likelihood with respect to the noise variance $\sigma^{2}$, and any parameters of the prior covariance function $k$. In the case of conjugate regression described above, the log marginal likelihood takes the form

The quadratic term measures how well the data $y$ lines up with degrees of variation that are allowed under the prior. The log-determinant term measures how much variation there is in the prior, and penalizes priors which are widely spread. The combination of these terms in the log marginal likelihood balances the ability of the model to fit the data with model complexity, which allows a suitable model to be chosen; see Rasmussen and Williams (2006) for more discussion of the marginal likelihood as a tool for model selection as well as an introduction to Gaussian processes.

Despite closed-form expressions for both the predictive posterior (Eq. 5) and marginal likelihood (Eq. 6), exact inference in Gaussian process regression models is impractical for large data sets due to the cost of storing and inverting the kernel matrix $\textup{K}_{\textup{ff}}$, leading to $\mathcal{O}(N^{2})$ memory and $\mathcal{O}(N^{3})$ time complexities. Sparse approximations have been widely adopted to address this issue.

3 Approximate Inference for Gaussian Processes

Approximate inference in Gaussian process regression is performed for a different reason than in most Bayesian models. Approximate inference is usually applied when the exact posterior is analytically intractable. In our case, we can analytically write down the posterior, but the cost of computation is often prohibitive. The methods we discuss here all approximate the posterior with a different Gaussian process which has more favorable computational properties. As this approximate posterior has a similar form to the exact posterior, and we can control the trade-off between accuracy and computation, it is plausible that our approximation may be very accurate.

The large cost of computing the posterior GP comes from needing to infer a Gaussian distribution for the function values at all $N$ input locations. Sparse approximations (Seeger et al., 2003; Snelson and Ghahramani, 2006; Titsias, 2009b) avoid this cost by instead computing an approximate posterior that only depends on the data through the process at $M\ll N$ locations.

The aim of these methods is to compress the combined effect of a large number of input and output pairs into a distribution over function values at a small set of inputs. In regions where data is dense, there is often redundant information about what the function is actually doing, so little is lost in performing this approximation. The selected input locations and their corresponding function values are named inducing inputs and outputs respectively, and together are named inducing points. Later, it was suggested that more general linear transformations of the process could also be used to compress knowledge into (Lázaro-Gredilla and Figueiras-Vidal, 2009). We generally refer to these approaches as inducing variable methods. In all of these methods, a low-rank matrix appears in place of $\textup{K}_{\textup{ff}}$ in the computation of the posterior predictive and log marginal likelihood. This matrix can be manipulated with a much lower computational cost than working with $\textup{K}_{\textup{ff}}$ directly.

The success of inducing variable methods depends heavily on which $M$ random variables are chosen to represent the knowledge about the function $f$. Because in this work we are concerned with characterizing how large $M$ should be, we need a good method for choosing the inducing variables, as well as a meaningful criterion for judging the quality of the resulting approximation. The variational formulation of Titsias (2009b) is of particular interest, as it uses a well-defined divergence for characterizing the quality of the posterior, which can also be used as a guide for selecting the inducing variables.

3.2 The Variational Formulation

Variational inference proceeds by defining a family of candidate distributions $\mathcal{Q}$, and then selecting the distribution $Q\in\mathcal{Q}$ that minimizes the KL-divergence between the approximation and the posterior. In practice, elements of $\mathcal{Q}$ are parameterized and the approximate posterior is selected by choosing an initial approximation which is then refined by finding a local minimum of the KL-divergence as a function of the variational parameters. In variational GP methods (Titsias, 2009b; Hensman et al., 2013) $\mathcal{Q}$ consists of GPs with finite dimensional marginal densities of the form

for any $X^{\prime}\subset\mathcal{X},|X^{\prime}|<\infty$, where $q$ is the density of the approximate posterior at this collection of points, and $p(f_{X^{\prime}}\,|\,U)$ is the density of the prior distribution of $\mathbf{f}_{X^{\prime}}$ at $f_{X^{\prime}}$ conditioned on the random variables $\mathbf{U}$ evaluated at $\mathbf{U}=U$. In inducing point approximations, we take the inducing variables to be point evaluations of $f$, i.e. $\mathbf{U}=\mathbf{f}_{Z}$, with inducing inputs $Z\subset\mathcal{X}$ and $|Z|=M$.

As discussed in the previous section, we can also define inducing variables as linear transformations of the prior process of the form

where we assume $\rho$ is a measure on $\mathcal{X}$ defined with respect to an appropriate $\sigma$-algebra and $g_{m}\in L^{1}(\mathcal{X},\rho)$. If $\rho$ is taken to be a discrete measure, then these features correspond to (weighted) sums of inducing points; while other forms of these inducing variables of this form have been explored (Lázaro-Gredilla and Figueiras-Vidal, 2009; Hensman et al., 2018).

The density $q(U)$ is chosen to be an $M$-dimensional Gaussian density. This choice of variational family induces a Gaussian process approximate posterior with mean and covariance functions

where $\mu_{\text{U}},\Sigma_{\text{U}}$ are the mean and covariance of $q(U)$, $\textup{K}_{\textup{uu}}$ is the $M\times M$ matrix with entries $[\textup{K}_{\textup{uu}}]_{m,m^{\prime}}=\text{cov}(\mathbf{u}_{m},\mathbf{u}_{m^{\prime}})$, $k_{f(x)\text{u}}$ is the row vector with entries $[k_{f(x^{\prime})\text{u}}]_{m}=\text{cov}(\mathbf{f}(x),\mathbf{u}_{m})$ and $k_{\text{u}f(x)}$ is a column vector defined similarly. The variational parameters consist of $Z$, which determines the random variables that are included in $\mathbf{U}$, and $\mu_{\text{U}}$ and $\Sigma_{\text{U}}$, which determine the distribution over $\mathbf{U}$.

where $P$ denotes the (exact) posterior process (Matthews et al., 2016).

When the likelihood is isotropic Gaussian, the unique optimum for the parameters $\{\mu_{\text{U}},\Sigma_{\text{U}}\}$ can be computed in closed form. Using these optimal values, we obtain the ELBO as it was introduced by Titsias (2009b),

3.3 Measuring the Quality of a Variational Approximation

Variational approximations using this KL-divergence have been criticized for failing to provide guarantees on important quantities such as posterior estimates of the mean and variance. Huggins et al. (2019) observed that there exist Gaussian distributions such that the (normalized) difference between the means of the distributions is exponentially large as a function of the KL-divergence between the two distributions, as is the ratio of the variances. This has been used to motivate variational approaches based on other notions of divergence, as well as a more careful assessment of the quality of the approximations obtained via variational inference (Huggins et al., 2020).

However, in our case of sparse Gaussian process regression, a sufficiently small KL-divergence between the approximate and true posterior implies bounds on the approximation quality of the marginal posterior mean and variance function. Proposition 1 states one such bound:

The proof (Section A) uses that the KL-divergence between any pair of joint distributions upper bounds the KL-divergence between marginals of these distributions. It then suffices to bound the difference between the mean and variance of univariate Gaussian distributions with a small KL-divergence between them.

Proposition 1 implies that in cases where we can prove the KL-divergence between the approximate posterior and the exact posterior is very small, we are guaranteed to obtain similar marginal predictions with the variational approximation to those we would obtain with the exact model. We note that direct approaches to bounding marginal moments may lead to tighter bounds on these quantities (e.g. Calandriello et al., 2019), but we prefer to consider the KL-divergence due to its connection to the variational objective function.

The consequences of a small KL-divergence for hyperparameter selection using the evidence lower bound are more subtle, as both the approximate posterior and exact posterior depend on model hyperparameters, and it is generally difficult to ensure that the KL-divergence is uniformly small. We will be discuss these issues in more detail in Section 7.

3.4 Computation and Accuracy Trade-Offs

The ELBO (Eq. 9) as well as the corresponding choices for $\mu_{\text{U}}$ and $\Sigma_{\text{U}}$ (needed for making predictions) can be computed in $\mathcal{O}(NM^{2})$ time, and with $\mathcal{O}(NM)$ space. If a good approximation can be found with $M\ll N$, the savings in computational cost are large compared to exact inference. From Eq. 9 we see that the approximation is perfect when choosing $Z=X$, as this leads to $\textup{Q}_{\textup{ff}}=\textup{K}_{\textup{ff}}$. However, no computation is saved in this setting. We seek a more complete understanding of the trade-off between accuracy and computational cost when $M<N$ by understanding how $M$ needs to grow with $N$ to ensure an approximation of a certain quality. We derive probabilistic upper and lower bounds on this rate that depend on kernel properties that can be analyzed before observing any data.

4 Spectrum of Kernels and Mercer’s Theorem

In the previous section, we noted that sparse methods imply a low-rank approximation $\textup{Q}_{\textup{ff}}$ to the kernel matrix $\textup{K}_{\textup{ff}}$. In order to understand the impact of sparsity on the variational posterior, it is necessary to understand how well $\textup{K}_{\textup{ff}}$ can be approximated by a rank-$M$ matrix. This depends on the behavior of the eigenvalues of $\textup{K}_{\textup{ff}}$.

For small data sets, an eigendecomposition of $\textup{K}_{\textup{ff}}$ allows direct empirical analysis. However, for problems where sparse approximations are actually of interest, eigendecompositions are not available within our computational constraints. However, even without access to a specific data set, we can reason that properties of the training inputs have a large impact on the properties of the eigendecomposition of the kernel matrix. For example, consider the case of a squared exponential kernel given by

If the kernel is continuous and bounded, then $\mathcal{K}$ has countably many eigenvalues. We denote these eigenvalues in non-increasing order, so that $\lambda_{1}\geq\lambda_{2}\geq\dots\geq 0$. Corresponding to each non-zero eigenvalue $\lambda_{m}$ there is an eigenfunction $\phi_{m}$ which can be chosen to be continuous.

where the sum on the left converges absolutely and uniformly.See Rasmussen and Williams (2006), section 4.3 for more discussion of Mercer’s theorem.

The bounds we derive in the remainder of this work will depend on how rapidly the eigenvalues $\{\lambda_{m}\}_{m=1}^{\infty}$ decay. As they are absolutely summable, they must decay faster than $1/m$. The decay of these eigenvalues is closely related to the complexity of the non-parametric model as well as the generalization properties of the posterior (Micchelli and Wahba, 1979; Plaskota, 1996). Generally, these eigenvalues decay faster for covariate distributions that are concentrated in a small volume, and for kernels that give smooth mean predictors (Widom, 1963, 1964). Therefore, the bounds we prove in Section 4 can be seen as verifying the intuition that sparse variational approximations can be successfully applied to models with smooth prior kernels, as well as data sets with densely clustered covariates.

5 Inducing Variable Selection and Related Bounds

While the kernel eigenvalues determine how well a kernel matrix can be approximated, the quality of an actual approximation depends on how the inducing variables are chosen. Inducing point selection has been widely studied for many methods that require constructing a Nyström approximation, like sparse Gaussian processes and kernel ridge regression (KRR). In the simplest case, a subset can be uniformly sampled from the training inputs. Bounds on the quality of the resulting matrix approximation, and downstream Kernel Ridge Regression predictor have been found for this case (Bach, 2013; Gittens and Mahoney, 2016) and depend heavily on assumptions about the covariate distribution and resulting kernel matrix. In the Gaussian process literature, some specific low-rank parametric approximations based on spectral information about the kernel operator or matrix have been proposed (Zhu et al., 1997; Ferrari-Trecate et al., 1999; Solin and Särkkä, 2020) together with analysis on the rate of decrease in error with additional features. However, these methods generally are either limited in the types of kernels they can be applied to or have higher computational complexity than inducing point methods.

Heuristic inducing point selection methods have also been proposed in the hope of improving performance, for instance approximately minimizing $\textup{tr}(\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}})$ (Smola and Schölkopf, 2000), approximating the information gain of including a data point in the posterior (Seeger et al., 2003), or using the k-means centres of the input distribution (Hensman et al., 2013, 2015).

Two methods from the KRR literature are of particular interest: sampling from a Determinantal Point Process (DPP) (Li et al., 2016), and ridge leverage scores (Alaoui and Mahoney, 2015; Rudi et al., 2015; Calandriello et al., 2017). Theoretical guarantees exist in the literature for these methods applied to KRR, as well as empirical evidence of their efficacy compared to uniform sampling. The initial version of this work (Burt et al., 2019) analyzed convergence of the sparse variational GP posterior and marginal likelihood using the DPP initialization. Concurrently, Calandriello et al. (2019) used ridge leverage scores to show the DTC approximation (Seeger et al., 2003; Quiñonero-Candela and Rasmussen, 2005) can be made similar to the true posterior, in terms of pointwise predictive means and variances. Given the similarity between the DTC and variational posteriors, we include an analysis of ridge leverage sampling in this extended work to also provide results of convergence of the ELBO, and of the posterior in terms of the KL, which also implies pointwise convergence of the predictive means and variances.

Assessing Variational Inference: a Posteriori Bounds on the KL-divergence

We begin our investigation by considering how to choose the number of inducing variables for a specific data set. The simplest approach to assessing whether sufficiently many inducing points are used is to gradually increase the number of inducing points, and assess how the evidence lower bound changes with each additional point. If the ELBO increases only slightly or not at all when an additional inducing point is added, it is tempting to conclude that the approximate posterior is very close to the exact posterior. However, this is not a sufficient condition for the approximation to have converged. It could be the case that the last inducing point placed was not placed effectively, or that increasing from $M$ to $M+1$ inducing points has little impact, but increasing to $M+c$, for some $c>1$, inducing points would lead to significantly better performance if these points are well-placed.

A more refined mechanism for assessing the quality of the variational posterior would be to consider an upper bound on the KL-divergence that can be computed in similar computational time to the ELBO. Such a bound was proposed by Titsias (2014) and discussed as a method for assessing convergence in Kim and Teh (2018). In order to state this bound, we first need to introduce some notation. Let $t\coloneqq\textup{tr}(\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}})$ denote the trace of $\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}}$ and $\|\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}}\|_{\textup{op}}$ denote the operator norm of $\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}}$, which in this case is equal to the largest eigenvalue of this matrix as it is symmetric positive semidefinite.

For completeness we give a brief derivation of Lemma 2 in Section B, which essentially follows the derivation of Titsias (2014).

In problems where sparse GP regression is applied, computing the largest eigenvalue of $\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}}$ in order to compute $\mathcal{U}_{1}$ is computationally prohibitive. However, $\textup{tr}(\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}})$ can be computed in $\mathcal{O}(NM^{2})$, so that $\mathcal{U}_{2}$ can be computed efficiently.

If the difference between the upper and lower bounds is small, we can therefore be sure that sufficiently many inducing points are being used for the KL-divergence to be small. This suggests a refinement of the method for selecting the number of inducing points discussed earlier: continue to place more inducing points until the difference between the upper and lower bounds is small.

This raises the question: how many inducing variables do we need for the KL-divergence to be small in a typical problem? The upper bounds discussed above assess the approximation a posteriori, i.e. for a given data set and a given approximation. We would like to characterize the required number of inducing variables for a whole class of problems, before observing any data. This allows us to understand a priori how much computation is needed to solve a particular problem. For example, if the number of inducing variables $M$ needs to grow linearly with the number of observations $N$, then the $\mathcal{O}(NM^{2})$ cost of the approximation effectively scales cubically in $N$, i.e. in the same way as the exact implementation. In Section 4, we show that under intuitive assumptions, the number of inducing points can be taken to be much smaller than the size of the data set, while still giving approximations with small KL-divergences.

Convergence of Sparse Variational Inference in Gaussian Processes

In this section, we prove upper bounds on the KL-divergence between the approximate posterior and the exact posterior that depend on the number of inducing points used in inference, properties of the prior and distributional assumptions on the training covariates. The proof proceeds in three parts:

Derive an upper bound on the KL-divergence for a fixed data set and fixed set of inducing points that only depends on the quality of the approximation of $\textup{K}_{\textup{ff}}$ by $\textup{Q}_{\textup{ff}}$. In order to do this we make assumptions about the data generating process for $y$.

Suggest a method for selecting inducing inputs that obtains a high quality low-rank approximation to $\textup{K}_{\textup{ff}}$. This yields an upper bound on the KL-divergence depending only on the eigenvalues of $\textup{K}_{\textup{ff}}$. We consider using a $k$-determinantal point process or ridge leverage scores as the initialization method.

Relate eigenvalues of the kernel matrix back to those of the corresponding kernel operator, Eq. 10, through assumptions on the distribution of the covariates.

The second step has precedent in the literature on sparse kernel ridge regression. For example, Li et al. (2016) consider using a $k$-DPP to select the sparse regressors. Meanwhile ridge leverage scores have been studied in the setting of sparse kernel ridge regression and Gaussian process regression (Alaoui and Mahoney, 2015; Rudi et al., 2015; Calandriello et al., 2017, 2019), and have been shown to lead to strong statistical guarantees.

The third step in our analysis is similar to the analysis carried out when studying generalization and approximation bounds for Gaussian processes and other kernel methods. We use a generalization of a lemma proven in Shawe-Taylor et al. (2005) for this step.

In order to carry out our analysis, especially steps 2 and 3, we will treat $X,y$ and $Z$ as realizations of random variables $\mathbf{X}$, $\mathbf{y}$ and $\mathbf{Z}$ and make distributional assumptions about these random variables. This will allow us to make statements about bounds that hold in expectation or with fixed probability.

In Section 3, we considered bounds on the KL-divergence that can be computed for a specific data set. In this section, we first derive an upper bound on the KL-divergence that only depends on the squared norm of $\mathbf{y}$, with no additional assumptions on the distribution of the $\mathbf{y}$ (Lemma 3). We then derive a second bound, given in Lemma 4, that improves on Lemma 3 in expectation, under the stronger assumption that $\mathbf{y}|\mathbf{Z},\mathbf{X}\sim\mathcal{N}(0,\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I})$. This assumption is satisfied if $\mathbf{y}$ is distributed according to the prior model and the distributions of $\mathbf{Z}$ and $\mathbf{y}$ are independent, i.e. the inducing inputs are chosen without reference to $y$. While our results are stated in terms of inducing points, the proofs generalize without modification to other inducing variables of the form discussed in Section 2.3.2.

We first consider the case where we make few assumptions on the distribution of $\mathbf{y}$.

with $t=\text{\emph{tr}}(\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}})$ and $\zeta=\|\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}}\|_{\textup{op}}$.

The first inequality has already been established (Eq. 13). The remainder of the proof, given in Section C relies on properties of symmetric positive semi-definite (SPSD) matrices.

1.2 Average Case Analysis for the Prior Model

In Lemma 3, we did not make any assumption on the distribution of $\mathbf{y}$. From the Bayesian perspective, it is natural to make stronger distributional assumptions on $\mathbf{y}|\mathbf{X}$. We will see that in some instances stronger assumptions can lead to a much tighter upper bound than Lemma 3 that holds in expectation.

The natural candidate distribution for $\mathbf{y}$ is the prior distribution, that is $\mathbf{y}|\mathbf{X}\sim\mathcal{N}(0,\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I})$; if we additionally assume that the distributions of $\mathbf{y}|\mathbf{X}$ and $\mathbf{Z}|\mathbf{X}$ are independent, then this implies $\mathbf{y}|\mathbf{X},\mathbf{Z}\sim\mathcal{N}(0,\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I})$. In this case we can derive upper and lower bounds on the conditional expectation of the KL-divergence conditioned on $\mathbf{X}$ and $\mathbf{Z}$.

Suppose $\mathbf{y}|\mathbf{X},\mathbf{Z}\sim\mathcal{N}(0,\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I})$. For any $X\in\mathcal{X}^{N}$ and $Z\in\mathcal{X}^{M}$,

where $t=\text{\emph{tr}}(\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}})$ and $\textup{K}_{\textup{ff}}$ and $\textup{Q}_{\textup{ff}}$ are defined with respect to this $X,Z$ as in Section 2.

Note that if $\mathbf{y}\sim\mathcal{N}(0,\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I})$,

Therefore, under the strong assumption that $\mathbf{y}$ is sampled from the prior model, Lemma 4 gives a significantly stronger bound on the expected KL-divergence as compared to Lemma 3.

Let $n(y;m,S)$ denote the density of a (multivariate) Gaussian random variable with mean $m$ and covariance matrix $S$ evaluated at $y$. Then,

2 Initialization of Inducing Points

We take a brief detour from discussing initializations of inducing inputs to discuss the set of inducing variables that minimize the upper bounds in Lemmas 4 and 3.

Both the trace and the operator norm are unitarily invariant, so $\textup{K}_{M}$ is the optimal rank-$M$ approximation to $\textup{K}_{\textup{ff}}$ according to either of these norms.While the trace is not generally a matrix norm, it agrees with the norm $\|\cdot\|_{1}$ as $\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}}$ is SPSD. In particular, for any rank $M$ $N\times N$ SPSD matrix A satisfying $\textup{A}\prec\textup{K}_{\textup{ff}}$ (i.e. $\textup{K}_{\textup{ff}}-A$ is SPSD), $\textup{tr}(\textup{K}_{\textup{ff}}-\textup{K}_{M})\leq\textup{tr}(\textup{K}_{\textup{ff}}-\textup{A})$ and $\|\textup{K}_{\textup{ff}}-\textup{K}_{M}\|_{\textup{op}}\leq\|\textup{K}_{\textup{ff}}-\textup{A}\|_{\textup{op}}$ (see Horn and Johnson, 1990, Theorem 7.4.9.1).

As any subset of $M$ inducing variables will lead to a rank-$M$ matrix $\textup{Q}_{\textup{ff}}\prec\textup{K}_{\textup{ff}}$ this implies

Consider the inducing features defined as linear combinations of the random variables associated to evaluating the latent function at each observed input location, with weights coming from the eigenvectors of $\textup{K}_{\textup{ff}}$, i.e.

2.2 M-Determinantal point processes

We now return to the more practical case of using inducing points for sparse variational inference. In order to derive non-trivial upper bounds on $\textup{tr}(\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}})$ and $\|\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}}\|_{\textup{op}}$, we need a sufficiently good method for placing inducing points. When using differentiable kernel functions, many practitioners select the locations of the inducing points with gradient-based methods by maximizing the ELBO. As this is a high-dimensional, non-convex optimization algorithm, directly analyzing the result of this procedure is beyond our analysis.

In this section, we assume $M$ inducing points are subsampled from data according to an approximate M-determinantal point process ($M$-DPP) (Kulesza and Taskar, 2011) and use known bounds on the expected value of $\textup{tr}(\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}})$.The standard terminology is $k$-DPP. We use $M$ as this determines the number of inducing points and to avoid confusion with the kernel function. We note that if this scheme is used as an initialization prior to a gradient-based optimization of the evidence lower bound with respect to the inducing inputs, the resulting KL-divergence will be at least as small, so our bounds still apply after optimization of variational parameters.

Given an SPSD matrix L, an $M$-determinantal point process (Kulesza and Taskar, 2011) with kernel matrix L defines a discrete probability distribution over subsets of the $N$ columns of L, with positive probability only assigned to subsets of cardinality $M$. The probability of any subset of cardinality $M$ is proportional to the determinant of the principal submatrix formed by selecting those columns and the corresponding rows, that is for any set $Z$ of $M$ columns of $L$

where $\textup{L}_{Z,Z}$ is the principal submatrix of L with columns in $Z$. For a thorough introduction to determinantal point processes, as well as an implementation of many sampling methods, see Gautier et al. (2019).

The next important question to address is whether a $M$-DPP can be sampled with sufficiently low computational complexity for this to be a practical method for selecting inducing inputs. Naively computing the probability distribution over all $\binom{N}{M}$ subsets of size $M$ is prohibitively expensive. Kulesza and Taskar (2011) gave an algorithm that runs in polynomial time and yields exact samples from an $M$-DPP. Unfortunately, this algorithm involves computing an eigendecomposition of the $N\times N$ kernel matrix ($\textup{K}_{\textup{ff}}$ in our case), which is computationally prohibitive.

Recently, Dereziński et al. (2019) gave an algorithm for obtaining an exact sample from an $M$-DPP in time that is polynomial in $M$ and nearly-linear in $N$. However, the polynomial in $M$ is high. We instead consider an approximate algorithm and therefore derive the following simple corollary of Lemma 7.

where $\eta_{m}$ is the $m^{th}$ largest eigenvalue of $L$.

The corollary is completed by noting that for all $Z$, $\textup{tr}(\textup{L}-\textup{L}_{Z})\leq\textup{tr}(\textup{L})\leq Nv$.

Corollary 8 shows that sufficiently accurate approximate sampling from an $M$-DPP only has a small effect on the quality of the resulting $\textup{Q}_{\textup{ff}}$. High quality approximate samples can be drawn using a simple Markov Chain algorithm described in Anari et al. (2016), given as Algorithm 1. This MCMC algorithm is well-studied in the context of $M$-DPPs and their generalizations, and is known to be rapidly mixing (Anari et al., 2016; Hermon and Salez, 2019).

Let $\rho$ be an $M$-DPP with $N\times N$ kernel matrix $L$. Fix $\epsilon\in(0,1)$. Then Algorithm 1 produces a sample from a distribution $\rho^{\prime}$ satisfying

in not more than T(\epsilon)=2MN\mathopen{}\mathclose{{}\left(\log\log\mathopen{}\mathclose{{}\left(\frac{1}{\rho(Z_{0})}}\right)+\log\frac{2}{\epsilon^{2}}}\right) iterations, where $Z_{0}$ is the subset of columns at which the Markov chain is initialized.

Since the determinant of a matrix is equal to the product of the determinant of a principal submatrix times the determinant of the Schur complement of this submatrix, the greedy initialization used in Algorithm 1 is equivalent to starting with $U=\emptyset$ and iteratively adding $\operatorname*{\arg\!\max}_{x\in X}k(x,x)-\textup{k}_{f(x)u}\textup{K}_{\textup{uu}}^{-1}\textup{k}_{uf(x)}$ to $U$. This can be performed in time $\mathcal{O}(NM^{2})$, for example by computing the pivot rules of a rank-$M$ incomplete Cholesky decomposition of $\textup{K}_{\textup{ff}}$ (Chen et al., 2018, Algorithm 1).

The per iteration cost of Algorithm 1 is dominated by computing the acceptance ratio, which can be performed in $\mathcal{O}(M^{2})$, by iteratively updating a Cholesky or QR factorization of the matrix associated to the current set of columns. This makes the total cost of obtaining an $\epsilon$-approximate sample \mathcal{O}\mathopen{}\mathclose{{}\left(NM^{3}\log\log\mathopen{}\mathclose{{}\left(1/\rho(Z_{\textup{greedy}})}\right)+NM^{3}\log 2/\epsilon^{2}}\right), where $Z_{greedy}$ denotes the set of columns selected by greedily maximizing the determinant of the submatrix. Moreover, the subset selected by the algorithm is known to have a probability at least $1/(M!)^{2}$ of the maximum probability subset (Çivril and Magdon-Ismail, 2009; Anari et al., 2016). By using the fact that the the maximum probability subset is more probable than the uniformly distributed probability, we obtain

We now take a brief detour to consider a different approach to initializing inducing inputs before completing the proof of a priori bounds on the KL-divergence.

2.3 Ridge Leverage Scores

Computing the ridge leverage scores exactly is too computationally expensive, as it involves inverting the kernel matrix. However, practical approximate versions of leverage sampling algorithms that retain strong theoretical guarantees have been developed.

Ridge leverage based sampling algorithms select a subset of training data to use as inducing points. Each point is sampled independently into the subset with probability proportional to its leverage score. Approximate versions of this algorithm generally rely on overestimating the ridge leverage scores, which lead to equally strong accuracy guarantees compared to using the exact ridge leverage scores, at the cost of sampling more points in the approximation.

We consider the application of Algorithm 3 in Musco and Musco (2017) to the problem of selecting inducing inputs for sparse variational inference in GP models. This algorithm comes with the following bounds on the quality of the resulting Nyström approximation.

While in Section 4.2.2 $M$ was fixed and the quality of the resulting approximation was random, in the algorithm discussed above $\mathbf{M}$ is additionally random.

An alternative approach to sampling using ridge leverage scores specifies a desired level of accuracy of the resulting approximation, and the number of points selected is chosen to obtain this approximation quality with fixed probability. This has the advantage of not requiring the user to manually select the number of inducing points, but may lead to a number of inducing points being used that exceeds a practical computational budget. We discuss the application of this approach to variational Gaussian process regression in Section G.

3 A-Priori Bounds on the KL-divergence

In the previous sections, the results on the quality of approximation depended on the eigenvalues of $\textup{K}_{\textup{ff}}$. As these eigenvalues depend on the covariates $X$, and we would like to make statements that apply to a wide-range of data sets, we assume $X$ is a realization of a random variable $\mathbf{X}$, and make assumptions about the distribution of $\mathbf{X}$.

If each $\mathbf{x}\in\mathbf{X}$ is i.i.d. distributed, according to some measure with continuous density $p(x)$, in the limit as the amount of data tends to infinity, the matrix $\frac{1}{N}\mathbf{K}_{\bf ff}$ behaves like the operator $\mathcal{K}$ (Koltchinskii and Giné, 2000) defined with respect to this $p$. For finite sample sizes, the large eigenvalues of $\frac{1}{N}\mathbf{K}_{\bf ff}$ tend to overestimate the corresponding eigenvalues of $\mathcal{K}$ and the small eigenvalues of $\frac{1}{N}\mathbf{K}_{\bf ff}$ tend to underestimate the small eigenvalues of $\mathcal{K}$. We make this precise through a minor generalization of a lemma of Shawe-Taylor et al. (2005).

where $\bar{c}=\frac{1}{N}\sum_{n=1}^{N}c_{n}$.

Consider the rank-$M$ approximation to $\mathbf{K}_{\bf ff}$ given by truncating this Mercer expansion, $[\bm{\Phi}]_{i,j}=\sum_{m=1}^{M}\lambda_{m}\phi_{m}(\mathbf{x}_{i})\phi_{m}(\mathbf{x}_{j})$. Then $[\mathbf{K}_{\bf ff}-\bm{\Phi}]_{i,j}=\sum_{m=M+1}^{\infty}\lambda_{m}\phi_{m}(\mathbf{x}_{i})\phi_{m}(\mathbf{x}_{j})$, so $\mathbf{K}_{\bf ff}-\bm{\Phi}\succ 0$.

For any covariates $\{\mathbf{x}_{n}\}_{n=1}^{N}$ satisfying the conditions of the lemma,

Taking expectations on both sides with respect to the covariate distribution,

The interchanging of integral and sum is justified by Fubini’s theorem as each $\phi_{m}$ is square integrable, each eigenvalue is non-negative, and the sum converges by Mercer’s theorem. We used the non-negativity of $\phi_{m}(x)^{2}$ in the second inequality to bound the expectation of $\phi_{m}(x)^{2}$ under $p_{n}$ in terms of its expectation under $q$.

Suppose the covariate distribution has identically distributed marginals, each with density $p(x)$, then

where $\lambda_{m}$ is the $m^{th}$ largest eigenvalue of the operator associated to the kernel and the distribution with continuous density $p(x)$.

This corollary follows from Lemma 11 by taking $q=p$ and $c_{n}=1$ for all $n$. For simplicity, we will state our main results using the assumptions of this corollary, though the generalization to cases with non-identical marginals satisfying the conditions of Lemma 11 is immediate. We have now accumulated the necessary preliminaries to prove our main theorems.

where the expectation is taken over the covariates, the mechanism for initializing inducing points and the observations.

Finally, taking expectation with respect to the covariate distribution over the covariate distribution and applying Lemma 11,

With the same assumptions on the covariates and inducing point distributions as in Theorem 13, but with the assumption that $\mathbf{y}|\mathbf{X}$ is conditionally Gaussian distributed with mean zero and covariance matrix $\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I}$,

where the expectation is taken over the covariate distribution, the observation distribution and the initialization mechanism.

The proof of Theorem 14 is nearly identical to the proof of Theorem 13, applying Lemma 4 instead of Lemma 3 in the first line.

Under the assumptions of Theorem 13, with probability at least $1-\delta$,

Under the assumptions of Theorem 14, with probability at least $1-\delta$,

3.2 Bounds for Ridge Leverage Score Sampling

We now state and derive statements similar to Corollaries 15 and 16 for a ridge leverage score initialization utilizing Musco and Musco (2017, Algorithm 3). In order to this we us that for any SPSD $A$, $\textup{tr}(A)\leq N\|A\|_{\textup{op}}$, so that Lemma 3 implies

Combining Lemmas 10 and 12 and using Markov’s inequality twice with Eq. 22 or Eq. 23 and a union bound respectively leads to the following bounds on the performance of sparse inference using ridge leverage scores:

when inducing points are initialized using Musco and Musco (2017, Algorithm 3).

when inducing points are initialized using Musco and Musco (2017, Algorithm 3).

3.3 Are these bounds useful?

Having established probabilistic upper bounds on the KL-divergence resulting from sparse approximation, a simple question is whether these bounds offer any insight into the efficacy of sparse inference. If in order for the upper bounds to be small, we need to take $M=N$, then they would not be useful, as it is already known that by taking $Z=X$, exact inference is recovered. In the next section, we discuss bounds on the eigenvalues of $\mathcal{K}$ for common kernels and input distribution. These bounds show that for many inference problems, the upper bounds in Theorems 13, 14, 17 and 18 imply that the KL-divergence can be made small with $M\ll N$ inducing points.

Bounds for Specific Kernels and Covariate Distributions

In this section, we consider specific covariate distributions and commonly used kernels, and investigate the implications of the upper bounds derived in Section 4. These results are summarized in Footnote 5. We begin with the case of the popular squared exponential kernel and Gaussian covariates in one-dimension. This kernel and covariate distribution are one of the few instances in which the eigenvalues of $\mathcal{K}$ have a simple analytic form. In Section 5.1.1, we consider the analogous multi-dimensional problem. In Section 5.2 we discuss implications for stationary kernels with compactly supported inputs, including the well-studied Matérn kernels.

and one-dimensional covariates distributed according to $\mathcal{N}(0,\beta^{2})$, the eigenvalues of $\mathcal{K}$ are (Zhu et al., 1997)

Let $k$ be a squared exponential kernel. Suppose that $N$ real-valued (one-dimensional) covariates are observed, with identical Gaussian marginal distributions. Suppose the conditions of Theorem 13 are satisfied for some $R>0$. Fix any $\gamma\in(0,1]$. Then there exists an $M=\mathcal{O}(\log(N^{3}/\gamma))$ and an $\epsilon=\Theta(\gamma/N^{2})$ such if inducing points are distributed according to an $\epsilon$-approximate $M$-DPP with kernel matrix $\textup{K}_{\textup{ff}}$,

Similarly, for any $\delta\in(0,1/32)$ using the ridge leverage algorithm of Musco and Musco (2017) and choosing $S$ appropriately, with probability $1-5\delta$, \mathbf{M}=\mathcal{O}\mathopen{}\mathclose{{}\left(\log\frac{N^{2}}{\delta^{2}\gamma}\log\frac{\log(N^{2}/\delta^{2}\gamma)}{\delta}}\right) and

The implicit constants depend on the kernel hyperparameters, the likelihood variance, the variance of the covariate distribution and $R$.

If we consider $\gamma$ and $\delta$ as fixed constants (independent of $N$), this implies that if inducing points are placed using an approximate $M$-DPP we can choose $M=\mathcal{O}(\log(N))$ inducing points leading to a computational cost of $\mathcal{O}(N(\log N)^{4})$ while for approximate ridge leverage scores sampling $\mathcal{O}(\log N\log\log N)$ inducing points suffice leading to a cost at most $\mathcal{O}(N(\log N)^{2}(\log\log N)^{2})$.

The generalization of Corollary 19 to the case of multi-dimensional input distributions is relatively straightforward. The multi-dimensional version of the squared exponential kernel can be written as a product of one dimensional kernels, i.e.

For any kernel that can be expressed as a product of one-dimensional kernels, and for any covariate distribution that is a product of one-dimensional covariate distributions, the eigenvalues of the multi-dimensional covariance operator is the product of the one-dimensional analogues. When obtaining rates of convergence, we lose no generality in assuming that the kernel is isotropic as is the covariate distribution. Otherwise, consider the direction with the shortest lengthscale, and the covariate distribution with the largest standard deviation and the eigenvalues of this operator are larger than a constant multiple of the corresponding eigenvalues of the non-isotropic operator.

In the isotropic case, each eigenvalue is of the form,

for some integer $m^{\prime}$ with $a,A$ and $B$ defined as in the one-dimensional case. Note that $m$ and $m^{\prime}$ are no longer equal. The number of times each eigenvalue with $m^{\prime}$ in the exponent is repeated is equal to the number of ways to write $m^{\prime}$ as a sum of $D$ non-negative integers. By counting the multiplicity of each eigenvalue, Seeger et al. (2008) arrived at the bound

In order to prove a multi-dimensional analogue of Corollary 19 we need an upper bound on $\sum_{m=M+1}^{\infty}\lambda_{m}$. This can be derived with following an argument made by Seeger et al. (2008, Appendix II).

The proof of Proposition 21 is in Section D.1.

The implicit constant depends on the kernel hyperparameters, the variance matrix of the covariate distribution, $D$ and $R$. With the same assumptions but applying the RLS algorithm of Musco and Musco (2017) to selecting inducing inputs, for any $\delta\in(0,1/32)$ there exists a choice of $S$ such that with probability $1-5\delta$, \mathbf{M}=\mathcal{O}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(\log\frac{N^{2}}{\delta\gamma}}\right)^{D}(\log\log\frac{N^{2}}{\delta\gamma}+\log(1/\delta)}\right) and

The proof follows from Proposition 21 and Theorem 13 or Theorem 17, by choosing parameters appropriately.

If we allow the implicit constant to depend on $\gamma$ and $\delta$, this implies that for inducing ploints distributed accoding to an approximate $M$-DPP we can choose $M=\mathcal{O}((\log N)^{D})$ inducing points leading to a computational cost of $\mathcal{O}(N(\log N)^{3D+1})$ while for approximate ridge leverage scores sampling $\mathcal{O}((\log N)^{D}\log\log N)$ inducing points suffice leading to a $\mathcal{O}(N(\log N)^{2D}(\log\log N)^{2})$ computational cost.

In order for the KL-divergence to be less than a fixed constant, the exponential scaling of the number of inducing points in the dimensions of the covariates is inevitable, as we will show in Section 6. However, practically the situation may not be quite so dire. First, many practioners use a SE-ARD kernel. If the data is essentially constant over many dimensions, then when training with empirical Bayes, the lengthscales of these dimensions tends to become large, effectively reducing the dimensionality of the inference problem. Additionally, in the case when covariates fall on a smooth, low-dimensional manifold, the decay of the eigenvalues only depends on the dimensionality and smoothness properties of this manifold, see Altschuler et al. (2019, Theorem 4). In addition, for a given problem, the dimensionality $D$ is fixed, meaning that the dependence of the number of inducing points $M$ depends polylogarithmically on $N$. This growth is slower than any polynomial, i.e. $(\log N)^{D}=o(N^{\epsilon})$ for $\epsilon>0$.

We also note that Corollary 22 can easily be adapted using Lemma 11 to show that if all of the $\mathbf{x}_{n}$ are drawn from any compactly supported distributions with continuous densities that are all bounded by some universal constant, the same asymptotic bound on the number of inducing points applies. This follows from noting that under these assumptions, $p_{n}(x)$, satisfies $p_{n}(x)<cq(x)$ where $q(x)$ is a Gaussian density for some $c>0$, so we can apply Lemma 11 to bound the expectation of the sum of the matrix eigenvalues associated to $p_{n}$ in terms of the eigenvalues associated to $q$.

2 Compactly Supported Inputs and Stationary Kernels

Stationary, continuous kernels can be characterized through Bochner’s theorem, which states that any such kernel is the Fourier transform of a positive measure, i.e.

We will refer to $s(\omega)$ as the spectral density of $k$.We assume $\kappa(x-x^{\prime})$ decays sufficiently rapidly so that such a continuous spectral density exists. The decay of the spectral density conveys information about how smooth the kernel function is.

Widom’s theorem (Widom, 1963) relates the decay of the eigenvalues of $\mathcal{K}$ to the decay of $s$. Widom’s theorem applies to input distributions with compact support and stationary kernels with spectral density satisfying several regularit conditions (stated in Section D). Seeger et al. (2008) give a corollary of Widom’s theorem, which is sufficient in many instances to obtain bounds on the number of inducing points needed for Theorems 14 and 13 to converge.

Let $k$ be an isotropic kernel (i.e. $\kappa(\alpha)=\kappa(\alpha^{\prime})$ if $\|\alpha\|=\|\alpha^{\prime}\|$). Suppose $k$ satisfies the criteria of Widom’s theorem, the covariate distribution has density zero outside a ball of radius $T$ around the origin, and is bounded above by $\tau$, then

Matérn kernels are widely applied to problems where the data generating process is believed to lead to non-smooth functions, and are known to satisfy the conditions of Widom’s theorem (Seeger et al., 2008). These kernels are defined as (Rasmussen and Williams, 2006),

where $K_{\nu}$ is a modified Bessel function. The spectral density of the Matérn kernel is

Lemma 24 tells us that for compactly supported covariates with bounded density and the Matérn kernel with smoothness paramater $\nu$

It follows that $\sum_{m=M+1}^{D}\lambda_{m}=\mathcal{O}(M^{\frac{-2\nu}{D}})$. From this, we can derive a result of the same form as Corollary 22 for Matérn kernels and compactly supported input distributions.

Under the same assumptions if inducing points are initialized using the RLS algorithm of Musco and Musco (2017) with $\delta\in(0,1/32)$ there exists an S=\mathcal{O}\mathopen{}\mathclose{{}\left(N^{\frac{2D}{2\nu+D}}(\gamma\delta^{2})^{\frac{-D}{2\nu+D}}}\right) such that with probability at least $1-5\delta$,

and $\mathbf{M}\leq S\log\frac{S}{\delta}$.

If we consider $\gamma$ and $\delta$ as fixed constants (independent of $N$), this implies that for an initialization with $M$-DPP we can choose $M=\mathcal{O}(N^{\frac{2D}{2\nu-D}})$ inducing points leading to a computational cost of $\mathcal{O}(N^{\frac{2\nu+5D}{2\nu-D}}\log(N))$ while for approximate ridge leverage scores sampling $\mathcal{O}(N^{\frac{2D}{2\nu+D}}\log N)$ inducing points suffice leading to a cost at most $\mathcal{O}(N^{\frac{2\nu+5D}{2\nu+D}}(\log N)^{2})$.

The first part corollary follows from Theorem 13, noting that we need to choose $M$ such that

for some constants $C,C^{\prime}$. The second part follows from similar considerations applied to Theorem 17.

These bounds on $M$ are vacuous (i.e. are no smaller than $M=\mathcal{O}(N)$) for Matérn kernels in high dimensional spaces or with low smoothness parameters. Additionally, the cost of sampling the $M$-DPP using Algorithm 1 makes this inference scheme less expensive than exact GP inference only when $\nu>2D$. If we instead make the stronger assumptions required by Theorem 14, we can choose $M=\mathcal{O}(N^{\frac{D}{2\nu-D}})$, which implies a computational complexity less than exact GP regression if $\nu>\frac{5}{4}D$.

The bounds for the RLS initialization are generally sharper, and are non-vacuous for all $\nu>D/2$ with the weaker assumptions on $\mathbf{y}$. Additionally, the computational complexity of choosing inducing points using the RLS algorithm is the same as the cost of inference up to logarithmic factors, so that for $\nu>D/2$ the cost of sparse inference with the RLS initialization is (asymptotically) smaller than the cubic cost of exact GP regression.

Lower Bounds on the Number of Inducing Points Needed

In Sections 4 and 5, we showed that for many problems the number of inducing points can grow sub-linearly with the number of data points, while maintaining a small KL-divergence between the approximate and exact posteriors. In this section we consider the inverse question, i.e. how many inducing points are necessary to avoid having the KL-divergence grow as the amount of data increases? In this section, we prove a-priori lower bounds on the KL-divergence under similar assumptions to those used in proving the upper bounds in Section 4.

Naively, it appears that the lower bound in Lemma 4 gives us a starting place for a lower bound on the KL-divergence. From this bound,

While lower bounding this quantity can be done using the approach taken in this section, it is not the most interesting quantity to study, as we average over $\mathbf{y}$ conditioned on $\mathbf{X}$ and $\mathbf{Z}$. This would not give a valid lower bound if the locations of the inducing points depend on $\mathbf{y}$, as illustrated in Fig. 1. This approach would establish a lower bound for initialization schemes considered in the previous sections, as well as any initialization scheme that does not take the observed $y$ into account, but not the common practice of performing gradient ascent on the ELBO with respect to inducing inputs.

In this section, we establish a lower bound on the number of inducing variables needed for the KL-divergence not to become large, which is valid regardless of the method for selecting inducing variables or the distribution of $\mathbf{y}$. These bounds assume that the covariates are independent and identically distributed (in contrast to the upper bounds, which do not require independence and require a slightly weaker condition than identical marginals). The independence assumption is necessary in order to lower bound the eigenvalues of the covariance matrix. For example, if all of the covariates were identically distributed and equal, the covariance matrix would be rank-1 and so a single inducing point could be used regardless of the size of the data set.

The proof of the lower bounds proceeds in two parts:

Second, we use a result on the concentration of eigenvalues of the kernel matrix to those of the corresponding operator due to Braun (2006) to derive a lower bound that holds with fixed probability under the assumption that the covariates are independent and identically distributed.

In the case of SE-kernel and Gaussian covariates, we establish a lower bound with the same dependence on $N$ as our upper bounds, that is we need $M=\Omega((\log N)^{D})$. In the case of Matérn kernels with uniform covariates and $\nu>1$, we establish a lower bound that increases as a power of $N$. However, there is a large gap between our upper and lower bounds for Matérn kernels, indicating room for improvement. These results are summarized in Table 2. While our results are stated in terms of inducing points, they hold for more general inducing variables.

In this section, we derive a lower bound on the KL-divergence that holds for any $y$ and $Z$ and depends on $\mathbf{X}$.

Given a kernel $k$, likelihood model with variance $\sigma^{2}$ and random covariates $\mathbf{X}$. Then,

In order to establish lower bounds on the number of inducing variables needed to ensure the KL-divergence does not grow as a function of $N$, it suffices to analyze the behavior of the lower bound in Lemma 27 for random covariates as a function of both $M$ and $N$.

2 Structure of the Argument

$N\gamma_{N}\lambda_{M+1}$ tends to infinity as $N$ tends to infinity,

it must be the case that the KL-divergence tends to infinity as a function of $N$ (at a rate $\Omega(N\gamma_{N}\lambda_{M+1})$).

3 Concentration of Eigenvalues

In order to complete the argument in the previous section, we need a more fine-grained understanding of the behavior of eigenvalues of $\mathbf{K}_{\bf ff}$ than given in Lemma 11. For this, we rely on the following result:

For the remainder of this section, we consider specific cases of kernel and input distributions for which we know properties of the spectrum, and derive lower bounds on the number of features needed so that the KL-divergence is not an increasing function of $N$.

4 Squared Exponential Kernel and Gaussian Covariates

For any fixed $\delta\in(0,1)$ the first term on the right hand side tends to zero with $N$ since $\eta<1$. For any $M\leq\log_{B}(\sqrt{A/(2av^{2})}N^{-\eta}\sqrt{\delta})$,

The second term is less than $1/2$ and the last term tends to for large $N$. We conclude that for any such $M$, the KL-divergence is bounded below by $c_{N}\frac{N\lambda_{M+1}}{8\sigma^{2}}=\Omega(N^{1-\eta})$, where $\lim_{N\to\infty}c_{N}=1$.

For univariate Gaussian kernels, if we choose $\eta=.01$ in the above argument, we get that for any $M\leq\frac{1}{\log(1/B)}\log(N^{.01}\sqrt{\frac{2av^{2}}{\delta A}})=\Omega(\log N)$, the KL-divergence is $\Omega(N^{.99})$, and will therefore be large as $N$ increases. We therefore need $M$ to grow faster than this to avoid this if we want the KL-divergence to be small for large $N$.

5 The Isotropic SE-kernel and Multidimensional Gaussian Covariates

In order to obtain lower bounds in the multivariate case, we first obtain a lower bound on the individual eigenvalues of the operator $\mathcal{K}$.

The proof (Section E) relies on a counting argument and standard bounds on binomial coefficients. We can now combine Lemmas 28, 21 and 30 in order to bound the multivariate SE-kernel with Gaussian inputs.

The proof follows from Lemma 27 by choosing $r$ appropriately in Lemma 28 to bound the empirical eigenvalues and using Propositions 21 and 30 to bound eigenvalues in the appropriate directions to control the error term. Details are given in Section E.

6 Lower Bounds for Kernels with Polynomial Decay

As discussed in Section 5 some popular choices of kernels lead to eigenvalues that decay polynomially instead of exponentially. For example, Widom (1963, Theorem 2.1) implies that the eigenvalues of the operator associated to the Matérn kernel with smoothness parameter $\nu$ and covariates uniformly distributed in the unit cube has eigenvalues satisfying $C_{1}m^{\frac{-2\nu+D}{D}}\leq\lambda_{m}\leq C_{2}m^{\frac{-2\nu+D}{D}}$ for some constant $C_{1}$ and $C_{2}$ independent of $m$ i.e. $\lambda_{m}=\Theta(m^{\frac{-2\nu+D}{D}})$.See Seeger et al. (2008) for more details on the derivation of this from Widom’s Theorem.

Proof We have $\lambda_{r}=\Theta(r^{-\eta})$ so $\sum_{s=r}^{\infty}\lambda_{s}=\Theta(r^{1-\eta})$. Choose $r=N^{\gamma}$ for some $\gamma\in(0,1)$ and $M+1=N^{\zeta}$. In this case the error term in Lemma 28 becomes:

Following the earlier proof sketch, we must show the RHS tends to a value less than $1$. To ensure the first of the three summands is small, we choose $\gamma\in(0,\frac{1}{4+\eta})$. Given this choice, for large $N$ the third summand in the error term is always smaller than the second, so that this entire term is $o(1)$ given the supposition that $\zeta\leq\gamma\frac{\eta-1}{\eta}$. We conclude that if $M=N^{-\zeta}$ for $\zeta\in(0,\gamma\frac{\eta-1}{\eta})$ with probability $1-\delta$, the lower bound in Lemma 27 is at least $N\lambda_{M+1}=\Omega(N^{1-\zeta\eta})$. In the case of $D$-dimensional Matérn kernels and a uniform covariate distribution ($\eta=\frac{2\nu+D}{D}$), by choosing $\zeta$ as large as possible, this means that for an arbitrary $\epsilon>0$, the KL-divergence is lower bounded by an increasing function of $N$ if fewer than \Omega\mathopen{}\mathclose{{}\left(N^{\frac{2\nu D}{(2\nu+5D)(2\nu+D)}-\epsilon}}\right) inducing variables are used. This lower bound on the number of inducing variables becomes vacuous (i.e. the exponent tends to ) as $\eta\to 1$ from above, meaning it is not useful when applied to many kernels that we expect would be very difficult to approximate. There is a large gap between the upper and lower bound, particularly when $\eta$ is near $1$ (i.e. for non-smooth kernels). The gap between the bounds is in part introduced by needing to choose $M$ so that the error term from Lemma 28 remains lower order. If we heuristically allow ourselves to replace matrix eigenvalues with the corresponding scaled operator eigenvalues and neglect the error term, we obtain a lower bound of $\Omega(N^{\frac{1}{\eta}})$, bringing the lower bound more closely in line with the upper bound of $\mathcal{O}(N^{\frac{1}{\eta-1}})$. The remaining gap between these bounds is essentially due to only bounding a single eigenvalue in the lower bound, while bounding the sum of eigenvalues in the upper bound. Improving the analysis to close the gap between the upper and lower bounds is important for better understanding the efficacy of sparse methods with non-smooth kernels.

Practical Considerations

Up to this point, we proved statements about the asymptotic scaling properties of variational sparse inference. Our results indicated which models could be well-approximated with relatively few inducing points for sufficiently large data sets. In this section, we investigate the limitations and practical implications of our results to real situations with finite amounts of data. We consider the applicability of our results to practical implementations, and perform empirical analyzes on how marginal likelihood bounds converge. Additionally, our proof suggests a specific procedure for choosing inducing points that differs from methods that are currently commonly applied. We empirically investigate this procedure, and provide recommendations on how to initialize inducing points.

Any practical implementation of a Gaussian process method will be influenced by the finite precision with which floating-point numbers are represented in a computer. These issues are not explicitly addressed in our mathematical analysis, which assume calculations are in exact arithmetic. Here, we briefly discuss the effects of this finite precision on 1) the implementation, 2) the precision to which we can expect convergence in practice compared to our analysis, and 3) the way that this is quantified by marginal likelihood bounds.

Finding various quantities for Gaussian process regression requires computing log determinants and matrix inverses. When the smallest and largest eigenvalues of the kernel matrix are many orders of magnitude apart, these computations become ill-conditioned, meaning that small changes on the input can lead to large changes to the output. For example, tiny changes in the elements of the vector $\mathbf{f}_{X}$ can lead to huge variations in the vector $\textup{K}_{\textup{ff}}^{-1}\mathbf{f}_{X}$ when $\textup{K}_{\textup{ff}}$ has an eigenvalue close to zero (see Deisenroth et al., 2019, §6.2 for a visual illustration). This typically occurs when considering many highly-correlated inputs to the GP (e.g clusters of nearby points with similar input values). These points have a high probability of having very similar function outputs under the prior. This ill-conditioning arises naturally in GPs when considering e.g. evaluating the prior density on function values: small differences in the function values result in huge changes to the value of the probability density. If the sensitivity of the calculations becomes too large, then the finite precision with which numbers are represented can lead to considerable error.

In the variational methods we consider, determinants and inverses are found based on the Cholesky decomposition of the kernel matrix: $\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I}$ for exact implementations, and $\textup{K}_{\textup{uu}}$ for sparse approximations.Conjugate gradient and Lanczos methods also give exact answers when they are run for sufficient iterations, and have been successfully applied in practice (Gibbs and Mackay, 1997; Davies, 2015; Gardner et al., 2018). When faced with a problem that is too ill-conditioned, most Cholesky implementations terminate with an exception. This can be seen as desirable from the point of view that a successful run usually indicates an accurate result.

Even in cases when the data set can be well-described by a GP model with hyperparameters that lead to reasonably well-conditioned matrices, conditioning problems frequently arise during training, when the log marginal likelihood or ELBO is values for other candidate hyperparameter values. For example, for stationary kernels, large lengthscales contribute to conditioning problems, as they increase the correlation between distant points. Hyperparameters are typically found by (approximately) maximizing the log marginal likelihood (Eqs. 6 and 9). Since these objective and their derivatives can be evaluated in closed-form, fast-converging quasi-Newton methods such as (L-)BFGS are commonly used. These methods often propose large steps, which lead to the evaluation of hyperparameter settings where the Cholesky decomposition raises an exception. Even though these hyperparameter settings are often of poor quality, (L-)BFGS still requires an evaluation of the objective function to continue the search. The Cholesky errors must therefore be avoided to successfully complete the entire optimization procedure.

1.2 Improving Matrix Conditioning

Increasing the smallest eigenvalue of the kernel matrix improves the conditioning. In exact implementations this can be done by increasing the likelihood noise variance, as we need to decompose $\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I}$ which has eigenvalues that are lower-bounded by $\sigma^{2}$This can be done by reparameterizing the noise to have a lower bound.. On the other hand, the sparse variational approximation requires inverting $\textup{K}_{\textup{uu}}$ without any noise. However, it is important to note that the conditioning of $\textup{K}_{\textup{uu}}$ is better than $\textup{K}_{\textup{ff}}$ for two reasons. Firstly, it is a smaller matrix, and often issues of conditioning are less severe for smaller matrices. Secondly, if inducing points are selected using a method that introduces negative correlations clusters of highly-correlated points are unlikely to appear in $\textup{K}_{\textup{uu}}$. Nevertheless, it is still possible for the Cholesky decomposition to fail, particularly when trying different hyperparameter settings when maximizing the ELBO.

To improve robustness in the sparse approximation, a small diagonal “jitter” matrix $\epsilon\textup{I}$, with $\epsilon$ commonly around $10^{-6}$, is added to $\textup{K}_{\textup{uu}}$, introducing a lower bound the on its eigenvalues. This change is often enough to avoid decomposition errors during optimization. While this modification changes the problem that is solved, the effect is typically small. Some software packages (e.g. GPy, since 2012) increase jitter adaptively by catching exceptions inside the optimization loop to only introduce bias where it is necessary.

1.3 Quantifying the Effect of Jitter

Let $\mathcal{L}_{\epsilon}$ denote the evidence lower bound computed with jitter $\epsilon\geq 0$ added to $\textup{K}_{\textup{uu}}$, that is

Then $\mathcal{L}_{\epsilon}$ is monotonically decreasing in $\epsilon$. Similarly if $\mathcal{U}_{\epsilon}$ denotes the upper bound Eq. 12 computed with added jitter to $\textup{K}_{\textup{uu}}$, that is

Then $\mathcal{U}_{\epsilon}$ is monotonically increasing in $\epsilon$. In particular, adding jitter can only make the upper bound on the log marginal likelihood larger and the ELBO smaller.

The proof (Section F) is a consequence of $\textup{Q}_{\textup{ff}}(\epsilon)+\sigma^{2}\textup{I}\succ\textup{Q}_{\textup{ff}}(\epsilon^{\prime})+\sigma^{2}\textup{I}$ for $\epsilon^{\prime}>\epsilon\geq 0$. Proposition 34 shows that even with jitter the upper and lower bounds are still valid. However, they are not exactly equal, even when $M=N$, due to the additional gap caused by the jitter. From a practical point of view, the impact is typically very small, with a gap being introduced on the order of a few nats.

To summarize, we saw 1) that jitter was needed to stabilize the computation of the hyperparameter objective functions using standard implementations of Cholesky decomposition, and 2) that jitter and finite floating-point precision prevented the approximate posterior and bounds from converging to their exact values. Notably, we can quantify the effect of the finite precision calculations using the same bounds as what is used to determine the effect of using an approximate inducing point posterior (Proposition 34). The variational bounds we analyzed in this work therefore provide a unified way of measuring the effect of both exact arithmetic approximate posteriors and the impact of finite precision arithmetic on the quality of the approximation.

2 Placement of Inducing Inputs

When training a sparse Gaussian process regression (Titsias, 2009b) model, we need to select the kernel hyperparameters as well as the inducing inputs, with the hyperparameters determining the generalization characteristics of the model, and the inducing inputs the quality of the sparse approximation. In the time since Snelson and Ghahramani (2006) and Titsias (2009b) introduced joint objective functions for all parameters, it has become commonplace to find the final set of inducing variables by optimizing the objective function together with the hyperparameters. Because this makes the inducing input initialization procedure less critical for final performance, less attention has been placed on it in recent years than in e.g. the kernel ridge regression literature.Kernel ridge regression lacks a joint objective function for the approximation and hyperparameters. Hyperparameters are commonly selected through cross-validation. However, the number of optimization parameters added by the inducing inputs is often large, and convergence can be slow, which makes the optimization cumbersome.

We set the free parameters for each of the methods as follows. For K-means, we run the Scipy implementation of K-means++ with $M$ centres. Gradient-based optimization is initialized using greedy variance selection. This choice was made since it was found to perform better than uniform selection, and our goal is to quantify how much can be gained by doing gradient-based optimization, and whether it is worth the cost. We ran $10^{4}$ steps of L-BFGS, at which point any improvement was negligible compared to adding more inducing variables. Approximate $M$-DPP sampling was done following Algorithm 1, using $10^{4}$ iterations of MCMC. For RLS, we use an adaptation of the public implementation of Musco and Musco (2017, Algorithm 3), which omits many of the constants derived in the proofs, and therefore loses theoretical guarantees.Their implementation is available at: https://github.com/cnmusco/recursive-nystrom. We additionally modify the algorithm to ensure that it selects exactly $M$ inducing points.

We consider 3 data sets from the UCI repository that are commonly used in benchmarking regression algorithms, “Naval” ($N_{train}=10740,N_{test}=1194,D=14$) , “Elevators” ($N_{train}=14939,N_{test}=1660,D=18$) and “Energy” ($N_{train}=691,N_{test}=77,D=8$). These data sets were chosen as near-exact sparse approximations could be found, so convergence could be illustrated.Not all data sets exhibit this property. For instance, the “kin40k” data set still isn’t near convergence when $M=\frac{N}{2}$ due to very short optimal lengthscales. A step functions being present would cause this, and would indicate that squared exponential kernels are inappropriate. Naval is the result of a physical simulation and the observations are essentially noiseless. To make statistical estimation more difficult, we add independent Gaussian noise with standard deviation $0.0068$ to each observation. For all experiments, we use a squared exponential kernel with automatic relevance determination (ARD), i.e. a separate lengthscale per input dimension.

We first consider regression with fixed hyperparameters to illustrate convergence in a situation that is directly comparable to our theoretical results. We investigate which of the inducing point selection methods recovers the exact model with the fewest inducing points. The hyperparameters are set to the optimal values for an exact GP model, or for “Naval” a sparse GP with 1000 inducing points. We find the hyperparameters by maximizing the exact GP log marginal likelihood using L-BFGS. This setting is for illustrative purposes only, as computing exact log marginal likelihood is not feasible in practical situations where sparse methods are of actual interest. In the next section we consider hyperparameters that are learned using the ELBO (Eq. 9).

Figure 5 shows the performance of various methods of selecting inducing points as we vary $M$, as measured by the evidence lower bound, test root mean squared error and per data point test negative log predictive density. From the results, we can observe the following:

For very sparse models where the ELBO is considerably lower than the true marginal likelihood, gradient-based tuning of the inducing inputs consistently performs best in all metrics.

The benefit of gradient-based tuning is small when many inducing points are added, provided they are added in the good locations. Greedy variance selection and $M$-DPP find these good locations, as they consistently recover the true GP’s performance with only a small number of additional inducing variables.

K-means, uniform subsampling, and RLS tend to underperform, and require far more inducing variables to converge to the exact solution. In our experiment, they never converge quicker than greedy variance selection.

In terms of the upper bound, greedy variance selection and $M$-DPP sampling both provide the best results.

Greedy variance selection seems to provide all the desirable properties in this case: convergence to the exact results with few inducing variables, simple to implement, and fast since it does not require as many expensive operations as optimization or sampling. The approximate ridge leverage score algorithm is also reasonably fast, and perhaps careful tuning of hyperparameters or different algorithms for approximate ridge leverage scores could lead to improved performance in practice.

2.2 Training procedure and hyperparameter optimization

Our results imply that for large enough data sets, an approximation with high sparsity can be found for the optimal hyperparameter setting that has a small KL-divergence to the posterior. This implies that the bias in the hyperparameter selection also is likely to be small. An impediment to finding the high-quality approximation for the optimal hyperparameters, is that our results depends on the inducing inputs being chosen based on properties of the kernel with the same hyperparameters that are used for inference. Here, we investigate several procedures for jointly choosing the hyperparameters and inducing inputs, in the regime where enough inducing variables are used to recover a close to exact model.

We propose a new procedure based on the greedy variance selection discussed in the previous section. To account for the changing hyperparameters, we alternately optimize the hyperparameters, and reinitialize the inducing inputs with greedy variance selection using the updated hyperparameters. This avoids the high-dimensional non-convex optimization of the inducing inputs, while still being able to tailor the inducing inputs to the kernel. In effect, the method behaves a bit like variational Expectation-Maximization (EM) (Beal and Ghahramani, 2003), with the inducing input selection taking the place of finding the posterior. When enough inducing points are used, the reinitialization is good enough to make the ELBO almost tight for the current setting of hyperparameters. We terminate when the reinitialization does not improve the ELBO. We note that reinitialization would not benefit K-means or uniform initializations (beyond random chance), as the inducing points that are selected do not depend on the setting of the kernel hyperparameters.

We start our evaluation by running all methods from the previous section in addition to greedy variance selection with reinitialization (Fig. 6). The initial inducing inputs are set with the untrained initialized hyperparameters, after which the hyperparameters are maximized w.r.t. ELBO (Eq. 9) using L-BFGS, with the reinitialization being applied for “Greedy variance (reinit.)”. We observe that the reinitialized greedy variance method provides consistent fast convergence to the exact model.

To evaluate the benefit of gradient-based optimization, we compare it to the reinitialized greedy variance method (the best from Fig. 6), as well as K-means. For the initial setting of the inducing inputs when optimizing inducing inputs, we use the greedy variance selection (denoted “gradient”). Since Fig. 6 shows that optimization of the inducing inputs is not needed to converge to the exact solution, the question becomes whether it is faster to perform gradient-based optimization. We choose $M$ to be the smallest value for which the ELBO given by the gradient method converges to within a few nats of the exact marginal likelihood based on Fig. 5. We plot the optimization traces in Fig. 7 for several runs to account for random variation in the initializations.

In this constrained setting, we see different behaviours on the different data sets. One constant is that placing inducing points using K-means leads to sub-optimal performance compared to the best method. For the Energy data set, “greedy var (reinit)” suffers from convergence to local optima. This is caused by the low sparsity, and disappears if more inducing points are used (see Fig. 6). For the Naval data set, we see very slow convergence when using gradient-based optimization initialized with greedy variance selection. K-means underperforms and also suffers from local optima, with reinitialization reliably reaching the best ELBO. For elevators, reinitialization reaches the optimal ELBO fastest.

We note that in the reinitialization method the hyperparameter optimization step was terminated when L-BFGS had determined convergence according to the default Scipy settings. This leads to a characteristic “step” pattern in the optimization traces, where progress halts for many iterations towards the end of a hyperparameter optimization phase, followed by large gains after a reinitialization of the inducing inputs. By terminating the hyperparameter optimization earlier after signs of stagnation, the reinitialization method could be significantly sped up. In addition, we measure computational cost through the number of function evaluations. This does not take into account the additional cost of computing the gradients for the inducing inputs, which make up the bulk of parameters that are to be optimized. As the amount of computation needed to reinitialization the inducing points is comparable to the computation required in as single iteration of gradient descent, Fig. 7 likely understates the computational savings of the reinitialization method.

2.3 Recommendation for Inducing Input Selection

The main conclusion from our empirical results is that well-chosen inducing inputs that are reinitialized during hyperparameter optimization give highly accurate variational approximations to the results of exact GPs. While performing gradient-based optimization of the inducing inputs may lead to improved performance in settings that are constrained to be very sparse, in some instances it is not worth the additional effort. We provide a GPflow-based (Matthews et al., 2017) implementation of the initialization methods and experiments that builds on other open source software (Coelho, 2017; Virtanen et al., 2020), available at https://github.com/markvdw/RobustGP.

It is important to note that we only considered data sets where sparse approximations were practically possible. The “kin40k” UCI data set is a notable example where a squared exponential GP regression model with learned hyperparameters could not accurately be approximated, due to a lengthscale that continuously decreased with increasing $M$. Given the underfitting and significant hyperparameter bias (Bauer et al., 2016), one can question whether variational approximations are appropriate. In cases where the covariates are less heavily correlated under the prior, conjugate gradient approaches (Gibbs and Mackay, 1997; Davies, 2015; Gardner et al., 2018) may be better. We choose to not make a recommendation for how to choose inducing variables in cases where the variational approximation is poor.

Conclusions

We provide guarantees on the quality of variational sparse Gaussian process regression when many fewer inducing variables are used than data points. We also consider lower bounds on the number of inducing variables needed in order to ensure that the KL-divergence between the approximate posterior and the full posterior is not large. These bounds provide insight into the number of inducing points that should be used for a variety of tasks, as well as suggest the sorts of problems to which sparse variational inference is well-suited. We also include an empirical results comparing the efficacy of different methods for selecting inducing inputs, which is of practical importance to the Gaussian process community. We believe that there is a great deal of interesting future research to be done on the role of sparsity in variational Gaussian process inference; both in refining the bounds given in this work and in better understanding non-conjugate inference schemes, such as those developed in Hensman et al. (2015).

We would particularly like to thank Guillaume Gautier for pointing out an error in the exact k-DPP sampling algorithm cited in an earlier version of this work, and for guiding us through recent work on sampling k-DPPs that led to an amended proof. MvdW would additionally like to thank James Hensman for his guidance, and PROWLER.io for providing an excellent research environment while this work was developed.

A Proof of Bound on Mean and Variance of One-dimensional Marginal Distributions

Proof By the chain rule of KL-divergence, we have

For any $x^{\star}\in\mathcal{X}$, the KL-divergence on the right hand side of Eq. 30 is a KL-divergence between one-dimensional Gaussian distributions, and has the form,

Define $r=\sigma_{1}^{2}/\sigma_{2}^{2}$, so Eq. 31 becomes $\gamma\geq r-1-\log r$. For $\gamma<\frac{1}{5}$, we have $r-\log(r)<1.2$, so $r\in[.493,1.78]$. For $r$ in this range, we have, $\gamma\geq r-1-\log r\geq(r-1)^{2}/3$. Solving, for $r$, we obtain the bound,

We now turn to the proof of the bound relating $\mu_{1}$ and $\mu_{2}$. From Eq. 31 and because $r-1-\log r>0$ for $r>0$, $\frac{(\mu_{1}-\mu_{2})^{2}}{\sigma_{2}^{2}}\leq\gamma$. Rearranging, $|\mu_{1}-\mu_{2}|\leq\sigma_{2}\sqrt{\gamma}$. The final bound on the mean follows from Eq. 32, which implies that,

B Proofs of A-Posteriori Bounds

In this section, we restate and prove the upper bound on the marginal likelihood given in Titsias (2014). See 2

This result relies on several properties of symmetric positive semi-definite (SPSD) matrices, which we state in Proposition 35.

Let $\succ$ denote the partial order on SPSD matrices induced by $A\succ B\iff A-B$ is SPSD. Then if $A\succ B$ are $N\times N$ SPSD matrices,

If $A^{-1},B^{-1}$ exist, then $A^{-1}\prec B^{-1}$.

If $\lambda_{1}(A)\geq\dotsc\geq\lambda_{N}(A)$, $\lambda_{1}(B)\geq\dotsc\geq\lambda_{N}(B)$ denote the eigenvalues of $A$ and $B$ respectively, $\lambda_{i}(A)\geq\lambda_{i}(B)$ for all $1\leq i\leq N$.

We also use that $\textup{Q}_{\textup{ff}}\prec\textup{K}_{\textup{ff}}$, which follows from properties of Schur complements of PSD matrices (Gallier, 2010, Proposition 2.1).

Proof of Lemma 2 This proof follows that of Titsias (2014). Recall Eq. 6,

The first inequality uses $\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I}\succ\textup{Q}_{\textup{ff}}++\sigma^{2}\textup{I}$, which implies $\log\det(\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I})\geq\log\det(\textup{Q}_{\textup{ff}}+\sigma^{2}\textup{I})$. The second inequality uses that $\textup{K}_{\textup{ff}}\prec\textup{Q}_{\textup{ff}}+\|\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}}\|_{\textup{op}}\textup{I}$ and the second part of Proposition 35.

In problems where sparse GP regression is applied, computing the largest eigenvalue of $\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}}$ is computationally prohibitive. However, we can use the upper bound $\|\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}}\|_{\textup{op}}\leq\textup{tr}(\textup{K}_{\textup{ff}}-\textup{Q}_{\textup{ff}})$, yielding

The bound $\mathcal{U}_{2}$ can be computed in time $\mathcal{O}(NM^{2})$ with memory $\mathcal{O}(NM)$ in much the same way as the ELBO is computed, as it only depends on the low-rank matrix $\textup{Q}_{\textup{ff}}$ and the diagonal entries of $\textup{K}_{\textup{ff}}$.

C Proofs for Results Leading to Upper bounds on the KL-divergence

In this appendix, we restate and provide proofs for the results in Section 4. See 3

Proof We apply the matrix identity $(A+B)^{-1}=A^{-1}-A^{-1}B(A+B)^{-1}$ to the expression

with $A=\textup{Q}_{\textup{ff}}+\sigma^{2}\textup{I}$ and $B=\zeta\textup{I}$. This gives

The matrix $\textup{Q}_{\textup{ff}}^{2}+(\zeta+\sigma^{2})\textup{Q}_{\textup{ff}}$ is SPSD, as it is the product of SPSD matrices that commute. This implies that the eigenvalues of $\textup{Q}_{\textup{ff}}^{2}+(\zeta+\sigma^{2})\textup{Q}_{\textup{ff}}+\sigma^{2}(\zeta+\sigma^{2})\textup{I}$ are bounded below by $\sigma^{2}(\zeta+\sigma^{2})$. As the eigenvalues of the inverse of a SPSD matrix are the inverse of the eigenvalues of the original matrix, the largest eigenvalue of $(\textup{Q}_{\textup{ff}}^{2}+(\zeta+\sigma^{2})\textup{Q}_{\textup{ff}}+\sigma^{2}(\zeta+\sigma^{2})\textup{I}))^{-1}$ is bounded above by $(\sigma^{2}(\zeta+\sigma^{2}))^{-1}$. Therefore,

This proves the second inequality in Lemma 3. The same argument using $\mathcal{U}_{2}$ in place of $\mathcal{U}_{1}$ yields,

See 4 We have already proven the lower bound in the main body. In order to prove the upper bound in Lemma 4, we use a Hölder-type inequality, Tao (2012, Exercise 1.3.26).

In particular, if $A$ and $B$ are SPSD (so that the singular values agree with the eigenvalues), taking $p=1,$ $q=\infty$,

The inequality uses that $\textup{Q}_{\textup{ff}}+\sigma^{2}\textup{I}\prec\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I}$, so $\det(\textup{Q}_{\textup{ff}}+\sigma^{2}\textup{I})\leq\det(\textup{K}_{\textup{ff}}+\sigma^{2}\textup{I})$ by Proposition 35. We can now apply Proposition 36 with $p=1,q=\infty$ to Section C giving,

Using this bound in Section C and combining with Eq. 15 completes the proof of the upper bound.

D Derivations of bounds for specific kernels and covariate distributions

In this appendix, we restate and provide proofs for the results in Section 5.

Proof of Corollary 19 Using Eq. 24 and applying the geometric series formula,

We can use this equation in Theorem 13 (a similar result could be obtained using Theorem 14) yielding,

Choose $\epsilon=\frac{\gamma\sigma^{2}}{2Nv(1+RN/\sigma^{2})}=\Theta(\gamma/N^{2})$. By Lemma 9, an $M$-DPP can be sampled to this level of accuracy using not more than $\mathcal{O}(NM(\log\frac{N^{2}}{\gamma\delta}))$ iterations of MCMC, making the computational cost of selecting inducing inputs $\mathcal{O}(NM^{3}(\log\frac{N^{2}}{\gamma\delta}))$. We may assume that $M<N$, otherwise by choosing $Z=X$ the KL-divergence is zero and nothing more needs to be shown. Then,

Take $M=\log_{B}\sqrt{\frac{A}{2a}}\frac{\gamma\delta\sigma^{2}(1-B)}{N^{2}(1+RN/\sigma^{2})}=\mathcal{O}(\log(N^{3}/\gamma\delta))$, then

In the case of ridge leverage score initializations, from Theorem 17 we have with probability $1-5\delta$,

D.2 Bounds for Multivariate Gaussian distributions and Squared Exponential Kernel

Proof of Proposition 21 The proof of this proposition is nearly identical to an argument in Seeger et al. (2008). Consider the upper bound,

where in the second to last line we make the substitution $t=\alpha s^{1/D}$ and in the final line we recognized the integral as an incomplete $\Gamma$-function.

From Gradshteyn and Ryzhik (2014, 8.352) for integer $D$ and $r>0$,

For fixed $D$ and $r$ large (which is satisfied by the condition $M\geq\frac{1}{\alpha}D^{D}+D-1$), we have that the final term in the above sum is the largest, so that

Proof Corollary 22 is a consequence of Theorem 13 and Proposition 21. In the case of the $M$-DPP, we take $\epsilon=\frac{\gamma\sigma^{2}}{2Nv(1+RN/\sigma^{2})}=\Theta(\gamma/N^{2})$ as in the proof of Corollary 19. It then remains to choose $M$ so that

From Proposition 21, there exists an $M=\mathcal{O}((\log\frac{N^{3}}{\gamma})^{D})$ that satisfies this criteria. In the case of ridge leverage scores, it is sufficient to choose S=\mathcal{O}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(\log\frac{N^{2}}{\delta\gamma}}\right)^{D}}\right), which means that with probability at least $1-5\delta$, \mathbf{M}=\mathcal{O}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left(\log\frac{N^{2}}{\delta\gamma}}\right)^{D}(\log\log\frac{N^{2}}{\delta\gamma}+\log(1/\delta))}\right).

D.3 Conditions for Widom’s Theorem

Widom’s Theorem (Widom, 1963), states that for stationary kernels on compact subsets of Euclidean space, the eigenvalues of the operator $\mathcal{K}$ are closely linked to the decay of the spectral density of the kernel function. The theorem applies to any compactly supported covariate distribution with Lebesgue density and stationary kernel with spectral density satisfying the following three conditions:

E Lower bounds on the number of features

In this appendix, we restate and prove the results stated in Section 6.

In the final line, we use that $\mathbf{Q}_{\bf ff}$ is at most rank $M$, so that $\bm{\psi}_{m}=0$ for all $m>M$. It follows from the inequality $\log(1+a)\leq a$ for $a\geq 0$ that each term in the first sum is non-negative. Hence,

E.2 Lower Bound on Eigenvalues of Multivariate Gaussian Inputs and Squared Exponential Kernel

Proof Recall from Section 5.1.1 that the eigenvalues of this operator are of the form,

where the number of times each eigenvalue is repeated is equal to the number of ways to write $s$ as a sum of $D$ non-negative integers, where the order of the summands matters. This is equal to $\binom{s+D-1}{D-1}$. The number of eigenvalues greater than \mathopen{}\mathclose{{}\left(2a/A}\right)^{D/2}B^{s} is therefore,

The equality follows from observing that the right hand side is equal to the number of way to write $s$ as a sum of $D+1$ non-negative integers. For each of these representations, the first $D$ integers sum to some $t\leq s$, and once these are fixed there is a unique choice for the final integer. This is equivalent to the left hand side. We therefore conclude \lambda_{\binom{s+D}{D}}=\mathopen{}\mathclose{{}\left(\frac{2a}{A}}\right)^{D/2}B^{s}. Define

Proof By Lemma 28 and Proposition 21, for $\delta\in(0,1)$, with probability $1-\delta$,

with $\alpha=-\log B$ for any $1\leq r\leq N$. Using Proposition 30 we have,

For $\gamma\in(0,1/2)$, choose $r=\lceil(\frac{1}{\alpha D}\log N^{\gamma})^{D}\rceil$, then noting that for this choice of $r$, the third term in the sum is smaller than the second term,

Applying Proposition 30, with $M+1=\lfloor(\frac{1}{D}\log_{B}N^{-\zeta/D})^{D}\rfloor$ with $\zeta\in(0,\gamma)$. We have

F Effect of Jitter on Bounds

The final inequality follows from $\textup{K}_{\textup{uu}}+\epsilon\textup{I}\prec\textup{K}_{\textup{uu}}+\epsilon^{\prime}\textup{I}$ and Proposition 35. Therefore, $\textup{Q}_{\textup{ff}}(\epsilon^{\prime})\prec\textup{Q}_{\textup{ff}}(\epsilon)$. From Proposition 35, we have

Let $A,B$ arbitrary $N\times N$ SPSD matrices with $A\succ B\succ\sigma^{2}\textup{I}$. Denote the eigenvalues of $A$ and $B$ respectively as $\lambda_{1}(A)\geq\dots\lambda_{N}(A)$ and $\lambda_{1}(B)\geq\dots\lambda_{N}(B)$. Then,

The first inequality follows applying $\log(1+a)\leq a$ to each term in the sum. The second inequality used that $\lambda_{i}(B)\geq\sigma^{2}$ since $B\succ\sigma^{2}\textup{I}$. Then,

where the final inequality follows from Eq. 40 with $A=\textup{Q}_{\textup{ff}}(\epsilon)+\sigma^{2}\textup{I}$ and $B=\textup{Q}_{\textup{ff}}(\epsilon^{\prime})+\sigma^{2}\textup{I}$. Combining Eq. 39 with Eq. 41 proves the monotonicity of the lower bound in $\epsilon$. The upper bound follows from Proposition 35 noting that in the quadratic form

G An Alternative Ridge Leverage Sampling Initialization

Many implementations of leverage score sampling allow for adaptively selecting the number of inducing points to achieve a desired level of accuracy. We briefly discuss the application of Algorithm 2 in Musco and Musco (2017) to the problem of sparse variational inference in Gaussian processes.

The number of points sampled by ridge leverage score methods to achieve a desired level of accuracy is closely related to the effective dimension of the kernel matrix, which can be thought of as measure of the complexity of the non-parameteric regression model. The effective dimension is defined as the sum of the ridge leverage scores,

and depends on the choice of kernel, the distribution of the covariates and the regularization parameter.

In order to compare such an adaptive method with the bounds discussed in Section 4, we need to consider the typical size of the effective dimension, assuming a fixed kernel and a random set of covariates with identical marginal distributions (or marginal distributions satisfying the conditions in Lemma 11).

For any fixed set of covariates, we can split the sum in Eq. 42 into two parts, yielding

where $S$ is an arbitrary positive integer. Upper bounds on the effective dimension can be obtained by choosing $S$ so that the two terms on the right hand side of Eq. 43 are of the same order of magnitude.

G.2 Adaptively Selecting the Number of Inducing Points with Leverage Scores

We consider the application of Musco and Musco (2017, Algorithm 2) to the problem of selecting inducing inputs for sparse variational inference in GP models. This algorithm comes with the following bounds on the quality of the resulting Nyström approximation.

Fix $\delta\in(0,\frac{1}{32})$. There exists an algorithm with run time $\mathcal{O}(NM^{2})$ and memory complexity $\mathcal{O}(NM)$ that with probability $1-3\delta$ returns $M<384d_{\textup{eff}}^{\omega}\log(d_{\textup{eff}}^{\omega}/\delta)$ columns of $\textup{K}_{\textup{ff}}$ such that the resulting Nyström approximation, $\textup{Q}_{\textup{ff}}$, satisfies

where $d_{\textup{eff}}^{\omega}$ denotes the effective dimension of the Gaussian process regressor with $\sigma^{2}=\omega$.Note that $\textup{K}_{\textup{ff}}$ and $\textup{Q}_{\textup{ff}}$ are both independent of the noise parameter, so there is no requirement that the ‘noise parameter’ used for initializing inducing points matches the noise parameter used in performing regression.

We can now consider the implications of this bound on sparse variational GP regression using Lemmas 3 and 4.

G.3 Ridge Leverage Scores and Sparse Variational Inference

Both of these bounds are small if $\|\mathbf{K}_{\bf ff}-\mathbf{Q}_{\bf ff}\|_{\textup{op}}\ll 1/N$.

For simplicity, we consider the case when $\mathbf{y}$ is assumed to have a conditional distribution that agrees with the GP prior. Fix $\delta\in(0,1/32)$ and $\gamma>0$. Applying Markov’s inequality to Eq. 44, with probability at least $1-\delta$,

We can apply the algorithm referred to in Lemma 37 with $\omega=\sigma^{2}\delta\gamma/N$, so that with probability at least $1-3\delta$ a set of inducing inputs is chosen such that,

We can then apply a union bound to conclude with probability at least $1-4\delta$,

By Corollary 12 and Markov’s inequality, with probability at least $1-\delta$,

for any $1\leq S\leq N$. On the event where this holds and recalling we chose the parameter $\omega=\sigma^{2}\delta\gamma/N$, Eq. 43 implies that,

We can again apply the union bound to lower bound the probability that both the effective dimension is less that the bound in Eq. 47 and that Eq. 46 holds. This yields the following probabilistic bounds on the quality of sparse VI in GP regression with inducing points placed according to approximate ridge leverage scores.

Fix $\delta\in(0,\frac{1}{32}),$ $\gamma>0$. Under the same assumptions on the covariate distribution and the distribution of $\mathbf{y}$ as in Theorem 14 if inducing points are placed according to Musco and Musco (2017, Algorithm 2) with $\omega=\sigma^{2}\delta\gamma/N$, then with probability $1-5\delta$, $\mathbf{M}<384d\log(d/\delta)$ and

A similar argument in the case when we do not assume $\mathbf{y}$ is distributed according to the prior model leads to the following result:

Fix $\delta\in(0,\frac{1}{32}),$ $\gamma>0$. Under the same assumptions on the covariate distribution and the distribution of $\mathbf{y}$ as in Theorem 13 if inducing points are placed according to Musco and Musco (2017, Algorithm 2) with $\omega=\frac{2\sigma^{2}\delta\gamma}{N(1+R/\sigma^{2})}$ then with probability $1-5\delta$, $\mathbf{M}<384d^{\prime}\log(d^{\prime}/\delta)$ and

Note that while the resulting bounds on $\mathbf{M}$ depend on the kernel and covariate distribution, the quality of the resulting approximation in both Theorems 38 and 39 does not.

The bounds implied by these results for various kernels are given in Table 3. Note that the asymptotic rates implied by both Theorems 38 and 39 are the same. This is because, unlike in the case of the $M$-DPP initialization in which the trace is bounded and this is used as an upper bound on the operator norm, the operator norm is bounded directly.

References