Rates of Convergence for Sparse Variational Gaussian Process Regression

Introduction

Gaussian processes (GPs) [Rasmussen & Williams, 2006] are distributions over functions that are convenient priors in Bayesian models. They can be seen as infinitely wide neural networks [Neal, 1996], and are particularly popular in regression models, as they produce good uncertainty estimates, and have closed-form expressions for the posterior and marginal likelihood. The most well known drawback of GP regression is the computational cost of the exact calculation of these quantities, which scales as \mathcal{O}\mathopen{}\mathclose{{}\left(N^{3}}\right) in time and \mathcal{O}\mathopen{}\mathclose{{}\left(N^{2}}\right) in memory where $N$ is the number of training examples. Low-rank approximations [Quiñonero Candela & Rasmussen, 2005] choose $M$ inducing variables which summarise the entire posterior, reducing the cost to \mathcal{O}\mathopen{}\mathclose{{}\left(NM^{2}+M^{3}}\right) time and \mathcal{O}\mathopen{}\mathclose{{}\left(NM+M^{2}}\right) memory.

While the computational cost of adding inducing variables is well understood, results on how many are needed to achieve a good approximation are lacking. As the dataset size increases, we cannot expect to keep the capacity of the approximation constant without the quality deteriorating. Taking into account the rate at which $M$ must increase with $N$ to achieve a particular approximation accuracy, as well as the cost of initializing or optimizing the inducing points, determines a more realistic sense of the costs of scaling Gaussian processes.

Approximate GPs are often trained using variational inference [Titsias, 2009], which minimizes the KL divergence from an approximate posterior to the full posterior process [Matthews et al., 2016]. We use this KL divergence as our metric for the approximate posterior’s quality. We show that under intuitive assumptions the number of inducing variables only needs to grow at a sublinear rate for the KL between the approximation and the posterior to go to zero. This shows that very sparse approximations can be used for large datasets, without introducing much bias into hyperparameter selection through evidence lower bound (ELBO) maximisation, and with approximate posteriors that are accurate in terms of their prediction and uncertainty.

The core idea of our proof is to use upper bounds on the KL divergence that depend on the quality of a Nyström approximation to the data covariance matrix. Using existing results, we show this error can be understood in terms of the spectrum of an infinite-dimensional integral operator. In the case of stationary kernels, our main result proves that priors with smoother sample functions, and datasets with more concentrated inputs admit sparser approximations.

We assume that training inputs are drawn i.i.d. from a fixed distribution, and prove bounds of the form

with probability at least $1-\delta$, where $\hat{P}$ is the posterior Gaussian process, $Q$ is a variational approximation, and $\mathbf{y}$ are the training targets. The function $g(M,N)$ depends on both the kernel and input distribution, and grows linearly in $N$ and generally decays rapidly in $M$. The quality of the initialization determines $\epsilon$, which can be made arbitrarily small (e.g. an inverse power of $N$) at some additional computational cost. Theorems 1 and 2 give results of this form for a collection of inducing variables defined using spectral information, theorems 3 and 4 hold for inducing points.

Background and notation

where ${\bf K}_{n}=\mathbf{K}_{\bf ff}+\sigma_{n}^{2}\mathbf{I}$, and \mathopen{}\mathclose{{}\left[\mathbf{K}_{\bf ff}}\right]_{i,j}=k(\mathbf{x}_{i},\mathbf{x}_{j}).

2 Sparse variational Gaussian process regression

While all quantities of interest have analytic expressions, their computation is infeasible for large datasets due to the \mathcal{O}\mathopen{}\mathclose{{}\left(N^{3}}\right) time complexity of the determinant and inverse. Many approaches have been proposed that rely on a low-rank approximation to $\mathbf{K}_{\bf ff}$ [Quiñonero Candela & Rasmussen, 2005; Rahimi & Recht, 2008], which allow the determinant and inverse to be computed in \mathcal{O}\mathopen{}\mathclose{{}\left(NM^{2}}\right), where $M$ is the rank of the approximating matrix.

We consider the variational framework developed by Titsias , which minimizes the KL divergence from the posterior process to an approximate GP

Titsias suggests jointly maximizing the ELBO (eq. 3) w.r.t. the variational and hyperparameters. This comes at the cost of introducing bias in hyperparameter estimation [Turner & Sahani, 2011], notably the overestimation of the $\sigma_{n}^{2}$ [Bauer et al., 2016]. Adding inducing points reduces the KL gap [Titsias, 2009], and the bias is practically eliminated when enough inducing variables are used.

3 Interdomain inducing features

Lázaro-Gredilla & Figueiras-Vidal showed that one can specify the distribution $q(\mathbf{u}),$ on integral transformations of $f(\cdot).$ Using these interdomain inducing variables can lead to sparser representations, or computational benefits [Hensman et al., 2018]. Interdomain inducing variables are defined by

