Rates of Convergence for Sparse Variational Gaussian Process Regression

David R. Burt; Carl E. Rasmussen; Mark van der Wilk

Source-linked AI summary

Rates of Convergence for Sparse Variational Gaussian Process Regression

David R. Burt, Carl E. Rasmussen, Mark van der Wilk

arXiv:1903.03571v3Mar 8, 2019stat.MLcs.LG

TL;DR

The paper addresses how the number of inducing variables must grow with dataset size to maintain accurate sparse Gaussian-process approximations. It derives a priori KL-divergence bounds and shows that sublinear growth of inducing variables can make the approximation arbitrarily accurate, including M = O(log N) under stated conditions.

Problem
Existing results characterize computational cost but provide limited guidance on how the number of inducing variables must grow with dataset size to preserve approximation quality.
Method
The paper analyzes a priori upper bounds on KL divergence, relating their behavior to kernel, input-distribution, inducing-variable, and dataset-size assumptions.
Results
Sublinear growth of inducing variables makes the KL divergence tend to zero under intuitive assumptions, and M = O(log N) suffices under the stated squared-exponential-kernel conditions.
Takeaways & Limitations
Sparse GP inference with M ≪ N can provide reliable marginal-likelihood and pointwise-posterior estimates while offering an a priori computation-growth rule for continual learning.
Takeaways & Limitations
Extensions to non-conjugate likelihoods remain a direction for future research.

Abstract

from arXiv · show

Excellent variational approximations to Gaussian process posteriors have been developed which avoid the $\mathcal{O}\left(N^3\right)$ scaling with dataset size $N$. They reduce the computational cost to $\mathcal{O}\left(NM^2\right)$, with $M\ll N$ being the number of inducing variables, which summarise the process. While the computational cost seems to be linear in $N$, the true complexity of the algorithm depends on how $M$ must increase to ensure a certain quality of approximation. We address this by characterising the behavior of an upper bound on the KL divergence to the posterior. We show that with high probability the KL divergence can be made arbitrarily small by growing $M$ more slowly than $N$. A particular case of interest is that for regression with normally distributed inputs in D-dimensions with the popular Squared Exponential kernel, $M=\mathcal{O}(\log^D N)$ is sufficient. Our results show that as datasets grow, Gaussian process posteriors can truly be approximated cheaply, and provide a concrete rule for how to increase $M$ in continual learning scenarios.

1. Introduction

The paper addresses how the number of inducing variables must grow with dataset size to preserve sparse GP approximation quality. It bounds KL divergence and shows that sublinear growth of M can make this divergence vanish.

Motivation: Sparse GP methods reduce computational cost with M inducing variables, but their true scaling depends on how M must grow with N.The approximation quality and the costs of initializing or optimizing inducing points determine realistic scaling.
Approach: The paper measures approximation quality using the KL divergence between the variational approximation and the full posterior process.The variational framework minimizes this divergence from the approximate posterior to the posterior process.
Main result: Under intuitive assumptions, M needs to grow only sublinearly with N for the KL divergence to approach zero.The result supports sparse approximations for large datasets while retaining accurate predictions and uncertainty estimates.
Proof strategy: The proof bounds KL divergence through the error of a Nyström covariance approximation and relates that error to an infinite-dimensional operator spectrum.For stationary kernels, smoother sample functions and more concentrated inputs admit sparser approximations.
Assumptions and bounds: The main bounds assume i.i.d. training inputs and depend on a function g(M, N) that grows linearly in N and generally decays rapidly in M.The initialization quality determines an additional error term that can be made arbitrarily small at extra computational cost.

2. Background and notation

This section introduces GP regression, sparse variational approximations, and spectral tools for analyzing their KL error. It motivates a priori bounds that specify how inducing-variable counts should scale before data are observed.

