Title: LOCALIZING AND AMORTIZING: EFFICIENT INFERENCE FOR GAUSSIAN PROCESSES

Abstract: The inference of Gaussian Processes concerns the distribution of the underlying function given observed data points. GP inference based on local ranges of data points is able to capture fine-scale correlations and allow fine-grained decomposition of the computation. Following this direction, we propose a new inference model that considers the correlations and observations of the K nearest neighbors for the inference at a data point. Compared with previous works, we also eliminate the data ordering prerequisite to simplify the inference process. Additionally, the inference task is decomposed to small subtasks with several technique innovations, making our model well suits the stochastic optimization. Since the decomposed small subtasks have the same structure, we further speed up the inference procedure with amortized inference. Our model runs efficiently and achieves good performances on several benchmark tasks.

1. INTRODUCTION:

Gaussian processes (GP) (Rasmussen & Williams, 2006) are flexible non-parametric models with a wide range of applications. GP poses a Gaussian prior over function values f and assumes observations y are generated independently given f . GP inference considers the calculation of the posterior of these function values (Matthews et al., 2016) given observations, namely p(f |y). Direct computation of the posterior is often intractable on large datasets, motivating people to consider its approximations. Variational inference (Jordan et al., 1999; Blei et al., 2017) for GP (Rasmussen & Williams, 2006) has achieved great successes recently. Variational inference constructs a variational distribution, which is usually a multivariate Gaussian distribution, to approximate the posterior. The approximation is done by minimizing the KL divergence from the posterior to the variational distribution (Blei et al., 2017) . The variational distribution is often constructed with some special structures to reduce the number of variational parameters and speed up the computation.

Inducing-point methods (Quiñonero-Candela & Rasmussen, 2005; Titsias, 2009; Hensman et al., 2013; 2015) define variational distributions on a small number M of inducing points and then derive the distribution of non-inducing points conditioned on these inducing points. Inducing points summarize the entire posterior distribution, and their number M balances the computational cost and the quality of the approximation. Inducing-point methods are further improved in several directions, such as generic inference for non-Gaussian likelihoods (Sheth et al., 2015; Dezfouli & Bonilla, 2015; Krauth et al., 2016; Hensman et al., 2015) , inter-domain and subspace inducing points (Hensman et al., 2017; Panos et al., 2018) , and decoupled approximation with two different sets of inducing points (Cheng & Boots, 2017; Salimbeni et al., 2018) . Burt et al. (2019) provide theoretical analysis to show that a relatively small M is sufficient to produce a reliable variational approximation when the input dimension is low.

While inducing-point methods capture global correlations among data points through inducing points, inference methods based on local neighbors focus more on correlation structures at local scales. These methods consider only local-range dependencies to save computation because localrange correlations are often much stronger than distant ones. Nguyen-Tuong et al. (2009) ; Park & Apley (2018) partition the input space into subregions, fit local models over subregions and then stitch local models into one. Other works examine neighbors of each data point directly. Gramacy & Apley (2015) investigate the properties of GP predictive equation and construct a local predictive approximator. Covariance tapering (Furrer et al., 2006; Kaufman et al., 2008) gains computational efficiency by constructing a sparse correlation matrix with zero correlations between distant data points. Methods based on Vecchia's approximation (Vecchia, 1988; Datta et al., 2016; Liu & Liu, 2019; Finley et al., 2019) decompose the joint probability of data points into conditionals according to a data ordering and then neglect far data points that are conditioned on.

Recently, Liu & Liu (2019) propose the AIGP method, which extends the idea of local inference to GP models with non-Gaussian likelihoods. They use directed graphical models to approximate both the prior and the posterior. With this construction, the inference task decomposes into local inference subtasks, then they introduce amortized inference and use inference networks to identify solutions to these subtasks (Kingma & Welling, 2013; Dai et al., 2015; Miao et al., 2016) . Amortization reduces the number of optimization parameters and greatly speeds up the inference procedure. However, this method has two drawbacks. First, the inference at a data point considers a few of its nearest neighbors but not all of them; therefore, it may lose some important correlations. Second, it depends on a data ordering. A bad ordering often deteriorates the performance, but it is hard to guard against such a bad situation. There are no easy fixes of the two issues, because all these designs in AIGP serve the purpose of decomposition.