When $g(\mathbf{x};\mathbf{z}_{m})=\delta(\mathbf{x}-\mathbf{z}_{m})$ the $u_{m}$ are inducing points. Interdomain features require replacing $\mathbf{k}_{\mathbf{u}\cdot}$ and $\mathbf{K}_{\bf uu}$ in eq. 2 with integral transforms of the kernel. In later sections, we investigate particular interdomain transformations with interesting convergence properties.

4 Upper bounds on the marginal likelihood

Combined with eq. 3, an upper bound on eq. 1 can show when the KL divergence is small, which indicates inference has been successful and hyperparameter estimates are likely to have little bias. Titsias introduced an upper bound that can be computed in \mathcal{O}\mathopen{}\mathclose{{}\left(NM^{2}}\right):

This gives a data-dependent upper bound, that can be computed after seeing the data, and for given inducing inputs.

5 Spectral properties of the covariance matrix

While for small datasets spectral properties of the covariance matrix can be analyzed numerically, we need a different approach for understanding these properties for a typical large dataset. The covariance operator, $\mathcal{K},$ captures the limiting properties of $\mathbf{K}_{\bf ff}$ for large $N$. It is defined by

where $p(\mathbf{x})$ is a probability density from which the inputs are assumed to be drawn. We assume that $\mathcal{K}$ is compact, which is the case if $k(\mathbf{x},\mathbf{x}^{\prime})$ is bounded. Under this assumption, the spectral theorem tells us that $\mathcal{K}$ has a discrete spectrum. The (finite) sequence of eigenvalues of $\frac{1}{N}\mathbf{K}_{\bf ff}$ converges to the (infinite) sequence of eigenvalues of $\mathcal{K}$ [Koltchinskii & Giné, 2000]. Mercer’s Theorem [Mercer, 1909] tells us that for continuous kernel functions,

where the \mathopen{}\mathclose{{}\left(\lambda_{m},\phi_{m}}\right)_{i=1}^{\infty} are eigenvalue-eigenfunction pairs of the operator $\mathcal{K},$ with the eigenfunctions orthonormal in $L^{2}(\mathcal{X})_{p}.$ Additionally, $\sum_{m=1}^{\infty}\lambda_{m}<\infty$.

6 Selecting the number of inducing variables

Ideally, the number of inducing variables should be selected to make the $KL(Q||\hat{P})$ small. Currently, the most common advice is to increase the number of inducing variables $M$ until the lower bound (eq. 3) no longer improves. This is a necessary, but not a sufficient condition for the ELBO to be tight and the KL to be small. Taking the upper bound (eq. 4) into consideration, we can guarantee a good approximation when difference between the upper and lower bounds converges to zero, as this upper bounds the KL.

Both these procedures rely on bounds computed for a given dataset, and a specific setting of variational parameters. While practically useful, they do not tell us how many inducing variables we should expect to use before observing any data. In this work, we focus on a priori bounds, and asymptotic behavior as $N\to\infty$ and $M$ grows as a function of $N$. These bounds guarantee how the variational method scales computationally for any dataset satisfying intuitive conditions. This is particularly important for continual learning scenarios, where we incrementally observe more data. With our a priori results we can guarantee that the growth in required computation will not exceed a certain rate.

Bounds on the KL divergence for eigenfunction inducing features

In this section, we prove a priori bounds on the KL divergence using inducing features that rely on spectral information about the covariance matrix or the associated operator. The results in this section form the basis for bounds on the KL divergence for inducing points (section 4).

We first consider a posteriori bounds on the KL divergence that hold for any $\mathbf{y},$ derived by looking at the difference between $\mathcal{L}_{\text{upper}}$ and $\mathcal{L}_{\text{lower}}.$ We will use these bounds in later sections to analyze asymptotic convergence properties.

Let $\widetilde{\mathbf{K}}_{\bf ff}=\mathbf{K}_{\bf ff}-\mathbf{Q}_{\bf ff},$ $t=\operatorname{Tr}(\widetilde{\mathbf{K}}_{\bf ff})$ and $\widetilde{\lambda}_{max}$ denote the largest eigenvalue of $\widetilde{\mathbf{K}}_{\bf ff}$. Then,

The proof bounds the difference between a refinement of $\mathcal{L}_{\text{upper}}$ also proven by Titsias and $\mathcal{L}_{\text{lower}}$ through an algebraic manipulation and is given in appendix A. The second inequality is a consequence of $t\geq\widetilde{\lambda}_{max}.$ We typically expect $\|\mathbf{y}\|_{2}^{2}=\mathcal{O}(N)$, which is the case when the variance of the observed $y$s is bounded, so if $t\ll 1/N$ the KL divergence will be small.

2 A priori bounds: averaging over 𝐲𝐲\mathbf{y}