Gaussian process regression: GP regression places a prior over functions and uses the posterior to predict new outputs while accounting for uncertainty from finite training data.The marginal likelihood is used for hyperparameter selection.
Computational motivation: Exact GP quantities become infeasible for large datasets because determinant and inverse computations have unfavorable time complexity.Low-rank approximations replace the full covariance with rank M, enabling reduced-cost computation.
Sparse variational approximation: The variational framework represents function values at M inducing inputs with a Gaussian q(u), whose parameters define the approximate posterior.The variational parameters are the inducing inputs and the mean and covariance of q(u).
Approximation effects: Joint ELBO maximization can bias hyperparameter estimation, including overestimation of the noise variance, although adding inducing points reduces the KL gap.The bias is practically eliminated when enough inducing variables are used.
Selecting inducing variables: An upper bound on the KL divergence provides a criterion for assessing approximation quality beyond checking whether the ELBO stops improving.A small difference between the upper and lower bounds guarantees a good approximation under the stated bound.
Spectral analysis: For large datasets, the covariance operator captures limiting covariance-matrix spectral properties under an input distribution and has a discrete spectrum when compact.Its eigenvalues and eigenfunctions support asymptotic analysis of approximation error.
A priori scaling: The paper focuses on a priori bounds describing asymptotic behavior as N grows and M varies with N, including guarantees relevant to continual learning.These bounds concern datasets satisfying the assumed conditions before the data are observed.

3. Bounds on the KL divergence for eigenfunction inducing features

The section bounds the variational KL divergence using spectral information and shows how eigenfunction features control approximation quality as datasets grow. These bounds imply kernel- and input-dependent growth rates for the number of inducing features, including logarithmic growth for the one-dimensional squared exponential case and sublinear growth for Matérn kernels.

The analysis bounds KL(Q||P̂) using spectral information about the covariance matrix or its associated integral operator.These a priori bounds underpin corresponding bounds for inducing points.
The KL upper-bound analysis depends on the Nyström residual trace t and largest residual eigenvalue, with t ≪ 1/N making KL small when ||y||2 = O(N).The bound is obtained by comparing refinements of variational lower and upper bounds.
Increasing N with fixed M requires more inducing points to ensure the sparse approximation converges, even when inputs lie in a fixed interval.For uniformly sampled inputs on [0, 5], this behavior is illustrated with 15 inducing points.
Eigenfunction inducing features are selected from eigenvectors of the training covariance matrix and produce the optimal rank-M approximation to Kff.Computing the first M eigenvalues and eigenvectors can be done using a truncated SVD, while the resulting residual is characterized by λM+1(Kff).
For one-dimensional squared exponential kernels with Gaussian inputs, M = O(log N) makes the KL bound tend to 0 as N tends to infinity.For one-dimensional Matérn k + 1/2 kernels, M = N^α with α > 1/(2k+1) suffices, which is sublinear for k > 0.

4. Bounds for inducing points

The paper extends KL-divergence bounds to inducing points, using determinant-based initialization and empirical input distributions. These bounds are only slightly weaker than interdomain approaches while remaining applicable to any kernel.

4. Bounds for inducing points: For squared-exponential kernels, determinant-based sampling produces more dispersed inducing points than uniform sampling for both tested lengthscales.Figure 2 compares ℓ = 2 and ℓ = .5 determinant-based samples with uniform samples.
4. Bounds for inducing points: Inducing points apply efficiently to any kernel, unlike the interdomain features whose practical use is limited by computational considerations and unavailable analytic quantities.The limitations include computational costs for eigenvector features, missing analytic expressions for K_uf, and unknown input density p(x).
4. Bounds for inducing points: Determinant-based sampling selects inducing-point subsets with probability proportional to the determinant of their principal kernel submatrix.This is the k-DPP construction used to initialize inducing points from the empirical covariance matrix.
4. Bounds for inducing points: M inducing points initialized by determinant-based sampling yield average bounds within a factor of M + 1 of eigenvector inducing features.The selection scheme also introduces negative correlations between inducing-point locations, making them well dispersed.
4. Bounds for inducing points: Approximate k-DPP sampling avoids eigendecomposition because exact sampling is computationally incompatible with the target sparse-GP cost.The method therefore uses cheaper “ϵ close” sampling procedures, with total variation distance controlling the approximation to the exact k-DPP.
4.1. A priori bounds for inducing points: Theorems 3 and 4 establish high-probability inducing-point bounds under i.i.d. inputs, bounded kernel diagonals, and, for Theorem 4, prior-generated targets.The bounds are derived through Nyström approximation error and the spectrum of the associated integral operator.