In this work, we propose a new GP inference method, Localized and Amortized Inference based on Nearest neighbors (LAIN). LAIN considers K nearest neighbors for the inference at each data point. Particularly, LAIN uses a variational distribution whose covariance is parameterized by a sparse decomposition. The decomposition focuses on the correlations between every data point and its K nearest neighbors1 . LAIN also eliminates the need for a data ordering. These nice properties come after several technical innovations. First, the new distribution does not admit a decomposable entropy calculation. We overcome this difficulty by using a decomposable lower bound of the entropy (Ranganath et al., 2016; Louizos & Welling, 2017) . Second, to decompose the logarithm of the prior, AIGP and previous methods use a directed graphical model as an approximation of the prior. We follow this idea, but we consider all possible orderings of data points and collapse them to local combinations, making the computation manageable. With these techniques, LAIN still decomposes the inference task into subtasks, so amortized inference can apply. It is worthing noting that subtasks in LAIN are generated from the same mechanism while those in AIGP are not. We argue that subtasks sharing the same "distribution" are more appropriate for amortization.

Our empirical evaluations show that the LAIN method outperforms baseline methods including AIGP in several learning tasks. Our investigation also indicates that LAIN can achieve decent performance even only a few neighbors are considered.

2. BACKGROUND:

Gaussian Processes. Suppose we have a dataset containing a feature matrix X = (x i ) N i=1 and observations y = (y i ) N i=1 . We assume there is a latent function f that generates y i from x i for each i. Particularly, each y i is generated by a likelihood model p(y i |f i ) with

f i = f (x i ). Denote f = (f i ) N i=1 , then p(y|f ) = N i=1 p(y i |f i ).

The likelihood p(y i |f i ) can be very general -here we only assume that log p(y i |f i ) is differentiable with respect to f i . This mild assumption allows a wide range of data distributions. For example, if y i is binary, p(y i |f i ) is a Bernoulli distribution with f i as the logit.

We put a GP prior with a mean function ν(•) and a kernel function κ(•, •) over the latent function f . The kernel function encodes the prior knowledge of the smoothness of f . One commonly used kernel function is the Radial Basis Function (RBF) kernel, κ(x i , x j ) = r 2 exp(-0.5 x i -x j 2 2 /σ 2 ), with r and σ as parameters. With this prior, function values in f follow a multivariate Gaussian, f ∼ N (ν, Σ), with the mean ν = (ν(x i )) N i=1 and the covariance matrix Σ with Σ i,j = κ(x i , x j ) ∀i, j. GP inference concerns the calculation of the posterior p(f |y) (Matthews et al., 2016) , from which we can infer the function value f for any new input x with integral f p(f |f )p(f |y)df . The posterior p(f |y) is generally not tractable, so we appeal to approximate inference.

Variational Inference for GP. Variational inference approximates the posterior p(f |y) with a variational distribution q(f ), which is defined as a multivariate Gaussian distribution, q(f ) ∼ N (µ, V). The inference is carried out by maximizing the Evidence Lower BOund (ELBO) with respect to q(f ) (Blei et al., 2017) .

log p(y|X) ≥ max q(f ) E q [log p(y|f )] L ell + E q [log p(f )] Lcross -E q [log q(f )] Lent (1)

Here we name the three terms in the ELBO for easy reference later. Typically the ELBO is maximized by gradient-based optimization, preferably stochastic gradient optimization when N is large. Direct optimization of the ELBO is challenging, since the kernel matrix Σ and the variational covariance V are both large and have size N × N .

Inducing-point methods define q(f ) = f I q(f I )p(f |f I ) df I , where q(f I ) is the distribution over inducing points I, and p(f |f I ) is derived from the prior. The computation is reduced mainly because only the small distribution q(f I ) is optimized, while the conditional p(f |f I ) is fixed when the prior is given.

AIGP parameterizes V by a Cholesky decomposition, V = LL . Here L is a sparse lower triangular matrix, and each row of L has at most K non-zero entries. AIGP uses a triangular L for easy entropy computation. It also approximates log p(f ) with a directed graphical model. Both the lower triangular matrix L and the directed graph require an ordering of data points.

3. METHOD:



3.1. THE VARIATIONAL DISTRIBUTION:

Following previous works, we also define the variational distribution q(f ) to be a multivariate Gaussian N (µ, V). We parameterize V = RR + δ 2 I with R being a sparse matrix and δ being a small constant. Note that we do not require R to be triangular. The sparse pattern of R is decided by the nearest neighbors: R ij = 0 only when j ∈ n(i).

Here n(i) is the neighbor set containing data points that have the largest covariance with i in the prior (by definition n(i) includes i). In this work, we fix the size of n(i) to be K, though our derivation works for varied sizes of n(i). The row R i can be viewed as a representation of f i in the variational distribution: R i informs f i 's correlation with other function values, just like a word embedding informs its relation with other words (Mikolov et al., 2013) .

Efficient sampling from the marginal is critical for the decomposition of the ELBO later. Owing to the sparse decomposition of the covariance matrix, we can cheaply draw marginal samples for an f i from q(f ) with a linear transformation of white noise. The sampling scheme is shown in (2) and pictured in the right part of Figure 1 .

EQUATION

The constructed distribution q(f ) well approximates the strong correlations in the prior. From (2), f i and f j correlate in q(f ) by sharing noise entries in n(i) ∩ n(j) when the intersection is not empty. In this case, either f i neighbors f j , or f j neighbors f i , or f i , f j share common neighbors. When the neighbor sets are large enough, most strong correlations will be approximated by some non-zero entries in V.

3.2. OPTIMIZATION OF THE ELBO:

We optimize the ELBO in (1) to find a good q(f ) to approximate the GP posterior. To apply stochastic optimization, we will decompose the three terms in the ELBO. We mainly consider the decomposition of L cross and L ent , as the decomposition of L ell is easy.

We first decompose the cross entropy L cross . By convention, the GP prior has a zero mean. Though there is a closed-form calculation of L cross with both q(f ) and p(f ) being multivariate Gaussian, it involves expensive calculations of det(Σ) and Σ -1 . Previous works approximate the prior with Vecchia's method for easy decomposition and good approximation (Vecchia, 1988; Stein et al., 2004; Datta et al., 2016; Liu & Liu, 2019; Finley et al., 2019) . The idea is to build a directed graphical model and approximate p(f )

≈ N i=1 p(f i |f α(i)

) with α(i) being a small parent set of i. In the original work, Vecchia (1988) first set an order to data points and then choose α(i) as the K nearest parents of i. But it is not easy to guarantee a good ordering of data points (Banerjee et al., 2014; Guinness, 2018 ). Here we consider all possible orderings and take the average of approximations to address the data ordering concern.

We estimate L cross as follows. First, we randomly sample a parent set n (i) ⊂ n(i) with i / ∈ n (i) for each i. Then we approximate the log-prior by log p(f )

≈ N i=1 log p(f i |f n (i) ), with the conditional distribution p(f i |f n (i) ) derived from the joint Gaussian p(f i , f n (i)

) in the prior. Then L cross is estimated by a random batch of terms. The complete calculation is given as

EQUATION

)

Here S is a random batch of data points.

Now we justify that this is an average over all data orderings. Suppose there is a data order π(•), such that we can define a directed graphical model over p(f ) by assigning every i a parent set n π (i) = {j : j ∈ n(i), π(j) < π(i)}. Denote Π as all permutations of N data points, with each permutation inducing a graphical model. The average of the log densities of all graphical models can be collapsed to the average computed from local neighborhoods. Denote Π n(i) as permutations of indices in the set n(i), then we have

EQUATION

)

Here we only need to consider permutations of data points within n(i) for each i. Then we obtain (3) by estimating the inner summation by a single random permutation of n(i) and the outer summation by a random batch S.

We then decompose the entropy L ent . The entropy of q(f ) requires the expensive computation of det(V). To circumvent this difficulty, we find a decomposable lower bound of the entropy by using an auxiliary distribution (Ranganath et al., 2016; Louizos & Welling, 2017) . Note that we always prefer a lower bound of the objective in this maximization problem. With an arbitrary distribution r( |f ), a lower bound of L ent is

EQUATION

)

The bound is tight when r( |f ) matches q( |f ). In this work, we try to let r( |f ) match q( |f ).

Particularly, we set r( |f

) = i q( i |f n(i) )