Lemma 1 is typically overly pessimistic, as it assumes $\mathbf{y}$ can be parallel to the largest eigenvector of $\widetilde{\mathbf{K}}_{\bf ff}.$ In this section, we consider a bound that holds a priori over the training outputs, when they are drawn from the model. This allows us to bound the KL divergence for a ‘typical’ dataset.

For any set of $\{\mathbf{x}_{i}\}_{i=1}^{N}$, if the outputs $\{y_{i}\}_{i=1}^{N}$ are generated according to our generative model, then

The lower bound tells us that even if the training data is contained in an interval of fixed length, we need to use more inducing points for problems with large $N$ if we want to ensure the sparse approximation has converged. This is shown in Figure 1 for data uniformly sampled on the interval $$ with 15 inducing points.

The second term on the right is a KL divergence between centered Gaussian distributions. The lower bound follows from Jensen’s inequality. The proof of the upper bound (appendix B) bounds this KL divergence above by $t/(2\sigma_{n}^{2}).$ ∎

3 Minimizing the upper bound: an idealized case

We now consider the set of $M$ interdomain inducing features that minimize the upper bounds in Lemmas 1 and 2. Taking into account the lower bound in Lemma 2, they must be within a factor of two of the optimal features defined without reference to training outputs under the assumption of Lemma 2. Consider

where $w_{i}^{(m)}$ is the $i^{th}$ entry in the $m^{th}$ eigenvector of $\mathbf{K}_{\bf ff}.$ That is, $u_{m}$ is a linear combination of inducing points placed at each data point, with weights coming from the entries of the $m^{th}$ eigenvector of $\mathbf{K}_{\bf ff}.$ We show in appendix C that

Inference with these features can be seen as the variational equivalent of the optimal parametric projection of the model derived by Ferrari-Trecate et al. .

4 Eigenfunction inducing features

We now modify the construction given in section 3.3 to no longer depend on $\mathbf{K}_{\bf ff}$ explicitly (which depends on the specific training inputs) and instead depend on assumptions about the training data. This construction is the a priori counterpart of the eigenvector inducing features, as it is defined prior to observing a specific set of training inputs.

Consider the limit as we have observed a large amount of data, so that $\frac{1}{N}\mathbf{K}_{\bf ff}\to\mathcal{K}.$ This leads us to replace the eigenvalues, $\{\lambda_{m}(\mathbf{K}_{\bf ff})\}_{m=1}^{M},$ with the operator eigenvalues, $\{\lambda_{m}\}_{m=1}^{M},$ and the eigenvectors, $\{\mathbf{w^{(m)}}\}_{m=1}^{M},$ with the eigenfunctions, $\{\phi_{m}\}_{m=1}^{M},$ yielding

Note that $p(\mathbf{x})$ influences $u_{m}$. In appendix C, we show

These features can be seen as the variational equivalent of methods utilizing truncated priors proposed in Zhu et al. , which are the optimal linear $M$ dimensional parametric GP approximation defined a priori, in terms of minimizing expected mean square error.

In the case of the SE kernel and Gaussian inputs, closed form expressions for eigenfunctions and values are known [Zhu et al., 1997]. For Matérn kernels with inputs uniform on $[a,b]$, expressions for the eigenfunctions and eigenvalues needed in order to compute $\mathbf{K}_{\bf uf}$ and $\mathbf{K}_{\bf uu}$ can be found in Youla . However, the formulas involve solving systems of transcendental equations limiting the practical applicability of these features for Matérn kernels.

5 A priori bounds on the KL divergence for eigenfunction features

Having developed the necessary preliminary results, we now prove the first a priori bounds on the KL divergence. We start with eigenfunction features, which can be implemented practically in certain instances discussed above.

Suppose $N$ training inputs are drawn i.i.d. according to input density $p(\mathbf{x}).$ For inference with $M$ eigenfunction inducing variables defined with respect to the prior kernel and $p(\mathbf{x}),$ with probability at least $1-\delta,$

where we have defined $C=N\sum_{m=M+1}^{\infty}\lambda_{m},$ and the $\lambda_{m}$ are the eigenvalues of the integral operator $\mathcal{K}$ associated to the prior kernel and $p(\mathbf{x}).$

With the assumptions and notation of Theorem 1 if $\mathbf{y}$ is distributed according to a sample from the prior generative model, with probability at least $1-\delta,$

We first prove a bound on $t$ that holds in expectation over input data matrices of size $N$ with entries drawn i.i.d. from $p(\mathbf{x}).$ A direct computation of $\mathbf{Q}_{\bf ff}$ shows that \mathopen{}\mathclose{{}\left[\mathbf{Q}_{\bf ff}}\right]_{i,j}=\sum_{m=1}^{M}\lambda_{m}\phi_{m}(\mathbf{x}_{i})\phi_{m}(\mathbf{x}_{j}). Using the Mercer expansion of the kernel matrix and subtracting,

The second equality follows from the eigenfunctions having norm 1. Applying Markov’s inequality and Lemmas 1 and 2 yields Theorems 1 and 2 respectively. ∎

