Gaussian process (GP) priors are popular models for unknown functions in a variety of settings, including geostatistics (e.g., Stein, 1999;Banerjee et al., 2004;Cressie & Wikle, 2011), computer model emulation (e.g., Sacks et al., 1989;Kennedy & O'Hagan, 2001;Gramacy, 2020), and machine learning (e.g., Rasmussen & Williams, 2006;Deisenroth, 2010). Latent GP (LGP) models, such as generalized GPs, assume a Gaussian or non-Gaussian distribution for the data conditional on a GP (e.g., Diggle et al., 1998;Chan & Dong, 2011). LGPs extend GPs to a large class of settings, including noisy, categorical, and count data. However, LGP inference is generally analytically intractable and hence requires approximations. In addition, direct GP inference is prohibitive for large datasets due to cubic scaling in the data size. There are two main challenges for (L)GPs in many applications: One is to specify or learn a suitable kernel for the GP, and the other is carrying out fast inference for a given kernel. In this paper, we make no contributions to the former and instead focus on the latter challenge: We assume that a parametric kernel form is given and propose an efficient approximation method for LGP inference via structured variational learning.
Many approaches to scaling GPs to large datasets were reviewed in Heaton et al. (2019) and Liu et al. (2020), including low-rank approaches with a small number of pseudo points that are popular in machine learning. Such low-rank GP approximations have been combined with variational inference for GPs (e.g., Titsias, 2009;Hensman et al., 2013) and LGPs (e.g., Hensman et al., 2015;Leibfried et al., 2020).
A highly promising approach to achieve GP scalability is given by nearest-neighbor Vecchia approximations from spatial statistics (e.g., Vecchia, 1988;Stein et al., 2004;Datta et al., 2016;Katzfuss & Guinness, 2021), which are optimal with respect to forward Kullback-Leibler (KL) divergence under the restriction of sparse inverse Cholesky (SIC) factors of the covariance matrix (Schäfer et al., 2021a). Such SIC approximations have several attractive properties (e.g., as reviewed by Katzfuss et al., 2022). They result in a valid joint density function given by the product of univariate conditional Gaussians, each of which can be independently computed in cubic complexity in the number of neighbors. This allows straightforward mini-batch subsampling with unbiased gradient estimators (Cao et al., 2022). For the ordering and sparsity pattern used here, the number of neighbors needs to grow only polylogarithmically with the data size to achieve ϵ-accurate approximations for Matérn-type kernels up to boundary effects (Schäfer et al., 2021a) due to the screening effect (Stein, 2011). Many existing GP approximations, including low-rank and partially-independent conditional approaches, can be viewed as special cases of SIC approximations corresponding to particular orderings and sparsity patterns (Katzfuss & Guinness, 2021). SIC p(f ) p(f )
Prior p(f |y) p(f |y)
Posterior q(f ) q(f )
Bayes' theorem
Figure 1. Double KL minimization for approximating the posterior distribution of a latent Gaussian field f given data y: Based on a forward-KL-optimal SIC approximation p(f ) of the prior, we obtain an SIC-restricted reverse-KL-optimal variational approximation q(f ) to the posterior. approximation using our ordering and sparsity pattern does not exhibit the same limitations as low-rank approximations (Stein, 2014) and can hence be significantly more accurate for non-latent (i.e., directly observed) GPs (Cao et al., 2022).
LGPs are more challenging. For
LGPs with Gaussian noise, applying SIC approximations to the noisy responses reduces accuracy, and SIC approximations of the latent field may not be scalable (e.g., Katzfuss & Guinness, 2021). Existing approaches addressing this challenge (Datta et al., 2016;Katzfuss & Guinness, 2021;Schäfer et al., 2021a;Geoga & Stein, 2022) do not consider estimation using stochastic gradient descent (SGD). For non-Gaussian LGPs, Laplace SIC approximations (Zilber & Katzfuss, 2021) are straightforward but can be inaccurate. Liu & Liu (2019) combined an SIC-type approximation to the prior with variational inference based on a variational family of Gaussians with a sparse Cholesky factor of the covariance matrix, but we are not aware of results guaranteeing that the covariance-Cholesky factor exhibits (approximate) sparsity under random ordering. Wu et al. (2022) combined SIC-type approximations of LGPs with mean-field variational inference, but the latter may be inaccurate when there are strong correlations in the GP posterior (MacKay, 1992).
To achieve scalable and accurate inference for LGPs, we propose a variational family of SIC Gaussian distributions and combine it with a SIC approximation to the GP prior (see Figure 1). Our approach is double-KL-optimal in the sense that variational approximation is reverse-KL-optimal for a given log normalizer (i.e., evidence) and our prior SIC approximation, which is available in closed form, is forward-KL-optimal for a given sparsity pattern (Schäfer et al., 2021a). Within our double-Kullback-Leibler-optimal Gaussian-process framework (DKLGP), we then focus on a particular ordering and nearest-neighbor-based sparsity pattern resulting in highly accurate prior and posterior approximations. We adopt a novel computational trick based on the concept of reduced ancestor sets for achieving efficient and scalable LGP inference. For this setting, our variational approximation can be computed via SGD in polylogarithmic time per iteration. While inducing-point methods assume that unobserved points depend on data only through inducing points (e.g., Frigola et al., 2014;Hensman et al., 2015), our method allows fast and accurate KL-optimal prediction based on the screening effect. Our numerical comparisons show that DKLGP can be vastly more accurate than state-ofthe-art alternatives such as inducing-point and mean-field approximations at a similar computational complexity.
with
. Throughout, we view the inputs x i as fixed (i.e., non-random) and hence do not explicitly condition on them.
Unless y|f follows a Gaussian distribution, inference (such as computing the posterior p(f |y)) generally cannot be carried out in closed form. In addition, even for Gaussian likelihoods, direct inference scales as O(n 3 ) and is thus computationally infeasible for large n. To address these challenges, we propose an approximation based on double KL minimization.
Consider a lower-triangular sparsity pattern S q ⊂ {1, . . . , n} 2 , with {(i, i) : i = 1, . . . , n} ⊂ S q and such that i ≥ j for all (i, j) ∈ S q . Our preferred choice of S q will be discussed in Section 2.5, but typically we will have (i, j) ∈ S q if x i and x j are "close." Corresponding to S q , define the family of distributions
where we write V ∈ S q if (i, j) ∈ S q for all V ij ̸ = 0. It is straightforward to show that any q ∈ Q can be represented in ordered conditional form as q(f ) = n i=1 q(f i |f s q i ), where s q i = {j > i : (j, i) ∈ S q } for i = 1, . . . , n -1 and s q n = ∅. We approximate the posterior p(f |y) by the closest distribution in Q in terms of reverse KL divergence:
We have KL(q(f )∥p(f |y)) = log p(y) -ELBO(q), where p(y) does not depend on q, and so q satisfies q(f ) = arg max q∈Q ELBO(q).
(2)
Proposition 2.1. The ELBO in (2) can be written up to an additive constant of n/2 as
where L is the inverse Cholesky factor of K such that K -1 = LL ⊤ , and L :,i denotes its ith column.
All proofs can be found in Appendix C.
Even for a sparse V, computing the ELBO in (3) is prohibitively expensive for large n, because computing L (or any of its columns) from K generally requires O(n 3 ) time.
To avoid this, we replace the prior p(f ) defined in (1) by a Gaussian distribution that minimizes a second KL divergence under an SIC constraint.
Specifically, consider a second lower-triangular sparsity pattern S p ⊂ {1, . . . , n} 2 , which may be the same as S q . We define the corresponding set of distributions
We approximate the prior p(f ) by the closest approximation in P in terms of forward KL divergence: p(f ) = arg min p∈P KL p(f ) p(f ) . (4) By a slight extension of Schäfer et al. (2021a, Thm. 2.1), we can show that this optimization problem has an efficient closed-form solution. Proposition 2.2. The solution to (4) is p(f ) = N n (f |µ, ( LL ⊤ ) -1 ), where the nonzero entries of the ith column of L can be computed in O(|S p i | 3 ) time as
and S p i = {j : (j, i) ∈ S p } is an ordered set with elements in increasing order (i.e., the first element is i).
Throughout, we denote by e i a vector whose ith entry is one and all others are zero, and we index matrices before inverting so that K -1 S p i ,S p i := (K S p i ,S p i ) -1 . The approximation in Proposition 2.2 is equivalent to an ordered conditional approximation (Vecchia, 1988) of the prior density p
where
Plugging p(f ) into (2), the ELBO in (3) becomes
with the ith summand depending on L only via its ith column L:,i , whose nonzero entries can be computed in O(|S p i | 3 ) time using (5). We need to compute V -1 L:,i and V -1 e i , the latter of which appears in E q log p(y i |f i ) (see Section 2.6). The nonzero entry of e i (i.e., {i}) is a subset of the nonzero entries of L:,i (i.e., S p i ), and hence we focus our discussion on computing V -1 L:,i . Solving this sparse triangular system in principle requires O(|S q |) time.
However, it is possible to speed up computation by omitting rows and columns of V that do not correspond to the ancestor set A i of S p i with respect to S q , which is defined as A i = j ≥ i : there exists a path L = {(j, l 1 ), (l 1 , l 2 ), . . . , (l a-1 , l a ), (l a , l)} ⊂ S q for some l ∈ S p i . Ancestor sets are properties of the directed acyclic graphs that can be used to represent our triangular sparsity structures, as illustrated in Appendix B.
Proposition 2.3. (V -1 L:,i ) j = 0 for all j / ∈ A i .
Thus, we have
where V -1 Ai,Ai LAi,i can be computed in O(|A i ||S q i |) time.
Schäfer et al. (2021a) proposed a sparsity pattern S based on reverse-maximum-minimum-distance (r-maximin) ordering (see Figure 2 for an illustration). R-maximin ordering picks the last index i n arbitrarily (often in the center of the input domain), and then the previous indices are sequentially selected for k = n -1, n -2, . . . ,
, where I k = {i k+1 , . . . , i n }. Throughout, we assume that our indexing follows r-maximin ordering (e.g., f k = f i k ). We can then define the sparsity pattern by S i = {j ≥ i : dist(x i , x j ) ≤ ρℓ i }, for some fixed ρ ≥ 1, where ℓ i = min j>i dist(x i , x j ). We can compute dist(x i , x j ) as Euclidean distance between the inputs, potentially in a transformed input space (see Section 2.6 for more details). The conditioning sets are all of approximately size |S i | = O(ρ d ) ≈ m = |S|/n under mild assumptions on the regularity of the inputs. Schäfer et al. (2021a) proved that an ϵ-accurate approximation of the prior can be obtained using
Reverse maximin ordering on a grid (small gray dots) of size n = 60 × 60 = 3,600 on a square. For three different indices i, we show the ith ordered input (▲), the subsequently ordered n -i inputs ( ), the distance ℓi to the nearest neighbor (-), the neighboring subsequent inputs Si (■) within a (yellow) circle of radius ρℓi (here, ρ = 2), the reduced ancestors Ãi (+), and the ancestors Ai (×).
kernels K that are Green's functions of elliptic boundaryvalue problems (similar to Matérn kernels up to boundary effects) and demonstrated high numerical accuracy of the posterior using S q = S for Gaussian likelihoods. For non-Gaussian likelihoods, this implies highly accurate approximations to the posterior when a second-order Taylor expansion can adequately approximate the posterior.
While this means that our DKLGP can achieve high accuracy by choosing S p = S q = S, the resulting ancestor sets can grow roughly linearly with n (e.g., see Figure 3a). Hence, even evaluating the ELBO based on the ancestor sets would often be prohibitively expensive for large n. However, it is possible to ignore most ancestors in (7) and only incur a small approximation error. Specifically, consider reduced ancestor sets Ãi = {j ≥ i : dist(x i , x j ) ≤ ρℓ j }, where the last subscript is now a j, not an i. As illustrated in Figure 2, we have S i ⊂ Ãi (because ℓ j ≥ ℓ i for j ≥ i) and approximately Ãi ⊂ A i . The reduced ancestor sets are of size
and can all be computed together in O(nm log 2 n) time (Schäfer et al., 2021b). Hence, reduced ancestor sets can be orders of magnitude smaller than full ancestor sets (see Figures 3a and 6). Claim 2.4. For Matérn-type LGPs with exponential-family likelihoods, (V -1 L:,i ) j ≈ 0 for all j / ∈ Ãi , where V minimizes the ELBO in (6), under mild conditions.
We provide a non-rigorous justification for this claim in Appendix C. Together, Proposition 2.3 and Claim 2.4 imply that ∥V 3b), and so replacing the former by the latter in the ELBO causes negligible error (Figure 3c). Similar numerical results were obtained for two other popular kernels in Figures 7 and 8 in Appendix A, suggesting that our approach is applicable to beyond the Matérn family.
The class of distributions
We propose to find the optimal q ∈ Q by minimizing our approximation of -ELBO(q) with respect to these O(nm) unknown parameters via minibatch stochastic gradient descent. For each minibatch B, this requires computing the gradient of
using automatic differentiation.
For more general distributions p(y i |f i ), we can use the Monte Carlo gradient estimator (Kingma & Welling, 2014) and approximate E q log p(y
The O(m 3 ) cost dominates, as we typically need m = O(log d n) for accurate approximations (Schäfer et al., 2021a); for example, in Figure 3a, | Ãi ||S i | is smaller than |S i | 3 . Also, L does not need to be pre-computed and stored, as each column L:,i can be computed "on-the-fly"; this is especially useful for hyperparameter estimation, for which p(f ) and hence L changes with the hyperparameters at each gradient-descent iteration.
0 8000 16000 24000 32000 n 0 2000 4000 6000 8000 sparsity reduced ancestor full ancestor (a) Average sizes of Si, Ãi, and Ai 0.000 0.050 0.100 0.150 0.200 ||V 1 L:,i|| 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 Ãi , Ãi L Ãi ,i ∥ with reduced ancestor sets versus ∥V -1 L:,i∥ for i = 1, . . . , n, where n = 500 and d = 2. (c) compares ELBO curves based on full (6) and reduced (8) ancestor sets, as functions of the range parameter with true value 0.1, for n = 500 and d = 2. In all plots, we set ρ = 2 and the n inputs are sampled uniformly on [0, 1] d .
We initialize the optimization using an estimate of ν and V based on a Vecchia-Laplace approximation (Zilber & Katzfuss, 2021) of p(f |y) combined with an efficient incomplete Cholesky (IC0) approximation (Schäfer et al., 2021a) of the posterior SIC factor. While this initialization itself provides a reasonable approximation to the posterior, hyperparameter estimation for this approach is more difficult, and it is less accurate than DKLGP even for known hyperparameters as shown in Figure 9 in Appendix A.
The ordering and sparsity pattern in Section 2.5 depend on a distance metric, dist(x i , x j ), between inputs. We have found that the accuracy of the resulting approximation can be improved substantially by computing the Euclidean distance between inputs in a transformed input space in which the GP kernel is isotropic, as suggested by Katzfuss et al. (2022); Kang & Katzfuss (2023). For example, consider an automatic relevance determination (ARD) kernel of the form K(x i , x j ) = K o (q(x i , x j )), where K o is an isotropic kernel (e.g., a Matérn kernel with smoothness 1.5 is used throughout this paper) and q(x i , x j ) = ∥x λ ix λ j ∥ is a Euclidean distance based on scaled inputs x λ = (x 1 /λ 1 , . . . , x d /λ d ) with individual ranges or length-scales λ = (λ 1 , . . . , λ d ) for the d input dimensions. In this example, we take dist(x i , x j ) = q(x i , x j ) when computing the sparsity pattern. When the scaled distance and hence the sparsity pattern depend on unknown hyperparameters (e.g., λ in the ARD case), we carry out a two-step optimization procedure: First, we run our ELBO optimization for a few epochs based on the sparsity pattern obtained using an initial guess of λ to obtain a rough estimate of λ, which we then use to obtain the final ordering and sparsity pattern and warm-start our ELBO optimization.
An important task for (L)GP models is prediction at unobserved inputs, meaning that we want to obtain the posterior distribution of latent GP variables f * at new inputs x * 1 , . . . , x * n * given the data y. To do so, we consider the joint posterior distribution of f = (f * , f ), from which any desired marginal distribution can be computed. Since working with the joint covariance matrix K is again computationally prohibitive, we make a joint SIC assumption on the posterior distribution of f (with the prediction variables ordered first) that naturally extends the SIC assumption for f in q(f ). For the exact posterior, we have
Similarly, we assume q( f ) = q(f * |f )q(f ), where q(f ) = N n (f |ν, (VV ⊤ ) -1 ) was obtained as described in previous sections, and q(f * |f ) is a sparse approximation of p(f * |f ). For i = 1, . . . , n * , let S * i ⊂ {i, i + 1, . . . , n * + n} denote the ith sparsity set relative to the joint posterior.
We define the approximation to the joint posterior by the minimizer of the expected forward-KL divergence between p(f * |f ) and q(f * |f ) for given ν and V, that is,
The resulting approximation can be obtained efficiently: Proposition 2.5. For given ν, V, and
and
The posterior distribution of a desired summary, say a ⊤f can then be computed as q(a ⊤f ) = N (a ⊤ ν, ∥ Ṽ-1 a∥ 2 ).
In particular, the marginal posterior of f * i can be obtained using a = e i as q(e ⊤
). We again consider an r-maximin ordering and nearestneighbor sparsity pattern similar to above, but now conditioned on the prediction points being ordered first, and the training points ordered after (in the same ordering as before). Once the prediction points are in this conditional r-maximin ordering, we can define
and
This ordering and sparsity pattern can be computed rapidly and was shown to lead to highly accurate approximations; more details can be found in Schäfer et al. (2021a, Section 4.2.1). Note that while computing the prediction variances can be expensive, we can again approximate ∥ Ṽ-1 e i ∥ ≈ ∥ Ṽ-1
where the last subscript is a j, not an i.
We compared the following approaches:
DKLGP: Our method with r-maximin ordering and nearest-neighbor sparsity pattern DKL-G: Same as DKLGP but with global sparsity pattern S p i = S q i = {1, . . . , m} DKL-D: Same as DKLGP but with diagonal sparsity pattern S q i = {i} SVIGP: Stochastic variational GP proposed by Hensman et al. (2013) VNNGP: Variational nearest neighbor GP proposed by Wu et al. (2022) In figures and tables, we use abbreviated acronyms DKL, SVI, and VNN to save space. SVIGP and VNNGP are two state-of-the-art variational GP methods, while DKL-G and DKL-D are variants of our DKLGP that resemble SVIGP and VNNGP, respectively. SVIGP assumes independence in f conditional on m global inducing variables. VNNGP scales up the inducing points to be equal to the observed input locations, ensuring computational feasibility by assuming that each conditions only on m others a priori, combined with a mean-field approximation to the posterior. We used the GPyTorch (Gardner et al., 2018) implementations of SVIGP and VNNGP. For DKL-G and DKL-D, one can easily see that A i = S p i , and so reduced ancestor sets are not necessary. For all methods, computing a term in the ELBO requires O(m 3 ) time per sample. (Reusing Cholesky factors for all samples in a minibatch is straightforward for SVIGP; similar savings may also be possible for the other methods based on the supernode ideas suggested by Schäfer et al., 2021a.) Hence, m can be viewed as a comparable complexity parameter that trades off computational speed (for small m) against accuracy (large m). Thus, for our numerical comparison, we aligned the m for all methods with the average size of S i for a given ρ.
Throughout, we assumed f (•) ∼ GP(0, K), where K is a Matérn1.5 ARD kernel whose variance (set to one for simulations) and range (i.e., length-scale) parameters λ were estimated. We considered three different likelihoods p(y i |f i ):
The noise variance σ 2 ϵ was estimated from the data; for simulations, we used σ 2 ϵ = 0.1 2 except where specified otherwise.
For estimation of hyperparameters, the initial values for λ, σ 2 ϵ , and the variance in K were all 0.25. DKLGP and its variants ran the Adam optimizer for 35 epochs. SVIGP and VNNGP used natural gradient descent and Adam, respectively, as their optimizer for 500 epochs as suggested in Wu et al. (2022). The minibatch size was 128 and a multistep scheduler with a scaling factor of 0.1 was used for all methods.
Figure 4 provides a visual comparison of SVIGP, VNNGP, and DKLGP predictions for a toy example in one dimension. We also included predictions from the exact GP (DenseGP) which cannot be obtained for large n. DKLGP approximated the DenseGP most closely, especially in terms of
0.0 0.2 0.4 0.6 0.8 1.0 -4 -3 -2 -1 0 1 2 0.64 0.66 0.68 0.70 0.72 0.74 0.76 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 y f DKL SVI VNN DenseGP the prediction intervals. SVIGP oversmoothed heavily and produced very wide prediction intervals. VNNGP assumes a diagonal covariance in the variational distribution q(f ), which appears to have caused sharply fluctuating predictions and narrow prediction intervals. Figure 11 in Appendix A shows similar comparisons for Student-t and Bernoulli likelihoods.
We also carried out a more comprehensive comparison for 10,000 inputs randomly distributed in the unit hypercube, [0, 1] 5 , with true range parameters λ = (0.25, 0.50, 0.75, 1.00, 1.25). We used n = 8,000 inputs for training and 2,000 for testing. Performance was measured in terms of the variational inference of the latent field f (•) at training and test inputs. For each scenario, results over five replicates were produced and averaged.
Figure 5 compares root mean squared error (RMSE) and negative log-likelihood (NLL) at test inputs. For the Gaussian and Student-t likelihoods, DKLGP produced the most accurate predictions, while for the Bernoulli-logit likelihood, SVIGP and DKLGP appeared to be similarly accurate in terms of RMSE. DKLGP outperformed the competing methods in terms of NLL. While DKLGP, DKL-G, and SVIGP improved with increasing ρ as expected, the mean-field approximations (VNNGP and DKL-D) generally did not. We performed the same comparison for the squared exponential and rational quadratic kernels in Figures 12 and 13 in Appendix A, which resulted in the same rankings as for the Matérn kernel, except that DKLGP was marginally outperformed by SVIGP in terms of RMSE for the Bernoulli-logit likelihood at ρ = 2.0.
We also computed RMSE and NLL scores at training inputs for the methods we considered in Figure 5, as presented in Figure 10 in Appendix A. Consistent with the results from Figure 5, DKLGP generally performed best, which is consistent with the results from Figure 5. VNNGP performed similarly to DKLGP for the Gaussian and Student-t likelihoods, but underestimated the variance at test inputs and so led to poor NLL scores. Note that variational methods are generally known to underestimate the posterior variance (Blei et al., 2017).
To provide a more comprehensive comparison of SVIGP, VNNGP and DKLGP, we considered datasets from the UCI data repository widely used for benchmarking purposes. For the UCI datasets we considered in this section, covariates were first standardized to [0, 1] and removed from analysis if the standard deviation after standardization was smaller than 0.01. Furthermore, inputs were filtered to ensure that the minimum distance between inputs was greater than 0.001 to prevent numerical singularity. Approximately 20% of each dataset was used for testing. We chose different ρ for different datasets and computed the corresponding m. Since Section 3.3 demonstrated the advantage of DKLGP over its variants DKL-G and DKL-D, we excluded the two variants here for ease of presentation. We included SVIGP with m = 32 and m = 512 inducing points as benchmarks for easier comparison with relevant works in the literature.
Table 1 summarizes the performance of the three methods across nine UCI datasets. DKLGP had better scores than VNNGP for all datasets except for COVTYPE, for which VNNGP ran out of memory on a 64GB node despite having reduced the data size to a subset of size 100K. Compared to the SVIGP with similar computation cost, DKLGP provided substantially better performance for the binary response data COVTYPE and for low-dimensional (d < 10) settings, and roughly similar performance for most high-dimensional datasets except the KEGGU data, for which SVIGP produced much lower RMSE than DKLGP. However, this does not appear to be due to DKLGP providing a less accurate approximation to the exact GP, but rather it appears to be due to the exact GP (with its simple ARD kernel) being severely misspecified for KEGGU. To explore this further, we fitted the exact GP (DenseGP) to KEGGU. The DenseGP's RMSE was 0.14 (same as for DKLGP), and the root average squared distance between the DenseGP predictions and the DKLGP and SVIGP predictions was 0.05 and 0.13, respectively, which implies that the DKLGP predictions were a much better approximation of the exact-GP predictions than the SVIGP predictions. VNNGP provided better point predictions than SVIGP for the low-dimensional datasets, which is consistent with the results in Wu et al. (2022); however, VNNGP's NLL was high due to its underestimation of posterior variance.
Table 2 summarizes the wall-clock times on an Intel Xeon E5-2680 v4 CPU with 14 cores and 28 threads for the methods under comparison, where the computation of sparsity and ancestor sets is only applicable to DKLGP or DKL. The DKL computation times were closer to those of SVI than to those of SVI 32 , indicating that DKLGP and SVIGP should be compared at the same m on the basis of comparable computation times. While increasing m significantly improved SVIGP's performance on low-dimensional datasets, even m = 512 inducing points made the training of SVIGP challenging on the workstation we used for comparison. The performance of DKLGP can also be improved by using a larger ρ; for example, DKLGP's RMSE for the KIN40K data was reduced to 0.27 for m = 21.
We have introduced a variational approach using a variational family and approximate prior based on SIC restrictions. The (r-)maximin ordering, nearest-neighbor sparsity pattern, and a computational trick based on reduced ancestor sets together result in efficient and accurate inference and prediction for LGPs. While the time complexity is cubic in the number of neighbors, quadratic complexity for the prior approximation can be achieved by grouping observations and re-using Cholesky factors (Schäfer et al., 2021a); we will investigate an extension of this idea to computing the ELBO in our variational setting. Although we here assume that the input domain is Euclidean, our method can be applied more generally; using a correlation-based distance instead of Euclidean distance (Kang & Katzfuss, 2023), one can use our method to perform LGP inference for large data on complex domains (cf. Tibo & Nielsen, 2022). We will also explore extensions to deep GPs (cf. Sauer et al., 2022). An implementation of our method, along with code to reproduce all results, is publicly available at https: //github.com/katzfuss-group/DKL-GP.
Our approach is applicable to irregularly spaced observations and in principle to any desired covariance structure.
Our method provides state-of-the-art performance when fine-scale structure in the function of interest can be discerned from the data; in contrast, if the data are highly noisy or sparse or the covariance model is severely misspecified, inducing-point methods such as SVIGP that produce smooth predictions and wide uncertainty intervals may be competitive with our approach.
This section contains additional figures not shown in the main paper. Complementing Figure 3a, Figure 6 shows that reduced ancestor sets Ãi are much smaller than full ancestor sets A i across a range of ρ values. Complementing Figure 3, Figures 7 and 8 show that the approximation error in computing the ELBO caused by using reduced ancestor sets is negligible even for the squared-exponential and rational-quadratic kernels, respectively. Figure 9 suggests that the initialization of ν using Vecchia-Laplace approximation and incomplete Cholesky (IC0) approximation provides reasonable starting values for ν, which can be further refined by optimizing the ELBO.
Figure 10 shows a comparison of RMSE and NLL scores for the posterior marginals of the entries of f at training inputs. In contrast to Figure 5, VNNGP performed similarly to DKLGP and outperformed SVIGP for Gaussian and Student-t likelihoods. Furthermore, the Vecchia-Laplace approximation with IC0 (used as the initialization for DKLGP) was usually the third best model, indicating an advantage of using the SIC restriction for L and V. 0.000 0.100 0.200 0.300 0.400 0.500 ||V 1 L:,i|| 0.000 0.100 0.200 0.300 0.400 0.500 ||V 1 i, i L i, i|| (a) ∥V -1 Ãi , Ãi L Ãi ,i ∥ vs ∥V -1 L:,i∥ 0.05 0.10 0.15 0.20 lengthscale (range) 400 500 600 700 800 ELBO full reduced (b) ELBO using reduced ancestors Figure 8. The rational-quadratic-kernel versions of Figures 3b (left) and 3c (right): The left figure compares ∥V -1 Ãi , Ãi L Ãi ,i ∥ with reduced ancestor sets versus ∥V -1 L:,i∥ for i = 1, . . . , n, where n = 500 and d = 2. The right figure compares ELBO curves based on full (6) and reduced (8) ancestor sets, as functions of the range parameter with true value 0.1, for n = 500 and d = 2. In all plots, we set ρ = 2 and the n inputs are sampled uniformly on [0, 1] d . 2 0 2 true posterior mean 2 0 2 initial posterior mean 2 0 2 true posterior mean 2 0 2 est'd posterior mean
We here illustrate the sparsity and ancestor sets, using their graph representations. As pointed out by Katzfuss & Guinness (2021), the sparsity patterns can be represented by directed acyclic graphs (DAGs), which also allows straightforward visualization of ancestor sets. Figure 14 presents sparsity and ancestor sets for three selected points (i = 12, 4, 1) of 16 grid points in the unit square. For example, x 1 = 1 3 , 1 and x 16 = 2 3 , 2 3 . One can easily see that
The edges of the graphs corresponding to the ancestor sets A 12 , A 4 and A 1 are denoted by the black curved arrows. Specifically, the sparsity set S 1 = {2, 7, 13}, the reduced ancestor set Ã1 = S 1 ∪ {9, 11, 12} and the (full) ancestor set A 1 = Ã1 ∪ {15, 16}. Note that A 1 contains Ã1 , which is a desirable property for leveraging the screening effect in GPs (Stein, 2011;Bao et al., 2020). This is not always the case for small-scale problems and it depends on distribution of the points, as shown in Figure 14b. Specifically, A 4 = {10, 11, 14, 15, 16}, but Ã4 \ A 4 = {13} ̸ = ∅. But our numerical studies suggest that Ãi \ A i are typically empty or very small for large-scale problems, for which computational issues are severe and hence our method is most likely to be used. For relatively large i = 12, S 12 = Ã12 = A 12 = {16}. As illustrated here, all the reduced ancestor sets include x 16 , since ℓ 16 = ∞. Otherwise, unlike Ã4 and Ã1 , Ã12 does not include x 15 since dist(x 15 , x 12 ) = √ 2 is larger than ρℓ 15 ≃ 1.226. 0.0 0.2 0.4 0.6 0.8 1.0 -4 -3 -2 -1 0 1 2 0.0 0.2 0.4 0.6 0.8 1.0 -3 -2 -1 0 1 2 y f DKL SVI VNN DenseGP Figure 11. Comparison of three variational approximations to the predictive GP posteriors for the Student-t (left) and Bernoulli-logit likelihoods (right). We show the means (solid lines) and 95% pointwise intervals of the posterior predictive distribution f * |y at 200 regularly spaced test inputs. Note that the noise variance σϵ = 0.3 and range (or length-scale) λ = 0.1 are used. The exact-GP result (DenseGP) is available only for the Gaussian likelihood.
This section contains the postponed proofs of technical statements in the main paper. A non-rigorous justification for Claim 2.4 can be also found here.
Proof of Proposition 2.1. We have
where E q log p(y|f ) = n i=1 E q log p(y i |f i ). Using a well-known expression for the KL divergence between two Gaussian distributions, we have
where log
Proof of Proposition 2.2. Using a well-known formula for the KL divergence between two Gaussian distributions (e.g., see ( 9)), we have The ith ordered input (▲), the subsequently ordered n -i inputs ( ), the distance ℓi to the nearest neighbor (-), the neighboring subsequent inputs Si (■) within a (yellow) circle of radius ρℓi, with ρ = 1.3, the reduced ancestors Ãi (+), and the ancestors Ai (×). The directed acyclic graphs of the sparsity patterns are denoted by arrows (↷ ↷ ↷).
claim is that for j
By the results on exponential screening in Schäfer et al. (2021b), the matrix L satisfies the above decay property for covariances that are Green's functions of elliptic PDEs. It satisfies even the stronger property with ℓ j replaced by ℓ i .
For a Gaussian likelihood, the matrix V satisfies
where R is a diagonal covariance matrix of the likelihood. Interpreted as a PDE, the diagonal matrix R -1 corresponds to a zero-order term. Thus, the associated covariance matrix ( LL ⊤ ) -1 behaves like a discretized elliptic Green's function and is therefore subject to an exponential screening effect (Schäfer et al., 2021a, Section 4.1). Let P ↕ denote the permutation matrix that reverts the order of the degrees of freedom. Since P ↕ V -⊤ P ↕ is lower triangular and P ↕ ΣP ↕ = P ↕ V -⊤ P ↕ P ↕ V -1 P ↕ = P ↕ V -⊤ P ↕ P ↕ V -⊤ P ↕ ⊤ , the matrix P ↕ V -⊤ P ↕ is the Cholesky factor of Σ in the maximin (as opposed to the reverse maximin) ordering. In Schäfer et al. (2021b), it is shown that the Cholesky factors of discretized Green's functions of elliptic PDEs in the maximin ordering have exponentially decaying Cholesky factors. In particular, the results of Schäfer et al. (2021b) suggest that ∀j ≥ i : log
As shown, for instance, in Schäfer et al. (2021b, Lemma 5.19), products of matrices that decay rapidly with respect to a distance function dist(•, •) on its index set, inherit this decay property. To this end, assume that lower triangular matrices A and B satisfy this property. We then have log (AB) ji = log By the triangle inequality, we have dist (x j , x k ) + dist (x k , x i ) ≥ dist (x j , x i ). Since the right hand is -∞ unless j > i and thus ℓ j ≥ ℓ i , we have thus log (AB) ji = log k A jk B ki ⪅ log(n)dist (x j , x i ) /ℓ j , proving the the result.
For a general exponential family likelihood, the matrix V does not necessarily satisfy (10). Instead, according to Nickisch & Rasmussen (2008), a quadratic approximation to the log-likelihood under mild conditions implies that
where W is the covariance of the effective likelihood obtained by dividing the approximate posterior by the prior. Assuming that W -1 corresponds to a zero-order term in the context of a PDE, one can also obtain the result from the justification for the Gaussian likelihood case above.
Proof of Proposition 2.5. Note that p(f * |f ) = p( f )/p(f ) = N n * µ * + K * o K -1 (fµ), K * |o , where K * |o = K * * -K * o K -1 K o * , and q(f * |f ) = q( f )/q(f ) = N n * ν * -(V * * ) -⊤ V o * ⊤ (fν), (V * * V * * ⊤ ) -1 . Then, since KL p(f * |f ) q(f * |f ) is a KL divergence between two Gaussian distributions, we have
Using the fact that the first term is quadratic in form, one can show that
Then, we can see that KL p(f * |f ) q(f * |f ) is minimized with respect to ν * by Gµ + h = 0. This implies that ν * = µ * -(V * * ) -⊤ V o * ⊤ (νµ). Plugging this in, we have From the definition of S * i , it can be easily shown that K-1