, where the conditional q( i |f n(i) ) is derived from the joint Gaussian distribution q( i , f n(i) ). Then all terms in the lower bound in (5) are Gaussian loglikelihoods and are decomposable over data points. We can then reach the estimation of the entropy lower bound with a batch of data points.

EQUATION

We finally decompose the likelihood L ell . The likelihood term log p(y|f ) naturally decomposes because y i -s are conditionally independent given f i -s.

EQUATION

Here for each term i in the summation, the expectation is estimated by a Monte Carlo sample fi from q(f i ). The gradients of variational parameters are propagated through fi via reparameterization (Kingma & Welling, 2013) .

Finally, the ELBO has a decomposable approximation Lell + Lcross + Lent to enable efficient stochastic optimization. From the derivations above, we see the objective can be decomposed by data points. The computation for a data point only involves itself and its K nearest neighbors. Therefore, each stochastic gradient calculation takes time only O(K 3 ). There are N (K + 1) parameters in µ and R to optimize, so the optimization takes at least O(N ) time. We further reduce the number of parameters by amortizing the cost through a shared inference model, taking advantage of the fact that the inference for each data point i only needs its K nearest neighbors.

3.3. AMORTIZED INFERENCE:

Following AIGP, we also apply amortized inference to GP inference. Particularly, we train an inference network to identify variational parameters (µ i and R i,n(i) ) for each data point i. Since node correlations at a neighborhood can be easily treated as a weighted graph, we use Graph Convolutional Networks (GCNs) (Kipf & Welling, 2017) as our inference network.

A GCN takes an adjacency matrix A ∈ R k×k of graph and the node features H (0) ∈ R k×d0 as the input and then makes predictions for all graph nodes. Let Ā be the normalized adjacency matrix, Ā = D -1 2 AD -1 2 , with D being the diagonal degree matrix. A GCN layer with the input H ( -1) is defined by H ( ) = g (H ( -1) , A) := σ ĀH ( -1) W ( ) . Here W (l) ∈ R d -1 ×d is the weight matrix of the layer . σ(•) is the activation function. An L-layer GCN computes its output by H = gcn(H 0 , A) := g L (. . . g 1 (H 0 , A) . . . , A). We use two GCNs for the inference task, gcn 1 for the calculation of µ i and gcn 2 for R i,n(i) :

EQUATION

Here we use Σ n(i),n(i) as the adjacency matrix and stack the observation y n(i) and the one-hot vector e i as the input feature. The vector e i indicates the element i for which the inference is running for. We choose the activation σ(•) to be ReLU for intermediate layers and identity for the last layer. The last layer of each GCN has size 1 to output a K × 1 vector. a is an averaging vector with all K elements as 1 K . The dashed box in Figure 1 shows the amortization. LAIN defines an inference subtask on a data point and its nearest neighbors, while AIGP defines a subtask on a data point and its parents. Due to this difference, LAIN has two advantages. First, the inference network of LAIN uses the observations from all the K nearest neighbors, while the inference network of AIGP uses observations from parents only but not children. Second, inference subtasks of LAIN are generated with the same mechanism because the nearest-neighbor relationship is homogeneous across all data points. However, the parent-child relationship in AIGP depends on the ordering of data points (e.g. the first one in the order does not have parents). As a learning model, the inference network prefers subtasks from the same "distribution". The computational cost of GCN is O(K 2 ) by treating the network size as constant. The complexity of one gradient calculation is O(K 3 ). The optimization procedure converges fast since it only optimizes a constant number of variational parameters. In practice, we often observe that the optimization procedure converges in less than one epoch, which is not possible for methods without amortization. Finding nearest neighbors is the only step with running time bounds to the data size, but it only needs one run and is often fast on medium to large data sizes. If the data has a very large size, we can use k-d trees for low-dimensional data and approximate algorithms (Arya et al., 1998; Datar et al., 2004) for high dimensional data.

3.4. PREDICTION:

For a new data point x with its K nearest neighbors n( ) in the prior, the predictive distribution is Here q(f |x , X n( ) , y n(

EQUATION

)

EQUATION

) with b = Σ ,n( ) Σ -1 n( ),n( ) .