6 Squared exponential kernel and Gaussian inputs

Using this bound with Theorems 1 and 2, we see that by choosing $M=\mathcal{O}(\log N),$ under the assumptions of either theorem, we can obtain a bound on the KL divergence that tends to as $N$ tends to infinity.

7 Matérn kernels and uniform measure

For the Matérn $k+1/2$ kernel in one dimension, $\lambda_{m}\asymp m^{-2k-2}$ [Ritter et al., 1995; Seeger et al., 2008], so $\sum_{m=M+1}^{\infty}\lambda_{m}=\mathcal{O}(M^{-2k-1}).$ In order for the bound in Theorem 2 to converge to $0,$ we need $\lim\limits_{N\to\infty}\frac{N}{M^{2k+1}}\to 0.$ This holds if $M=N^{\alpha}$ for $\alpha>\frac{1}{2k+1}.$ For $k>0,$ this bound indicates the number of inducing features can grow sublinearly with the amount of data.

Bounds for inducing points

We have shown that using spectral knowledge of either $\mathbf{K}_{\bf ff}$ or $\mathcal{K}$ we obtain bounds on the KL divergence indicating that the number of inducing features can be much smaller than the number of data points. While mathematically convenient, the practical applicability of the interdomain features used is limited by computational considerations in the case of the eigenvector features and by the lack of analytic expressions for $\mathbf{K}_{\bf uf}$ in most cases for the eigenfunction features, as well not knowing the true input density, $p(\mathbf{x})$.

In contrast, inducing points can be efficiently applied to any kernel. In this section, we show that with a good initialization based on the empirical input data distribution, inducing points lead to bounds that are only slightly weaker than the interdomain approaches suggested so far.

[Belabbas & Wolfe, 2009] Given a symmetric positive semidefinite matrix, $\mathbf{K}_{\bf ff},$ if $M$ columns are selected to form a Nyström approximation such that the probability of selecting a subset of columns, $Z,$ is proportional to the determinant of the principal submatrix formed by these columns and the matching rows, then,

This means that on average well-initialized inducing points lead to bounds within a factor of $M+1$ of eigenvector inducing features.

The selection scheme described introduces negative correlations between inducing points locations, leading the $\mathbf{z}_{i}$ to be well-dispersed amongst the training data, as shown in fig. 2. The strength of these negative correlations is determined by the particular kernel.

The proposed initialization scheme is equivalent to sampling $Z$ according to a discrete k-Determinantal Point Process (k-DPP), defined over $\mathbf{K}_{\bf ff}$. Belabbas & Wolfe suggested that sampling from this distribution, which has support over $\binom{N}{M}$ subsets of columns, may be computationally infeasible. Kulesza & Taskar provided an exact algorithm for sampling from k-DPPs given an eigendecomposition of the kernel matrix. In our setting, we require our initialisation scheme to have similar computational cost to computing the sparse GP bounds, which prohibits us from performing eigendecomposition. Instead, we rely on cheaper “$\epsilon$ close” sampling methods. We therefore provide the following corollary of lemma 3, proven in appendix D.

Suppose the inducing points, $Z,$ are sampled from an $\epsilon$ k-DPP, $\nu,$ i.e a distribution over subsets of $\mathbf{X}$ of size $M$ satisfying, $d(\mu,\nu)_{TV}\leq\epsilon$ where $d(\cdot,\cdot)_{TV}$ denotes total variation distance and $\mu$ is a k-DPP on $\mathbf{K}_{\bf ff}$. Suppose the $k(\mathbf{x},\mathbf{x})<v$ for all $\mathbf{x}\in\mathcal{X}.$ Then

We show analogues of theorems 1 and 2 for inducing points.

Suppose $N$ training inputs are drawn i.i.d according to input density $p(\mathbf{x}),$ and $k(\mathbf{x},\mathbf{x})<v$ for all $\mathbf{x}\in\mathcal{X}.$ Sample $M$ inducing points from the training data with the probability assigned to any set of size $M$ equal to the probability assigned to the corresponding subset by an $\epsilon$ k-DPP with $k=M$. With probability at least $1-\delta,$

where $C=N\sum_{m=M+1}^{\infty}\lambda_{m},$ and $\lambda_{m}$ are the eigenvalues of the integral operator $\mathcal{K}$ associated to kernel, $k,$ and $p(\mathbf{x}).$

With the assumptions and notation of theorem 3 and if $\mathbf{y}$ is distributed according to a sample from the prior generative model, with probability at least $1-\delta,$

We prove theorem 4. Theorem 3 follows the same argument, replacing the expectation over $\mathbf{y}$ with the bound given by lemma 1.

