E.g. you have a large dataset and a kernel you consider suitable, but do not want to run another dense Cholesky factorisation. Someone suggests a sparse GP. This sounds reasonable, but we still need to decide what kind of approximation to use.

A stochastic variational GP (SVGP) and a Vecchia approximation can both reduce the computational cost. However, they do not discard the same information. Understanding the geometry of this difference helps explain when each approximation may be suitable.

For the discussion below, assume a zero-mean GP and, unless stated otherwise, independent Gaussian observation noise. The distinction also applies beyond regression, but this setting makes the algebra easier to follow.

What the inducing points are doing

Take M inducing inputs Z and their function values u = f(Z). The usual construction gives us

QXX = KXZ KZZ−1 KZX

as the covariance explained through those inducing values. In the kernel's feature space, this is the Gram matrix of projections onto the span of k(·, zj). Two inducing points may appear well separated in the input plot but contribute almost the same direction in that space. With a short lengthscale, even nearby points can provide different information.

This is sometimes described as replacing the GP by a rank-M process. For the standard variational construction, that description is not fully correct. We use q(f, u) = p(f | u)q(u), retaining the prior conditional p(f | u). If q(u) is Gaussian with mean m and covariance S, and A = KXZKZZ−1, then

Eq[f] = Am
Covq(f) = KXX − QXX + ASAT

The residual KXX − QXX is still there. What cannot be freely updated is the conditional distribution of the remaining function values given u. That restriction is where the approximation enters. Increasing M helps only if the added inducing variables capture something the existing representation misses. [1]

The “stochastic” part of SVGP concerns optimisation. A factorised likelihood lets us estimate the expected log-likelihood term using minibatches, while learning a shared q(u). This is useful for large datasets and non-Gaussian likelihoods, but minibatching does not remove the representation bottleneck. Optimising a small inducing set for longer is not the same as giving it more capacity. [2]

INDUCING REPRESENTATIONuu
RESTRICTED CONDITIONING123456
Shared inducing variables versus selected predecessor connections. Conceptual diagrams, not fitted models; the left diagram omits residual covariance.

With Vecchia, the sparsity is elsewhere

Order the observations and write their joint density as a product of conditionals. This is still exact. In the Vecchia approximation, we retain only a subset N(i) of the predecessors in each factor:

p(y) ≈ ∏i p(yi | yN(i))

For Gaussian observations, each factor is a small conditional regression. Its residual variance is Cii − CiNCNN−1CNi, where C includes the observation noise. Collecting these regressions gives a sparse triangular factor for the approximate precision. We have restricted the conditioning graph, not projected everything through a single small set of variables. There are latent-process and mixed variants too; the expression above is the observation-level version. [3]

A nearby observation is useful if conditioning on it makes omitted observations less informative. That is a stronger claim than saying nearby points have high correlation. With noisy observations, for example, one close neighbour may leave substantial uncertainty about the underlying function.

We also need to define what “nearby” means. Under an anisotropic kernel, distance in the original coordinates may not be appropriate. For an axis-aligned kernel with lengthscales ℓd, scaling each coordinate as xd/ℓd gives a more relevant notion of neighbourhood. Scaled Vecchia methods use this idea. The neighbourhoods then depend on covariance parameters being estimated, so we need to account for this in the implementation. [4]

The lengthscale relative to the domain

Consider observations spread over an interval. If the correlation length is comparable to the interval's width, a relatively small inducing representation may be sufficient. Now keep that lengthscale fixed and make the interval much wider. There are more regions that can vary almost independently. We should not expect the same number of inducing variables to describe them equally well.

So knowing that “the kernel is smooth” is not enough. The kernel spectrum over the input distribution, the dimension, and the accuracy required all influence how M needs to grow. Convergence results for sparse variational regression describe these dependencies under specific assumptions. They do not mean a few hundred inducing points will be sufficient for every dataset. [5]

A large spatial domain with short-range dependence is a reasonable setting to consider Vecchia. Adding observations in a distant region does not necessarily make every local conditional more complicated. This is useful if the selected neighbours make most omitted observations conditionally less informative. Smooth processes, noise, unusual sampling patterns or longer-range components can make this less straightforward.

SVGP is useful when minibatch training is important and a suitable inducing representation fits within the computational budget. Vecchia can also be used when the inputs are not geographic, and inducing methods do not have to be global. The distinction between these families is therefore not simply local versus global.

Comparing them on a new dataset

Fixed hyperparameters and an exact GP on a problem small enough to fit provide a useful initial comparison. Otherwise, a method that learns a longer lengthscale or a larger noise variance can appear to be a better approximation partly because the fitted problem has become easier. This is also worth studying, but it is a different comparison.

The next comparison allows both methods to fit the model, using several inducing budgets and conditioning-set sizes. Using the same number for both is not a fair comparison: M inducing variables and k neighbours represent different computational costs and approximation choices. Standard dense SVGP algebra typically involves O(B M2) minibatch work plus inducing-matrix operations that can cost O(M3); a straightforward Vecchia implementation uses many small conditional solves. Reuse, grouping, hardware and optimisation all change the practical cost. Measured runtime and peak memory therefore give a more useful comparison than these complexity expressions alone.

RMSE only describes part of the result. Predictive log density, coverage and interval width help show how the approximation handles uncertainty, especially in thinly sampled regions. For an application that integrates over a surface or uses joint draws, those quantities also need to be evaluated. Getting the marginal variances roughly right does not establish that the joint dependence is right.

The holdout also needs to reflect the intended application. If the task is to predict into a gap, a random split with training points close to almost every test point may give an optimistic result. A blocked holdout is more representative of that task.

No benchmark results are presented here. The choice depends on whether a small kernel-space representation is sufficient, or whether a small conditioning set retains most of the useful information. Either assumption can be wrong, but both can be tested on the dataset rather than accepted in advance.

References

  1. Titsias (2009). Variational Learning of Inducing Variables in Sparse Gaussian Processes. AISTATS.
  2. Hensman, Fusi & Lawrence (2013). Gaussian Processes for Big Data. UAI.
  3. Katzfuss & Guinness (2021). A General Framework for Vecchia Approximations of Gaussian Processes. Statistical Science.
  4. Katzfuss, Guinness & Lawrence. Scaled Vecchia approximation for fast computer-model emulation.
  5. Burt, Rasmussen & van der Wilk (2020). Convergence of Sparse Variational Inference in Gaussian Processes Regression. JMLR.

A different experience with either approach? Send a comment or correction.