F is a set of Monte Carlo samples from q(f |x , X n( ) , y n( ) ) . The Monte Carlo estimation is accurate since the integral is one-dimensional.

4. EXPERIMENT:

We compare our method with five state-of-the-art methods: SVGP (Hensman et al., 2015) , SAVIGP (Dezfouli & Bonilla, 2015) , DGP (Cheng & Boots, 2017) , VFF (Hensman et al., 2017) , and AIGP (Liu & Liu, 2019) . The first three methods are based on inducing points, VFF uses inter-domain inducing points, and AIGP uses local neighbors. Through all experiments, we use RBF as the default kernel, except for VFF we use Ma tern-3 2 kernel (the code does not provide RBF kernel). We use the implementation of SVGP from GPFlow (Matthews et al., 2017) , the implementation of DGP from Faust (2018) , and implementations of all other algorithms from their authors.

For SVGP, SAVIGP, and VFF, we vary the number of inducing points, M ∈ {200, 1000, 2000}, to check their performances. DGP has separate inducing points for mean approximation and those for variance approximation. We use 256 inducing points for variance approximation and vary the number of inducing points for mean approximation from 200 to 2000. We vary the number of neighbors, K ∈ {10, 20, 40}, for AIGP and LAIN. GCNs used in these two methods have three hidden layers with dimensions [20, 10, 1]. We randomly split each dataset into training (75%) and testing (25%) and report both the predictive performance on the test set and the inference running time. To save the space, we report results from two settings for each competing method: one setting is M = 200 or K = 10, with which all methods have their fastest speed (marked by ), and another setting giving the best predictive performance (marked by ).

4.1. A TOY EXAMPLE:

In this section, we test different methods on a one-dimensional toy example studied in (Snelson & Ghahramani, 2006) . The dataset contains 200 data points, shown as black dots in figure 2 . We assume Gaussian likelihood in this experiment and run exact inference as the baseline. A smaller GCN (hidden dimensions [10, 5, 1] ) is used in this task.

The predictive mean and variance from the exact inference and LAIN with K = 10 are shown in the first two plots of Figure 2 . The result of LAIN is very similar to that of exact inference, except that the mean curve of LAIN is less smooth, which does not really hurt the predictive performance.

We test different methods with very small M and K and observe how they behave. We are likely to face this situation when we work on large datasets in high-dimensional spaces. The last three plots of Figure 2 exhibit predictive distributions of SVGP with M = 2 inducing points, AIGP with K = 2 parents, and LAIN with K = 2 nearest neighbors. When there are not enough inducing points, SVGP over-smooths the prediction and performs poorly for a good fraction of data points. AIGP does not have a good predictive mean either, because under a random ordering the directed graph constructed by AIGP cannot well capture neighboring relations. The predictive mean of LAIN does not deviate far from the ground-truth in the area with training instances, though the curve is rugged due to local variations. In this experiment, we estimate the spatial abundance of a bird species (Savannah Sparrow) using eBird dataset (Munson et al., 2015) . The dataset has 14,393 observations, each of which is a reported bird count at a GPS location. We model the observed counts with GPS locations as the input. We set the likelihood to be a Poisson distribution, with rate given by λ i = exp(f i ).

We compare different inference methods in terms of Negative predictive Log-Likelihood (NLL, the smaller the better predictive performance). Table 1 shows the results. We can see that LAIN achieves the best predictive performance at K = 40. Methods based on inducing points generally perform worse. In this dataset, observations have strong correlations in local areas, but inducing points are not efficient to capture the posterior at such a fine scale. In terms of running speed, LAIN is comparable to AIGP and DGP but faster than other methods. In our experiment, we have also tried to increase inducing points for DGP, but it does not improve its performance.

In this experiment, we also investigate whether inference networks work correctly. We run LAIN without inference networks and optimize µ and L for the variational distribution directly. Then we compare LAIN models with and without inference networks by checking their optimization procedure. In this task, we fix hyperparameters, so the two methods solve a pure inference problem. Figure 3 is the trace plot of the negative ELBO versus training epochs. The figure shows the ELBO of the two LAIN models eventually converge to very similar values, though the ELBO without inference networks is slightly better after 50 epochs (likely due to the amortization gap). LAIN with inference networks significantly reduces the number of optimization epochs -the inference networks are well trained after only 0.01 epoch (about 100 iterations). In summary, the result indicates that inference networks can effectively identify the variational parameters using local information.