The first two inequalities use lemma 2 and corollary 1. The third follows from noting that the sum inside the expectation is the error in trace norm of the optimal rank $M$ approximation to the covariance matrix for any given $\mathbf{X}$, and is bounded above by the error from the rank $M$ approximation due to eigenfunction features. We showed that this error is in expectation equal to $N\sum_{m=M+1}^{\infty}\lambda_{m}$ so this must be an upper bound on the expectation in the second to last line. Shawe-Taylor et al. [2005, Proposition 4] gives a different proof of the final inequality.

We apply Markov’s inequality, yielding for any $\delta\in(0,1)$ with probability at least $1-\delta,$

Figure 3 compares the actual KL divergence, the a posteriori bound derived by $\mathcal{L}_{\text{upper}}-\mathcal{L}_{\text{lower}},$ and the bounds proven in theorems 3 and 4 on a dataset with normally distributed training inputs and $\mathbf{y}$ drawn from the generative model.

Consequences of theorem 3 and theorem 4

We now investigate implications of our main results for sparse GP regression. Our first two corollaries consider Gaussian inputs and the squared exponential (SE) kernel, and show that in $D$ dimensions, choosing $M=\mathcal{O}(\log^{D}(N))$ is sufficient in order for the KL divergence to converge with high probability. We then briefly summarize convergence rates for other stationary kernels. Finally we point out consequences of our definition of convergence for the quality of the pointwise posterior mean and uncertainty.

Using the explicit formula for the eigenvalues given in section 3.6, we arrive at the following corollary:

Suppose that $\|\mathbf{y}\|_{2}^{2}\leq RN.$ Fix $\gamma>0,$ and take $\epsilon=\frac{\delta\sigma_{n}^{2}}{vN^{\gamma+2}}.$ Assume the input data is normally distributed and regression in performed with a SE kernel. Under the assumptions of theorem 3, with probability $1-\delta,$

when inference is performed with $M=\frac{(3+\gamma)\log(N)+\log D}{\log(B^{-1})},$ where $D=\frac{v\sqrt{2a}}{2\sqrt{A}\sigma_{n}^{2}\delta(1-B)}.$

The proof is given in appendix E. If the lengthscale is much shorter than the standard deviation of the data then $B$ will be near 1, implying that $M$ will need to be large in order for the bound to converge.

The assumption $\|\mathbf{y}\|_{2}^{2}\leq RN$ for some $R$ is very weak. For example, if $\mathbf{y}$ is a realization of an integrable function with constant noise,

The first sum is asymptotically $N\int f(\mathbf{x})p(\mathbf{x})dx,$ and the second is asymptotically $N\sigma_{n}^{2}.$

The consequence of corollary 2 is shown in fig. 4, in which we gradually increase $N,$ choosing $M=C\log(N)+C_{0},$ and see the KL divergence converges as an inverse power of $N.$ The training outputs are generated from a sample from the prior generative model. Note that as theorem 4 assumes $\mathbf{y}$ is sampled from the prior and is not derived using the upper bound (eq. 4); it may be tighter than the a posteriori bound in cases when this upper bound is not tight.

For the SE kernel and Gaussian inputs, the rate that we prove $M$ must increase for inducing points and eigenfunction features differs by a constant factor. For the Matérn $k+1/2$ kernel in one dimension, we need to choose $M=N^{\alpha}$ with $\alpha>1/(2k)$ instead of $\alpha>1/(2k+1).$ This difference is particularly stark in the case of the Matérn $3/2$ kernel, for which our bounds tell us that inference with inducing points requires $\alpha>1/2$ as opposed to $\alpha>1/3$ for the eigenfunction features. Whether this is an artifact of the proof, the initialization scheme, or an inherent limitation for inducing points is an interesting area for future work.

2 Multidimensional data, effect of input density and other kernels

The proof uses ideas from Seeger et al. and is given in appendix E. While for the SE kernel and Gaussian input density $M$ can grow polylogarithmically in $N,$ and the KL divergence still converges, this is not the case for regression with other kernels or input distribution.

Closed form expressions for the eigenvalues of operators associated to many kernels and input distributions are not known. For stationary kernels and compactly supported input distributions, the asymptotic rate of decay of the eigenvalues of $\mathcal{K}$ is well-understood [Widom, 1963, 1964; Ritter et al., 1995]. The intuitive summary of these results is that smooth kernels, with concentrated input distributions have rapidly decaying eigenvalues. In contrast, kernels such as the Matérn-1/2 that define processes that are not smooth have slowly decaying eigenvalues. For Lebesgue measure on $[a,b]$ the Sacks-Ylivasker conditions of order $r$ (appendix F), which can be roughly thought of as meaning that realizations of the process are $r$ times differentiable with probability 1 [Ritter et al., 1995], implies an eigendecay of $\lambda_{m}\asymp m^{-2r-2}.$ Table 1 summarizes the spectral decay of several stationary kernels, as well as the implications for the number of inducing points needed for inference to provably converge with our bounds.

3 Computational complexity