5. Consequences of theorem 3 and theorem 4

The results characterize how inducing-feature counts must grow for sparse GP approximations to converge, with rates determined by kernels and input distributions. For SE kernels with Gaussian inputs, polylogarithmic growth in N suffices, yielding substantially reduced computational complexity.

SE kernels and Gaussian inputs: The KL divergence converges as an inverse power of N when M = C log(N) + C0 for one-dimensional SE regression with Gaussian inputs.Figure 4 increases N while generating training outputs from the prior generative model.
SE kernels and Gaussian inputs: M = O(log^D N) inducing points suffice for SE-ARD kernels with Gaussian inputs in D dimensions, with KL convergence holding with high probability.The result assumes fixed ϵ′, δ > 0 and ϵ = O(N^-γ) for some fixed γ > 2.
Kernel and input dependence: Smooth kernels and concentrated input distributions admit sparser approximations because their covariance-operator eigenvalues decay more rapidly.For other kernels or input distributions, the inducing-feature count need not grow polylogarithmically in N.
Computational complexity: O(N log^(4D+1) N) computational complexity follows for SE kernels when k-DPP initialization uses M = O(log^D N) inducing points.This is compared with O(N^3) exact inference cost.
Computational complexity: For one-dimensional Matérn k + 1/2 kernels, M = O(N^(1/(2k)+ϵ′)) yields O(N^(1+2/k+4ϵ′) log(N)) complexity, improving on full inference for k > 1.The stated rate uses the average-case analysis of theorem 4.
Pointwise posterior quantities: Small process-level KL divergence implies pointwise convergence of posterior means and variances using M ≪ N inducing features.The implication follows from Gaussian mean-and-variance bounds together with the main convergence theorems.

6. Related work

The paper relates its variational-inference guarantees to earlier approximation results for parametric GP models and kernel methods. Its contribution is framed as a non-parametric analogue with Nyström-selection analysis.

Parametric GP approximations: Earlier parametric GP approximation results established related conclusions about the required approximating rank.A rank-M truncated SVD can produce almost no change in expected mean squared error when λ_M ≪ σ_n^2/N.
Parametric GP approximations: The paper interprets its guarantees as an equivalent result for variational inference, extending theoretical convergence guarantees to non-parametric approximate inference.The comparison is made with prior SVD-based parametric results.
Nyström approximations: Nyström approximation analyses for kernel methods show that the column-selection method strongly affects guarantees, including approximate k-DPP and ridge-leverage-score sampling.The paper analyzes an initialization nearly identical to the approximate k-DPP method used in related work.

7. Conclusion

The paper proves KL-divergence bounds for sparse variational GP regression based on covariance-operator eigenvalue decay. These bounds support reliable marginal-likelihood and pointwise-posterior estimates with M ≪ N, while non-conjugate likelihoods remain a future extension.

Main conclusion: The KL bounds depend only on the decay of covariance-operator eigenvalues and imply that smooth kernels with concentrated training inputs admit high-quality sparse approximations.The conclusion states that these approximations can use M ≪ N inducing variables.
Main conclusion: Sparse non-parametric inference with M ≪ N can provide reliable estimates of the marginal likelihood and pointwise posterior.This is the paper's stated practical consequence of the KL bounds.
Future work: Extensions to models with non-conjugate likelihoods, particularly within the Hensman et al. framework, are identified as future research.

A. Proof Of Lemma 1

The proof develops Gaussian KL bounds through spectral properties of the covariance approximation. It diagonalizes the relevant positive-definite matrix and bounds the resulting terms using its eigenvalues.