4.3. PRECIPITATION LEVEL ESTIMATION:

In this task, we evaluate LAIN on a rainfall dataset. We process the precipitation dataset (Climate Data Online) and obtain the average precipitation level in May at 8,832 stations that are spatially distributed in the US. The GP inputs are GPS locations of these stations, and the observations are the average precipitation levels. We use the log-normal distribution as the likelihood, with its mean as function value f from GP and variance as a hyperparameter learned from the data.

Table 2 summaries the experimental results. LAIN has better predictive performance, and its running speed is comparable to or faster than other methods.

We also analyze the goodness of our prior approximation since we can compute the exact L cross on this dataset. We compute L cross with the optimized q(f ) distribution as well as Lcross . The true value L cross and the approximation Lcross are: 5,465 versus 5,396 when K = 10, 9,905 versus 8,490 when K = 20, and 10,086 versus 9,009 when K = 40. This result indicates that the approximation Lcross is relatively accurate. Furthermore, Lcross tends to be smaller than the true value and can be considered as a lower bound in such cases.

4.4. HAND-WRITTEN DIGIT CLASSIFICATION:

In this experiment, we explore a high-dimensional inference problem, GP classification of MNIST digits (LeCun & Cortes, 2010) . We consider a binary classification on handwritten images of 5 and 8.

To make performance values of different methods more differentiable, we randomly choose a subset of size 7,858 from the original dataset. Pixel values are normalized to [0,1] in the preprocessing step. In the results, we also report the accuracy obtained by different methods.

The results are shown in Table 3 . We see that LAIN performs the best in terms of classification accuracy. Its predictive NLL and running speed also overperform competing methods, though not very significant. We also observe that AIGP makes less confident predictions than other methods, which accounts for its worse predictive NLL but high accuracy. We do not report results from VFF due to memory issues.

We also examine KNN in this experiment. From the results, we notice that a small number of neighbors are often sufficient for KNN and LAIN models to perform well. By checking the running time of KNN, we also see that the time of finding nearest neighbors is only a small fraction of the total inference time on this dataset. There are slight differences regarding the test accuracy between KNN and LAIN, presumably due to different weighting schemes: LAIN weights different nearest neighbors according to their correlations, while KNN treats all nearest neighbors uniformly.

5. CONCLUSION:

In this work, we propose a novel approach for GP inference. We construct a variational distribution that has a sparse decomposition on its covariance matrix. With this distribution, function value at a data point is inferred from its nearest neighbors, encouraging the inference efficiently focuses on approximating strong correlations posed by the prior. The proposed variational distribution is expressive to approximate the GP posterior and also provides a decent structure for efficient ELBO optimization. We further decompose the ELBO into homogeneous subtasks and therefore enable stochastic optimization. Finally, we devise inference networks to perform these subtasks and significantly reduce the number of variational parameters. Our proposed method performs well in terms of predictive performance and running speed on a series of benchmark tasks.

:

We use the term "nearest neighbors" for easy reference, but we actually consider K most correlated neighbors in the prior. When the kernel is based on some distance metrics, then they are K nearest neighbors.

Figures and Tables:

Figure1: The structure of the variational distribution. The left box shows the amortization, which fits µi-s and Rij-s from their related prior kernel and observations. The right part shows the generation process of fi-s.

Figure 2: The first two plots compare predictive distributions of full GP and LAIN with K = 10. The right three plots show how SVGP, AIGP, and LAIN perform with a very small number of inducing points/neighbors. Data points in blue circles are not well fitted.

Comparison on the eBird dataset.

Table content: <table><tr><td colspan="2">Method Config</td><td>Pred NLL Time</td></tr><tr><td>SVGP</td><td>M=200 M=1000</td><td>1.90±.03 107s 1.88±.02 4.5ks</td></tr><tr><td>SAVIGP</td><td>M=200 M=2000</td><td>2.04±.03 167s 1.99±.03 50ks</td></tr><tr><td>VFF</td><td>M=200 M=2000</td><td>1.91±.02 1.3ks 1.91±.02 13ks</td></tr><tr><td>DGP</td><td>M=200 M=2000</td><td>1.82±.02 96s 1.80±.02 213s</td></tr><tr><td>AIGP</td><td>K=10 K=20</td><td>1.79±.05 45s 1.71±.05 125s</td></tr><tr><td/><td>K=10</td><td>1.69±.03 55s</td></tr><tr><td>LAIN</td><td>K=20</td><td>1.65±.03 384s</td></tr><tr><td/><td>K=40</td><td>1.60±.02 1.3ks</td></tr></table>