We now have all the components necessary for analyzing the overall computational complexity of finding an arbitrarily good GP approximation. To understand the full computational complexity, we must consider the cost of initializing the inducing points using an exact or approximate k-DPP, as well as the $\mathcal{O}(NM^{2})$ time complexity of variational inference. Recent work of Dereziǹski et al. indicate that an exact algorithm for sampling a k-DPP can be implemented in $\mathcal{O}(N\log(N)\text{poly}(M)).$ We base our method on Anari et al. , who show that an $\epsilon$ k-DPP can be sampled via MCMC methods in $\mathcal{O}(NM^{4}\log(N)\!+\!NM^{3}\log(\frac{1}{\epsilon}))$ time with memory $\mathcal{O}(N\!+\!M^{2})$ (see appendix D).Open source implementations of approximate k-DPPs are available (e.g. [Gautier et al., 2018]). We can choose $\epsilon$ to be any inverse power of $N$ which only adds a constant factor to the complexity. For the SE kernel, taking $M\!=\!\mathcal{O}(\log^{D}\!N)$ inducing points leads to a complexity of $\mathcal{O}(N\log^{4D+1}N),$ a large computational saving compared to the $\mathcal{O}(N^{3})$ cost of exact inference. For the Matérn $k\!+\!\frac{1}{2}$ kernel in one-dimension and the average case analysis of theorem 4 we need to take $M\!=\!\mathcal{O}(N^{1/(2k)+\epsilon^{\prime}})$ implying a computational complexity of $\mathcal{O}(N^{1+2/k+4\epsilon^{\prime}}\log(N))$ which is an improvement over the cost of full inference for $k\!>\!1.$ Improvements in methods for sampling exact or approximate k-DPPs (e.g. recent bounds on mixing time [Hermon & Salez, 2019]) or bounds on other selection schemes for Nyström approximations directly translate to improved bounds on the computational cost of convergent sparse Gaussian process approximations through this framework.

4 Pointwise approximate posterior

In many applications, pointwise estimates of the posterior mean and variance are of interest. It is therefore desirable that the approximate variational posterior gives similar estimates of these quantities as the true posterior. Huggins et al. derived an approximation method for sparse GP inference with provable guarantees about pointwise mean and variance estimates of the posterior process and showed that approximations with moderate KL divergences can still have large deviations in mean and variance estimates. However, if the KL divergence converges to zero, estimates of mean and variance converge to the posterior values. By the chain rule of KL divergence [Matthews et al., 2016],

Therefore, bounds on the mean and variance of a one-dimensional Gaussian with a small KL divergence imply pointwise guarantees about posterior inference when the KL divergence between processes is small.

The proof is in appendix B. If $\epsilon\to 0,$ proposition 1 implies $\mu_{1}\to\mu_{2}$ and $\sigma_{1}\to\sigma_{2}.$ Using this and theorems 3 and 4, the posterior mean and variance converge pointwise to those of the full model using $M\ll N$ inducing features.

Related work

Statistical guarantees for convergence of parametric GP approximations [Zhu et al., 1997; Ferrari-Trecate et al., 1999], lead to similar conclusions about the choice of approximating rank. Ferrari-Trecate et al. showed that given $N$ data points, using a rank $M$ truncated SVD of the prior covariance matrix, such that $\lambda_{M}\ll\sigma_{n}^{2}/N$ results in almost no change in the model, in terms of expected mean squared error. Our results can be considered the equivalent for variational inference, showing that theoretical guarantees can be established for non-parametric approximate inference.

Guarantees on the quality Nyström approximations have been used to bound the error of approximate kernel methods, notably for kernel ridge regression [Alaoui & Mahoney, 2015; Li et al., 2016]. The specific method for selecting columns in the Nyström approximation plays a large role in these analyses. Li et al. use an approximate k-DPP, nearly identical to the initialization we analyze; Alaoui & Mahoney sample columns according to ridge leverage scores. The substantial literature on bounds for Nyström approximations [e.g. Gittens & Mahoney, 2013] motivates considering other initialization schemes for inducing points in the context of Gaussian process regression.

Conclusion

We proved bounds on the KL divergence between the variational approximation of sparse GP regression to the posterior, that depend only on the decay of the eigenvalues of the covariance operator. These bounds prove the intuitive result that smooth kernels with training data concentrated in a small region admit high quality, very sparse approximations. These bounds prove that truly sparse non-parametric inference, with $M\ll N,$ can provide reliable estimates of the marginal likelihood and pointwise posterior.

Extensions to models with non-conjugate likelihoods, especially within the framework of Hensman et al. , pose a promising direction for future research.

Acknowledgements

We would like to thank James Hensman for providing an excellent research environment, and the reviewers for their helpful feedback and suggestions. We would also 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.

References

Appendix A Proof Of Lemma 1

We can rewrite the second term (ignoring the factor of one half) in Appendix A as,

The last inequality comes from noting that the fraction in the sum attains a maximum when $\gamma_{i}$ is minimized. Since $\sigma_{n}^{2}$ is a lower bound on the smallest eigenvalue of ${\bf Q}_{n},$ we have,