Spectral decomposition: Writing Q_n = UΛU^T diagonalizes the positive-definite matrix, whose eigenvalues satisfy γ_1 ≥ γ_2 ≥ ... ≥ γ_N ≥ σ_n^2.The matrix Q_n is positive definite because Q_ff is symmetric positive semidefinite.
Spectral decomposition: The change of variables z = U^T y preserves the Euclidean norm because U is unitary.
Bounding the KL terms: The proof bounds the relevant sum by minimizing the eigenvalue-dependent denominator, using σ_n^2 as a lower bound on the smallest eigenvalue.The inequality is applied after rewriting the second term in Equation 18.
Gaussian KL calculation: The proof uses the formula for KL divergence between multivariate Gaussian distributions, together with the univariate case for Proposition 1.
Bounding the KL terms: The inequality used in the bound is a special case of Jensen’s inequality.

B.2. Proof of Upper Bound in lemma 2

The proof bounds the second term using matrix inequalities and properties of the covariance matrices, completing the upper-bound argument.

The log-determinant term is nonpositive because K_n − Q_n is positive definite.
Tr(AB) ≤ Tr(A)λ_1(B) bounds the trace term for positive semidefinite symmetric matrices.This is identified as a special case of Hölder’s inequality for Schatten norms.
The largest-eigenvalue bound for Q_n^-1 completes the required estimate when substituted into Equation 20.

B.3. Proof of Proposition 1

The proposition is proved by bounding variance ratios and deriving covariance expressions for eigenvector and eigenfunction inducing features under stated integrability conditions.

A cubic-formula bound gives x ∈ [0.493, 1.77] when ϵ < 1/5.The bound follows from x − log(x) < 1.2 over the relevant range.
The resulting variance-ratio bound completes the proof of Proposition 1.
Eigenvector inducing-feature covariances are computed using the eigenvector relation K_ff w^(m) = λ_m(K_ff)w^(m) and orthonormality.
Eigenfunction inducing-feature covariances use Fubini’s theorem, the eigenfunction property, and eigenfunction orthonormality.The interchange is justified because the integrals converge absolutely under the stated boundedness and L2 conditions.

D.1. Proof of Corollary 1 from Lemma 3

The corollary combines total-variation control with a determinant-based inducing-point sampler whose initialization and updates admit explicit computational bounds.

The sampler proposes replacing one selected point with one outside the current set and accepts with probability 1/2 min{1, det(K_T)/det(K_Sr)}.
Choosing ϵ = N^-γ makes the second term O(NM log(N)), while sufficiently large γ makes 2Nvϵ small.
The total algorithm cost is greedy initialization plus R(ϵ) times the per-iteration cost.
Each iteration maintains a determinant and Cholesky factorization, extending the factorization after a point replacement.
Rank-one Cholesky updates and related calculations can be performed in O(M^2) time per iteration.

E. Proof of Corollaries

The corollaries derive feature-growth conditions from eigenvalue bounds, including logarithmic scaling for fixed dimension and explicit dimensional limitations.

The general case can be bounded using an isotropic kernel with the shortest lengthscale and the largest marginal input standard deviation.
The eigenvalue tail analysis uses the substitution t = αs^(1/D) and an incomplete gamma function.
For fixed D, the selected feature-growth choice yields the stated asymptotic bound for any ε > 0 and sufficiently large N.
The isotropic analysis requires at least α log^D(N) features for some constant α to achieve λ_M = O(1/N).The lower-bound constant decays rapidly with D, whereas the upper-bound constant does not.
If data lies on a lower-dimensional manifold, the authors conjecture that scaling depends mainly on manifold dimension rather than ambient dimension.For a linear axis-aligned manifold, the kernel matrix depends on distances along the manifold.
The one-dimensional uniform-input example assumes Sacks–Ylivasker smoothness conditions, which are satisfied by Matérn half-integer kernels of the corresponding order.