Comparison on the precipitation dataset.

Table content: <table><tr><td colspan="2">Method Config</td><td>Pred NLL Time</td></tr><tr><td>SVGP</td><td>M=200 M=2000</td><td>1.57±.03 2.5ks 1.28±.03 42ks</td></tr><tr><td>SAVIGP</td><td>M=200 M=2000</td><td>1.70±.02 2.8ks 1.58±.02 50ks</td></tr><tr><td>VFF</td><td>M=200 M=2000</td><td>1.54±.03 9.1ks 1.53±.03 32ks</td></tr><tr><td>DGP</td><td>M=200 M=2000</td><td>1.07±.05 402s 1.00±.05 889s</td></tr><tr><td>AIGP</td><td>K=10 K=10</td><td>0.96±.03 155s 0.96±.03 155s</td></tr><tr><td/><td>K=10</td><td>0.74±.05 129s</td></tr><tr><td>LAIN</td><td>K=20</td><td>0.72±.05 903s</td></tr><tr><td/><td>K=40</td><td>0.69±.04 2.3ks</td></tr></table>

Comparison on the MNIST dataset.

Table content: <table><tr><td>Method</td><td>Config</td><td>Pred NLL</td><td colspan="2">Accuracy Time</td></tr><tr><td>SVGP</td><td>M=200 M=1000</td><td>0.053±.004 0.051±.004</td><td>98.4 98.5</td><td>623s 23ks</td></tr><tr><td>SAVIGP</td><td>M=200 M=200</td><td>0.339±.008 0.339±.008</td><td>51.7 51.7</td><td>6.5ks 6.5ks</td></tr><tr><td>DGP</td><td>M=200 M=2000</td><td>0.059±.005 0.052±.005</td><td>98.1 98.3</td><td>292s 2.1ks</td></tr><tr><td>AIGP</td><td>K=10 K=40</td><td>0.293±.002 0.215±.003</td><td>98.0 98.2</td><td>3.9ks 24ks</td></tr><tr><td/><td>K=10</td><td>0.053±.003</td><td>98.9</td><td>128s</td></tr><tr><td>LAIN</td><td>K=20</td><td>0.050±.003</td><td>99.0</td><td>632s</td></tr><tr><td/><td>K=40</td><td>0.051±.003</td><td>99.1</td><td>2.9ks</td></tr><tr><td/><td>K=9</td><td>N.A.</td><td>98.6</td><td>24s</td></tr><tr><td>KNN</td><td>K=19</td><td>N.A.</td><td>98.3</td><td>26s</td></tr><tr><td/><td>K=39</td><td>N.A.</td><td>97.6</td><td>28s</td></tr></table>