Appendix B KL Divergence Gaussian Distributions

We make use of the formula for KL divergences between multivariate Gaussian distributions in our proof of lemma 2, and the univariate case in proposition 1.

Recall the KL divergence from p_{1}\sim\mathcal{N}\mathopen{}\mathclose{{}\left(\mathbf{m_{1}},\mathbf{S_{1}}}\right) to p_{2}\sim\mathcal{N}\mathopen{}\mathclose{{}\left(\mathbf{m_{2}},\mathbf{S_{2}}}\right) both of dimension $N$ is given by

The inequality is a special case of Jensen’s inequality.

B.2 Proof of Upper Bound in lemma 2

In order to complete the proof, we need to show that the second term on the right hand side is bounded above by $t/(2\sigma_{n}^{2})$. Using Equation 19:

The inequality follows from noting the log determinant term is negative, as ${\bf K}_{n}\succ{\bf Q}_{n}$ (i.e. ${\bf K}_{n}-{\bf Q}_{n}$ is positive definite). Simplifying the last term,

The first inequality uses that for positive semi-definite symmetric matrices $\operatorname{Tr}(AB)\leq\operatorname{Tr}(A)\lambda_{1}(B)$ which is a special case of Hölder’s inequality for Schatten norms. The final line uses that the largest eigenvalue of ${\bf Q}_{n}^{-1}$ is bounded above by $\sigma_{n}^{-2}.$ Using this in Equation 20 finishes the proof.

B.3 Proof of Proposition 1

where we have defined $x=\frac{\sigma^{2}_{1}}{\sigma^{2}_{2}}.$

Applying the lower bound $x-\log(x)-1\geq(x-1)^{2}/2-(x-1)^{3}/3,$

A bound on $\lvert x-1\rvert$ that holds for all $\epsilon$ can then be found with the cubic formula. Under the assumption that $\epsilon<\frac{1}{5},$ we have $x-\log(x)<1.2$ which implies $x\in[0.493,1.77].$ For $x$ in this range, we have

Using our bound on the ratio of the variances completes the proof of proposition 1.

Appendix C Covariances for Interdomain Features

We compute the covariances for eigenvector and eigenfunction inducing features.

Recall we have defined eigenvector inducing features by,

This is the $i^{th}$ entry of the matrix vector product $\mathbf{K}_{\bf ff}\mathbf{w}^{(m)}=\lambda_{m}(\mathbf{K}_{\bf ff})\mathbf{w}^{(m)}_{i}.$

C.2 Eigenfunction inducing features

Recall we have defined eigenfunction inducing features by,

The expectation and integration may be interchanged by Fubini’s theorem, as both integrals converge absolutely since $p(\mathbf{x})$ is a probability density, the $\phi_{m}(\mathbf{x})$ are in $L^{2}(\mathcal{X})_{p}\subset L^{1}(\mathcal{X})_{p}$ and $k$ is bounded.

We may then apply the eigenfunction property to the inner integral and orthonormality of eigenfunctions to the result yielding,

Appendix D Discrete k-DPPs

The first inequality follows from lemma 3. The second uses the triangle inequality replace $t(Z)$ with a bound on its maximum. The final line uses one of the definitions of total variation distance for discrete random variables.

D.2 Sampling Approximate k-DPPs

Belabbas & Wolfe proposed using the Metropolis method for approximate sampling from a k-DPP. Several recent works have shown that a natural Metropolis algorithm on k-DPPs mixes quickly. In particular, Anari et al. considers the following algorithm:

Let A denote algorithm 1. Let $\nu^{R}$ denote the distribution induced by $R$ steps of A. Let $R(\epsilon)$ denote the minimum $R$ such that $\|\mu-\nu^{R}\|_{TV}<\epsilon,$ where $\mu$ is a k-DPP on some kernel matrix $\mathbf{K}_{\bf ff}.$ Then

Taking $\epsilon$ to be any fixed inverse power of $N,$ (i.e. $\epsilon=N^{-\gamma},$ will make the second term $\mathcal{O}(NM\log(N)),$ while by taking $\gamma$ large (e.g. greater than 2), we can make $2Nv\epsilon$ small.

The computation of $L_{T}$ proceeds in two steps: first we compute $\mathbf{L}_{S\setminus i}$ using $\mathbf{L}_{S}.$

We now need to extend a Cholesky factorization from $S\setminus i$ to $T,$ which involves adding a row.

Appendix E Proof of Corollaries

From theorem 3, with probability $1-\delta$ and this choice of $\epsilon$

Take $M=\frac{(3+\gamma)\log(N)+\log D}{\log(B^{-1})}.$ If $M\geq N$ the KL-divergence is zero and we are done. Otherwise, $C(M+1)<N^{2}\sum_{i=M+1}^{\infty}\lambda_{i}.$ By the geometric series formula,

implying $\frac{C(M+1)}{2\delta\sigma_{n}^{2}}<N^{-1-\gamma}.$ Using this in section E.1 completes the proof.

E.2 Corollary 3

It is sufficient to consider the case of isotropic kernels and input distributions.For the general case, the eigenvalues can be bounded above by constant times the eigenvalues of an operator with an isotropic kernel with all lengthscales equal to the shortest kernel lengthscale and the input density standard deviation set to the largest standard deviation of any one-dimensional marginal of $p(\mathbf{x}).$ From [Seeger et al., 2008] in the isotropic case (i.e. $B_{i}=B_{j}=:B$ for all $i,j\leq D),$

In the second line, we use that $B<1,$ so $B^{s^{1/D}}$ obtains its minimum on the interval $s\in[s^{\prime},s^{\prime}+1]$ at the right endpoint (i.e. monotonicity). We now define $\alpha=-\log(B).$ So,

In the second line we made the substitution $t=\alpha s^{1/D},$ so $ds=\alpha^{-D}Dt^{D-1}.$ We now recognize,

as an incomplete gamma function, $\Gamma(D,\alpha M^{1/D}).$ From Gradshteyn & Ryzhik [2014, 8.352],

As $M$ grows as a function of $N$ and $D$ is fixed, for $N$ large $D\leq\alpha M^{1/D}.$ This implies that that the largest term in the sum on the right hand side is the final term, so

Choose M=\frac{1}{\alpha}\log\mathopen{}\mathclose{{}\left(N^{\gamma^{\prime}}\mathopen{}\mathclose{{}\left(\frac{2a}{A}}\right)^{D/2}D^{2}\alpha^{-1}}\right)^{D}. Then

For any fixed $D,$ for this choice of $M$ for any $\varepsilon>0$ for $N$ large,

By choosing $\gamma^{\prime}>3+\epsilon^{\prime},$ for some fixed $\epsilon^{\prime}>0$ the proof is complete, using a similar argument as the one used in the proof of the previous corollary.

Note that using the bound proven in Theorem 2 (the tightest of our bounds) the exponential scaling in dimension is unavoidable. If both $k$ and $p(\mathbf{x})$ are isotropic, then the eigenvalue $(\frac{2a}{A})^{D}B^{m}$ appears $\binom{m+D-1}{D-1}$ times. This follows from noting that this is the number of ways to write $m$ as a sum of $D$ non-negative integers. Using an identity and standard lower bound for binomial coefficients \sum_{i=1}^{K}\binom{m+D-1}{D-1}=\binom{K+D}{D}\geq\mathopen{}\mathclose{{}\left(\frac{K+D}{D}}\right)^{D}>C(D)K^{D}. for some constant depending on D, $C(D).$ Using Theorem 2 we need to choose $M$ such that $\lambda_{M}=\mathcal{O}(1/N).$ This means choosing $K\gg\log(N)$ in the sum above, leading to at least $\alpha\log^{D}(N)$ features being needed for some constant $\alpha$. The constant in this lower bound decays rapidly with $D,$ while the constant in the upper bound does not. Better understanding this gap is important for understanding the performance of sparse Gaussian process approximations in high dimensions.

If the data actually lies on a lower dimensional manifold, we conjecture the scaling depends mainly on the dimensionality of the manifold. In particular, if the manifold is linear and axis-aligned, then the kernel matrix only depends on distances along the manifold (not in the space it is embedded in) so the eigenvalues will not be effected by the higher dimensional embedding. We conjecture that similar properties are exhibited when the data manifold is nonlinear.

Appendix F Smoothness and Sacks-Ylivasker Conditions

In many instances the precise eigenvalues of the covariance operator are not available, but the asymptotic properties are well understood. A notable example is when the data is distributed uniformly on the unit interval. If the kernel satisfies the Sacks-Ylivasker conditions of order $r$:

$k(x,x^{\prime})$ is $r$-times continuously differentiable on $^{2}$ Moreover, $k(x,x^{\prime})$ has continuous partial derivatives up to order $r+2$ times on $(0,1)^{2}\cap(x>x^{\prime})$ and $(0,1)^{2}\cap(x<x^{\prime}).$ These partial derivatives can be continuously extended to the closure of both regions.

Let $L$ denote $k^{(r,r)}(x,x^{\prime})$, $L_{+}$ denote the restriction of $L$ to the upper triangle and $L_{-}$ the restriction to the lower triangle, then $L_{+}^{(1,0)}<L_{-}^{(1,0)}$ on the diagonal $x=x^{\prime}.$

$L^{(2,0)}_{+}(s,\cdot)$ is an element of the RKHS associated to $L$ and has norm bounded independent of $s.$

Notably, Matérn half integer kernels of order $r+1/2$ meet the S-Y condition of order $r.$ See Ritter et al. for a more detailed explanation of these conditions and extensions to the multivariate case.