Peer review (iclr2020):
PAPER SUMMARY: This paper proposes a fast inference method for Gaussian processes (GPs) that imposes a sparse decomposition on the VI approximation of the posterior GP (for computational efficiency) using the KNN set of each data point. This is further coupled with armortized inference for better scalability. NOVELTY & SIGNIFICANCE: This paper adopts a different approach of characterizing the VI approximation of a GP posterior than original VI approximation that was developed in Titsias (2009): Instead of characterizing the surrogate q(f_I) of p(f_I | Y) for a small collection of inducing inputs, the proposed method characterize q(f) directly where q(f) = int_f_I q(f_I) p(f | f_I)df_I. This is, however, a somewhat strange direction which, to me, seems to raise extra issues that could have been avoided if one follows the conventional VI approximation: (1) As the posterior surrogate is now directly over f instead of f_I, the number of variational parameters is now proportional to the data size which requires several (redundant) extra approximations including armortized inference & the lower-bound on the entropy term that admits a sparse decomposition. (2) This at least creates the armortized and entropy gaps that decrease the expressiveness of the original VI. While I understand that this is in exchange for the ability to encode local information (via KNN) within the surrogate posterior, it is not clear to me why do we need to incur all these computational issues to incorporate such local information. (3) For example, instead of forcing such local information in the posterior surrogate q(f), we could alternatively let it be reflected in the test conditional p(f_* | f_I, Y_n(*)) such that the test output depends on both the inducing output and a local partition of data (e.g., via K-mean), which has been previously explored in [*] and later incorporated in the conventional VI paradigm of Titsias (2009) without incurring extra intractability [**]. (4) This maintains the dense correlation between data points within the same neighborhood while allowing the VI surrogate to be more concisely specified and independent of the no. of training data points. Furthermore, it also helps avoid the data-bound overhead of computing a KNN per test point. [*] Local and global sparse Gaussian process approximations (AISTAT-07) [**] A distributed variational inference framework for unifying parallel sparse gaussian process regression models (ICML-16) To summarize, the practical significance of placing such a VI approximation directly on q(f) to encode such (indirectional) local information is, given the above, questionable to me. Please note that I am not disputing the potential use of this VI form here, which could have been the only way to encode a different (directional) type of information. For encoding KNN information, however, it only seems to create more troubles than it solves. Minor point: The above references, especially [*], should have been cited. TECHNICAL SOUNDNESS: [A] Optimization of the ELBO: (1) The ordering of data (i.e., the directional information) was mentioned repeatedly in the paper but its importance to the fast approximation was neither explained nor discussed. (2) The decomposition form of Eq. (6) also raises a question: How do we know that the term inside the log is positive? There seems to be missing information on the constraint of R. [B] Amortized Inference: (1) The choice of the GCN seems arbitrary here. I am in fact not sure why GCN is necessary for the inference network & furthermore, GCN also brings to the table another heuristic choice of A. (2) How do we set the adjacency graph A? (3) How do we know what is the right complexity for the GCN? [C] Complexity: The complexity analysis is too informal and lacking fine-grained information. Please include a detailed complexity analysis of the training and inference cost in terms of the input dimension, the no. of data points, the size of the neighborhood and the batch size. It is also necessary to factor in the KNN overhead (e.g., the cost of building the K-D tree for low-dimensional embedding of data & the approximation cost of projecting that information to high-dimensional data) EXPERIMENT: The experiment results only show marginal improvement over the baselines, and the size of the dataset for regression is too small. If I read correctly, both have fewer than 20000 data points. SVGP in particular has been tested on a much larger datasets (AIRLINE, UK Housing) comprising millions of data points -- comparison on such dataset should have been reported. Note that the largest dataset used to evaluate the efficiency of fast approximation of GP is on the scale of 6M data points [****]. On that note, eBird and precipitation should not even be considered mid-sized. To demonstrate the efficiency of local information encoding, comparison with [*] should be reported. There is another class of inducing-point methods that use expectation propagation that should have been discussed and/or compared with: [***] A Unifying Framework for Gaussian Process Pseudo-Point Approximations using Power Expectation Propagation (JMLR-18) [****] Distributed Gaussian Processes (ICML-15) CLARITY: The paper is clearly written. REVIEW SUMMARY: This paper adopts a VI approximation that deviates from the conventional form of (Titsias, 2009) to encode the KNN information, which causes extra computational issues (that incurs extra approximations). I find this deviation redundant seeing that the same information could have also been accounted for using the old VI form, which is a lot more computational efficient. I also find the experiment lacking as comparison with fast approximation method such as [*] that incorporate local information is not included. There are also a few technical ambiguities that need to be clarified. ------ Post-Rebuttal Update ------ Thank you for the rebuttal & I have read it in detail. However, it still does not address my concern, which I re-summarize here: I do not dispute the beneficial of exploiting neighborhood information but my point is we could still leverage the same amount of neighborhood information without going through the trouble of incurring extra steps of approximation due to approximating q(f) instead of q(f_I) -- I think I have elaborated this in points (1) - (4) in my original review -- which also creates the amortized gap. I am also not sure what model reuse means (in this context) and whether it is relevant since it appears somewhat noncentric to the objective of this paper. Also, apparently, the experiment is still lacking since to me, comparing with [*] is important in substantiating the contribution claim of this paper.