Scaling Gaussian Processes with Derivative Information Using Variational Inference
Year: 2022

Abstract:
Gaussian processes with derivative information are useful in many settings where derivative information is available, including numerous Bayesian optimization and regression tasks that arise in the natural sciences. Incorporating derivative observations, however, comes with a dominating O(ND) computational cost when training on N points in D input dimensions. This is intractable for even moderately sized problems. While recent work has addressed this intractability in the low-D setting, the high-N , high-D setting is still unexplored and of great value, particularly as machine learning problems increasingly become high dimensional. In this paper, we introduce methods to achieve fully scalable Gaussian process regression with derivatives using variational inference. Analogous to the use of inducing values to sparsify the labels of a training set, we introduce the concept of inducing directional derivatives to sparsify the partial derivative information of a training set. This enables us to construct a variational posterior that incorporates derivative information but whose size depends neither on the full dataset size N nor the full dimensionality D. We demonstrate the full scalability of our approach on a variety of tasks, ranging from a high dimensional stellarator fusion regression task to training graph convolutional neural networks on Pubmed using Bayesian optimization. Surprisingly, we find that our approach can improve regression performance even in settings where only label data is available.

heading: 1 Introduction
text: Gaussian processes (GPs) are a popular tool for probabilistic machine learning, widely used in scenarios where uncertainty quantification for regression is necessary [27, 38, 14]. When used for Bayesian optimization (BO) [18, 30], or in some regression settings found in the physical sciences like estimation of arterial wall stiffness, derivative information may be available [37, 34]. In these settings, we have not only noisy function values y = f(x) + but also noisy gradients ∇y = ∇xf(x) + at some set of training points X ∈ RN×D. On paper, GPs are ideal models in these settings, because they allow for training on both labels y and gradients∇y in closed form. Though analytically convenient, Gaussian process inference with derivative information scales poorly: computing the marginal log likelihood and predictive distribution for an exact GP in this setting requires O(N3D3) time and O(N2D2) memory. Recent work has addressed this scalability in certain settings, e.g. for many training points in a low-dimensional space [5] or for few training points in a high-dimensional space [3]. Despite these advances, training and making predictions for a GP with derivatives remains prohibitively expensive in regimes where both N and D are on the order of hundreds or even thousands.
We introduce a novel method to scale Gaussian processes with derivative information using stochastic variational approximations. We show that the expected log likelihood term of the Evidence Lower
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Bound (ELBO) decomposes as a sum over both training labels and individual partial derivatives. This lets us use stochastic gradient descent with minibatches comprised of arbitrary subsets of both label and derivative information. Just as variational GPs with inducing points replace the training label information with a set of learned inducing values, we show how to sparsify the derivative information with a set of inducing directional derivatives. The resulting algorithm requires only O(M3p3) time per iteration of training, where M N and p D. We demonstrate the quality of our approximate model by comparing to both exact GPs with derivative information and DSKI from [5] on a variety of synthetic functions and a surface reconstruction task considered by [5]. We then demonstrate the full scalability of our model on a variety of tasks that are well beyond existing solutions, including training a graph convolutional neural network [16] on Pubmed [28] with Bayesian optimization and regression on a large scale Stellarator fusion dataset with derivatives. We then additionally show that, surprisingly, our variational Gaussian process model augmented with inducing directional derivatives can achieve performance improvements in the regression setting even when no derivative information is available in the training set.

heading: 2 Background
text: In this section we review the background on Gaussian processes (GP) (Section 2.1), Gaussian processes with derivative observations (Section 2.2), and variational inference inducing point methods for training scalable Gaussian processes (Section 2.3).
Derivative notation. Throughout this paper for compactness we abuse notation slightly and use ∂jyi to refer to the jth element of ∇yi. In this particular case, this would correspond to the partial derivative observation in dimension j for training example xi. We also use ∂vyi to refer to the directional derivative in the direction v, i.e. ∇y>i v.

heading: 2.1 Gaussian processes
text: A Gaussian process (GP) is a distribution over functions f ∼ GP(µ(x), k(x,x′)) specified by mean and covariance function µ, k [23]. Given data points X = {x1, ...,xN} and function observations f = {f(x1), ..., f(xN )}, placing a GP prior assumes the data is normally distributed with f ∼ N (µX ,KXX) where µX is the vector of mean values atX andKXX ∈ RN×N is a covariance matrix. Conditioning on noisy observations y = f + where ∼ N (0, σ2I) induces a posterior distribution p(f∗|y) over the value of f at points x∗, which is Gaussian with mean µ(x∗) − Kx∗X(KXX + σ2I)−1(f − µX) and covariance k(x∗,x∗) − Kx∗X(KXX + σ2I)−1KXx∗ . Thus, standard GP inference takes O(N3) time. Hyperparameters such as σ, θ are generally estimated by Maximum Likelihood. The log marginal likelihood
L(X, θ, σ|y) = −1 2 (y−µX)T (KXX +σ2I)−1(y−µX)− 1 2 log |KXX +σ2I| − n 2 log(2π) (1)
can be optimized with methods like BFGS [19] at a complexity of O(N3) flops per iteration.

heading: 2.2 Gaussian processes with derivatives
text: GPs can leverage derivative information to enhance their predictive capabilities. Notably, as differentiation is a linear operator, the derivative of a GP is a GP [20]. Derivative observations can then be naturally included in a GP by defining a multi-output GP over the tuple of function observations and partial derivative observations (y,∇y) [24]. The GP has mean and covariance functions
µ∇(x) = [ µ(x) ∇xµ(x) ] , k∇(x,x′) = [ k(x,x′)
( ∇x′k(x,x′) )T ∇xk(x,x′) ∇2k(x,x′) ] . (2)
While including partial derivative observations can enhance prediction of f , and vice versa, a price is paid in the computational cost, as training and inference of GPs with derivatives scale as O(N3D3). This scalability issue has been addressed in the low D setting, and is discussed in section 3.

heading: 2.3 Stochastic Variational Gaussian Processes
text: Inducing point methods [29, 22, 33, 10] achieve scalability by introducing a set of inducing points: an “artificial data set” of points Z = [zj ]Mj=1 with associated inducing values, u = [uj ] M j=1. Stochastic
Variational Gaussian Processes (SVGP) [9] augment the GP prior p(f | X)→ p(f | u,X,Z)p(u | Z) and then learn a variational posterior q(u) = N (m,S). Inference for an observation y∗ at x∗ takes time O(M3):
q(y∗) = N (y∗;Kx∗ZK−1ZZm, σf (x ∗)2 + σ2) (3)
where σf (x)2 = Kxx − KxZK−1ZZKZx + KxZK −1 ZZSK −1 ZZKZx is the data-dependent variance. Using Jensen’s inequality and the variational ELBO [10, 11], SVGP develops a loss that is separable in the training data and amenable to stochastic gradient descent (SGD) [25], as the Kullback-Leibler (KL) divergence regularization only depends on u
ELBOSVGP = N∑ i=1 { logN (yi|µf (xi), σ2)− σf (xi) 2 2σ2 } − KL [q(u)||p(u)] . (4)
In equation (4), µf (xi), σf (xi)2 are the predicted mean and variance, respectively. The ELBO is maximized over the variational parameters m,S and the GP hyperparameters θ. Training with SGD on mini-batches of B data points brings the time per iteration to O(BM2 +M3).
While SVGP scales well, its predictive variances are often dominated by the likelihood noise [13]. Modeling derivatives necessarily involves heteroscedastic noise, or at least different noise for the function values and gradients, which may make SVGP with a Gaussian likelihood ill-suited to the task. The Parametric Gaussian Process Regressor (PPGPR) achieves heteroscedastic modeling by using the latent function variances without modifying the likelihood by symmetrizing the dependence of the loss on the data-dependent variance σf (xi)2 term
ELBOPPGPR = N∑ i=1 logN (yi|µf (xi), σ2 + σf (xi)2)− KL[q(u)||p(u)]. (5)
In Section 5, we evaluate our approach as an extension to both SVGP and PPGPR, and find that PPGPR gives significant performance gains.

heading: 3 Related Work
text: DSKI and DSKIP [5], derivative extensions of SKI [36] and SKIP [7], are among the first methods to address scaling Gaussian processes with derivative information in a low dimensional setting. DSKI and DSKIP approximate derivative kernels by differentiating interpolation kernels ∇k(x,x′) ≈∑
i∇wi(x)k(xi,x′) where wi(x) are interpolation weights used in SKI. Like SKI, DSKI suffers from the curse of dimensionality, and matrix-vector products cost O(ND6D + M logM) time. DSKIP improves the dependence on D, but still costs O(D2(N +M logM + r3N logD)) to form the approximate kernel matrices, where r N is the effective rank of the approximation. Thus while these methods exhibit high model fidelity, they are limited to low dimensional settings.
Recently [3] introduced an exact method for training GPs with derivatives in time O(N2D+ (N2)3), which improves on the naive O(N3D3) when N < D. This method is not applicable as N grows moderately large, while our paper chiefly focuses on the high-N and high-D setting.
Bayesian optimization with derivatives was considered in [37]. Here, the authors consider conditioning on directional derivatives to achieve some level of scalability, but the dataset sizes considered never exceed N of around 200 or D of around 8. Distinct from their consideration of directional derivative information, we will be equipping each inducing point in a sparse GP model with its own set of distinct directional derivatives, allowing the model to learn derivatives in many directions in regions of space where there are a large number of inducing points.

heading: 4 Methods
text: Our goal is to enable training and inference on data sets with large N and D when derivatives are available. We will present our method in three steps. First, we describe a naive adaptation of stochastic variational Gaussian processes to the setting with derivatives. Second, we argue that this adaptation again scales poorly in D. Finally, we show that using additional sparsity on the derivatives gives us scalability in both N and D.

heading: 4.1 Variational Gaussian processes with derivatives.
text: As described in Section 2.3, SVGP creates a dataset of inducing points Z = [zj ]Mj=1 with labels (or inducing values) u = [uj ]Mj=1. Assume we are given a dataset X = [xi] N i=1 with labels y = [yi] N i=1 and derivative observations ∇y = [∇yi]Ni=1. A natural extension of SVGP to this data is to augment the standard inducing dataset with inducing derivatives,∇u = [∇uj ]Mj=1, each of length D, so that each inducing point becomes a triple (zj , uj ,∇uj). This corresponds to a new augmented GP prior:
p(f ,∇f | X)→ p(f ,∇f | u,∇u,X,Z)p(u,∇u | Z). (6) Analogous to SVGP, we introduce a variational posterior:
q(u,∇u) = N ( m∇,S∇ ) = N ([ m ∇m ] , [ S ∇S ∇S> ∇2S ]) . (7)
Here, m and∇m are trainable parameters learned by maximizing the ELBO. We abuse notation and call the second portion of the vector∇m because these variational mean parameters correspond to the M ×D inducing derivative values. This also holds for the matrices∇S and ∇2S. With this augmented variational posterior, the ELBO becomes:
Eq(f ,∇f) [log p(y,∇y | f ,∇f)] − KL(q(u,∇u)||p(u,∇u)). (8) Assuming the typical iid Gaussian noise likelihood for regression and expanding the first term further:
Eq(f ,∇f) [log p(y,∇y | f ,∇f)] = N∑ i=1 Eq(fi) [log p(yi | fi)] + N∑ i=1 D∑ j=1 Eq(∂jfi) [log p(∂jyi | ∂jfi)] . (9)
Here, we have used linearity of expectation and the conditional independence between yi and ∂jfi given fi to show that the term of the ELBO that depends on training data decomposes as a sum over labels yi and partial derivatives ∂jyi. Thus, minibatches can contain an arbitrary subset of labels (xi, yi) and partial derivatives (xi, ∂jyi), and the minibatch size B remains independent of N and D.
The moments of q(f ,∇f) = ∫ p(f ,∇f | u,∇u)p(u,∇u) du d∇u are similar to those in SVGP, but the kernel matrices have been augmented with derivatives (i.e., using the kernel k∇(x,x′)):
µf ,∇f = K ∇ XZ ( K∇ZZ )−1 m∇, Σf ,∇f = K ∇ XX +K ∇ XZK ∇−1 ZZ (S ∇ −K∇ZZ) ( K∇ZZ )−1 K∇ZX . (10)
Here, K∇XX is a B ×B matrix that corresponds to a randomly sampled subset of label and partial derivative information. K∇XZ is B ×M(D + 1), and both S∇ and K∇ZZ are M(D + 1)×M(D + 1). Similarly, the KL divergence KL(q(u,∇u)||p(u,∇u)) involves multivariate Gaussians with covariance matrices of size M(D + 1)×M(D + 1). As a result, the running time complexity of an iteration of training under this framework is O(M3D3) which, grows rapidly with dimension.

heading: 4.2 Variational Gaussian processes with directional derivatives.
text: The procedure above is deceptively expensive despite the asymptotic complexity of a single iteration. Because a minibatch of size B contains an arbitrary subset of the N labels and ND partial derivatives rather than simply a subset of the N labels, each epoch in the above procedure must process roughly N+ND
B minibatches, rather than the usual N B . Additionally, because K ∇ ZZ is of size M(D + 1) × M(D + 1), the above procedure is also analogous to SVGP using M(D + 1) inducing points rather than using M . While minibatch training adapts readily to N(D + 1) training examples, it is rare to use significantly more than 1000 inducing points, which can require specialized numerical tools to make scale even to M = 10000 [21]. In practice, M(D + 1) would rapidly result in matrices K∇ZZ that make training infeasibly slow.
To make the matrix K∇ZZ not directly scale with the input dimensionality, we replace the inducing derivatives from equation (7) with inducing directional derivatives. Rather than the triplet (zi, ui,∇ui) with ∇ui having dimension D, each inducing point is now equipped with a set of p distinct directional derivatives (zi, ui, ∂Vi1ui, ..., ∂Vipui) in the directions vi1, ...,vip. We include the inducing directions V = [V1 · · ·VM ] ∈ RMp×D as trainable parameters.
GPs with Directional Derivatives. Similar to how we built the derivative kernel matrix in Section 2.2, we may define a multi-output GP over an unknown function and its directional derivatives. For a point zi and some direction vi and another point and direction zj and vj the directional-derivative covariance function is:
k∂vi∂vj (zi, zj) =
[ k(zi, zj) ∇zjk(zi, zj)>vj
v>i ∇zik(zi, zj) v>i ∇2zizjK(zi, zj)vj
] , (11)
which is of size 2× 2 rather than (D + 1)× (D + 1) as with k∇(·, ·).
GivenMp inducing directions V, p per each of theM inducing points, the relevant kernel matrices (1) between all pairs of inducing values and directional derivatives, KZZ , and (2) between all inducing values and training examples with full partial derivative observations, KXZ , are:
KZZ =
[ KZZ ∇ZKZZV
V >∇ZKZZ V >∇2ZZKZZV
] , KXZ = [ KXZ ∇ZKXZV ∇XKXZ ∇2ZKXZV ] , (12)
the first of which has shape M(p+ 1)×M(p+ 1). Constructing KZZ ,KXZ is inexpensive as we compute them directly from the directional derivative kernel (11), rather than computing the full gradient kernel k∇ and multiplying by the directions V which would incur a cost of O(M2D2).
Variational inference with this model is nearly identical to inference with full inducing gradients. We define a variational posterior, this time over the M(p+ 1) inducing values and directional derivatives:
q(u, ∂Vu) = N ( m,S ) (13)
where m ∈ RM(p+1) and S ∈ RM(p+1)×M(p+1). The model is trained by optimizing the variational ELBO
Eq(f ,∇f) [log p(y,∇y | f ,∇f)] − KL(q(u, ∂Vu)||p(u, ∂Vu)). (14)
over the variational parameters m,S, the inducing points Z, the inducing directions V, and the hyperparameters. Because the structure of the ELBO remains unchanged, the training labels and partial derivatives can again be subsampled to form minibatches of size B, yielding KXZ ∈ RB×M(p+1). Inference proceeds by computing q(f ,∇f) from (10) by replacing the kernel matrices K∇XZ and K∇ZZ with our directional derivative variants KXZ and KZZ .
Learning inducing directions The above algorithm requires the selection of a set V = [V1 · · ·VM ] of inducing directions. In the setting where all inducing points have a shared set of p global inducing directions, there is an optimal fixed choice for the inducing directions [2, 5]. However, for a variational GP with directional derivatives, sharing inducing directions does not improve scalability as the size of the kernel matrices KXZ ,KZZ is unchanged from the case where each inducing point has p distinct inducing directions. Thus we may improve model flexibility and performance by allowing each inducing point to have distinct directions at no additional computation complexity (see the supplementary materials for a comparison against shared inducing directions). In this case, a principled approach to setting the inducing directions is to include them as a set of mpd trainable parameters, and learn them when maximizing the ELBO (14). Learning inducing directions allows nearby inducing points to balance the directions from which they capture information, and encourages the model to capture the locally most informative directions. We adopt this approach in our experiments in section 5.
Derivative modeling with p D. A key feature of this framework is that it allows for the use of a different number p of directional derivatives per inducing point than the number of partial derivative observations per training point. Particularly for kernel matrices involving training examples with full partial derivative information, using p D directional derivatives keeps the matrix dimension small and independent of D. Nevertheless, allowing each inducing point to have its own set of learnable directions enables the model to learn many derivative directions where necessary in the input space by placing multiple inducing points with different directions nearby. A notable case is when each inducing point zi has the p = D canonical inducing directions Vi = I , through which we recover the full variational GP with derivatives as described in section 4.1.
Complexity. For a minibatch size B, when learning p directional derivatives per inducing point, the matrices KXZ and KZZ become B ×M(p+ 1) and M(p+ 1)×M(p+ 1) respectively. As a result, the time complexity of variational GP inference with directional derivatives is O(M3p3). When using p directions per inducing point, this is computationally equivalent to running SVGP with p+ 1 times as many inducing points. To counteract the additional matrix size, one may use the whitened formulation of variational inference [17] for GPs when computing equation (10) and use contour integral quadrature as in [21].

heading: 5 Experiments
text: In this section we evaluate the empirical performance of variational GPs with directional derivatives. In sections 5.1 and 5.2 we benchmark the performance of our method on small regression problems where we can compare to prior work such as exact GPs with derivatives and DSKI. In sections 5.3, 5.4, and 5.5 we compare to variational GPs without derivatives on high dimensional regression and Bayesian optimization (BO) tasks which are well beyond the scalability means of all prior work we are aware of. In section 5.6 we perform an ablation study to understand the effect of increasing the number of inducing directions. In section 5.7 we investigate the value of learning directional derivative information even when derivative observations are not available. All of our GP models use a constant prior and Gaussian kernel (or associated directional derivative kernel) and were accelerated through GPyTorch [8] on a single GPU. Code is available at https://github.com/ mishapadidar/GP-Derivatives-Variational-Inference.

heading: 5.1 Synthetic functions
text: In this section, we consider low-dimensional regression with derivatives on test functions including Branin (2D), SixHumpCamel (2D), Styblinksi-Tang (2D) and Hartmann (6D) from [32], a modified 20D Welch test function [1] (Welch-m)1, and a 5D sinusoid f(x) = sin(2π||x||2) (Sin-5). We compare variational GPs without derivatives (SVGP, PPGPR) to variational GPs with derivatives (GradSVGP, GradPPGPR), exact GPs with derivatives (GradGP), non-variational GPs with derivatives (DSKI), and variational GPs with p = 2 directional derivatives per inducing point (DSVGP2,DPPGPR2). Exact and variational GPs with full derivatives are only tractable in lowdimensional settings due to the scalability issues mentioned in sections 3 and 4.1. Therefore, to apply GradSVGP and GradPPGPR on the 20D Welch-m function, we first perform dimension reduction onto a low dimensional active subspace [2], similar to [5]. To show the limitation of GradSVGP and GradPPGPR, we modified the Welch function to have a low-quality low-dimensional active subspace.
In this low-dimensional setting, we find that variational GPs with directional derivatives, DSVGP2 and DPPGPR2, perform comparably to the methods that incorporate full derivatives (DKSI, GradSVGP, GradPPGPR, GradGP); see Table 1. In Figure 1 we compare the negative log likelihood of each method as the inducing matrix size grows on the Sin-5 and Hartmann test functions. We find that DSVGP2 and DPPGPR2 often outperform other methods due to their ability of incorporating derivative information while only modestly increasing the inducing matrix size.

heading: 5.2 Implicit Surface Reconstruction
text: In order to further validate the fidelity of our method’s derivative modeling, we consider the surface reconstruction task considered in [5]. We compare to DSKI with the goal of achieving comparable performance, as DSKI is nearly exact for this problem. In Figure 2, we reconstruct the Stanford Bunny by training DSVGP with p = 3 inducing directions for 1200 epochs and DSKI on 11606 noisy observations of 34818 locations and corresponding noise-free surface normals (gradients of the bunny level sets). DSVGP smoothly reconstructs the bunny and is comparable to DSKI.

heading: 5.3 Training Graph Convolutional Neural Networks with Bayesian Optimization
text: In this section, we demonstrate the full scalability of our approach by training the D = 4035 parameters of a two layer graph convolutional neural network (GCN) [16] on the node classification
1The Welch test function has intrinsically a 6D active space. We modified it to have a low-quality 6D active subspace and to show the limitation of GradSVGP and GradPPGPR.
task of the Pubmed citation dataset [28] using Bayesian optimization. The Bayesian optimization setting compounds the need for scalability, as the GP model must be retrained after each batch of data is acquired. For example, in the last 500 of 2500 optimization iterations with a batch size of 10, a GP must be fit 50 times to datasets with N(D + 1) ≈ 2500(4035 + 1) > 106 combined function and partial derivative labels. Any one of these datasets would be intractable to existing methods for training GPs with gradient observations.
For this experiment, we make no effort to modify the Bayesian optimization routine itself to account for the derivative information (e.g., as in [37]), as this would confound the performance improvements achieved by higher fidelity modelling by incorporating derivative information. Instead, we focus only on swapping out the underlying Gaussian process model. We consider TuRBO [6] as a base Bayesian optimization algorithm which we run with an exact GP, PPGPR, DPPGPR1, DPPGPR2 surrogate models. See supplementary material for p = 3 results. We additionally include traditional BO with the Lower Confidence Bound (LCB) [31] acquisition function, Adam [15] and random search. All algorithms were initialized with 400 random evaluations, and the TuRBO variants were run with a batch size of 20 and retrained over 150 steps. Figure 3 summarizes results averaged over 6 trials. We observe that TuRBO with DPPGPR significantly outperform traditional BO and other TuRBO variants. While all Bayesian optimization methods under-perform compared to Adam, we conjecture that this performance gap could be narrowed by incorporating the gradient information into the Bayesian optimization algorithm itself.

heading: 5.4 Stellarator Regression
text: In this experiment, we show the capacity of variational GPs with directional derivatives to extend GP regression with derivatives to massive datasets in a high dimensional settings. We perform regression on N = 500000 function and gradient observations gathered from a D = 45 dimensional optimization objective function through the FOCUS code [39]: a code for evaluating the quality of magnetic coils for a Stellarator, a magnetic confinement based fusion device for generating renewable energy [12]. The dataset is available upon request.
We compare variational GPs with directional derivatives using p = 1, 2 directions (DSVGP1, DPPGPR1, DSVGP2, DPPGPR2) to variational GPs without directional derivatives (SVGP, PPGPR). See supplementary material for p = 3 results. While D is too large to use variational or exact GPs with derivatives, we can apply variational GPs with derivatives to a projection of the data set onto a low-dimensional active subspace as in section 5.1. Variational GPs with derivatives trained on reduced datasets of dimension two and three performed poorly compared to all other methods tested. The results of this experiment are shown in Figure 4: variational GPs with directional derivatives significantly enhance regression performance. Even the inclusion of one directional derivative is enough the enhance the predictive capabilities of the regressor. The experiments were averaged over 5 trials, using an Adam optimizer with a multi-step learning rate scheduler and 1000 epochs.

heading: 5.5 Rover Trajectory Planning
text: The rover trajectory planning problem [6, 35] is a D = 200 dimensional optimization problem with gradients. This experiment validates the use of variational GPs with directional derivatives in Bayesian optimization by leveraging derivative information in a setting where no other methods can. We solve a variant of the rover problem in which the goal is to find an open-loop controller that minimizes the energy of guiding a rover through a series of waypoints in the xy-plane. The rover trajectory is integrated over 100 steps at which forces in x and y directions are applied to the rover,
making a total of D = 200 decision variables. We compare the performance of TuRBO, TuRBO with DPPGPR using p = 1, 2, traditional Bayesian optimization with the LCB acquisition function, and random search. All algorithms were initialized with a 100-point experimental design, and the TuRBO variants were run with a batch size of 5, and retrained over 300 steps. The results, averaged over 5 trials, are summarized in Figure 4. We observe that the TuRBO variants that leverage derivative information outperform the other algorithms almost immediately.

heading: 5.6 Ablation Study Over Number of Inducing Directions
text: Increasing the number of inducing directions p, should increase model flexibility and allow for the model to capture the full set of derivative information. However, as inducing directions are distinct to each inducing point, and are trainable parameters, each inducing directional derivative can encode information from multiple partial derivative directions. Thus variational GPs with directional derivatives may only need a small set of inducing directions to accurately encode the full set of derivative information. In table 2 we show the effect of increasing p on the performance of DPPGPR in the GCN Bayesian optimization (Section 5.3), stellarator regression (5.4), and rover Bayesian optimization (5.5) experiments. We find that using p = 1 is usually sufficient for capturing much
of the benefit of derivative information, and that increasing p can improve performance but at diminishing returns. See supplementary material for more p = 3 results on the GCN, Stellarator, and Rover tasks.

heading: 5.7 UCI Regression
text: Increasing the number of inducing points for Gaussian process models often results in diminishing returns on final model performance, with 500 ≤M ≤ 2000 often proving sufficient [9, 26, 13, 21]. This saturation is likely due in part to the ability of sparse Gaussian processes to represent variation
in the data, but may also be due to increasingly challenging optimization dynamics as more inducing points are added.
One hypothesis worth exploring is that, in some cases, it may be beneficial to augment a smaller set of inducing points with additional descriptive variables rather than to simply increase the number of inducing points. To that end, we test our method on a number of UCI benchmark regression datasets for which no derivative information is available. To accomplish this, a variational GP with directional derivatives is initialized with inducing directional derivatives as normal, but during training, minibatches of data only correspond to function observations. In other words, rows of the matrix KXZ never correspond to partial derivative observations because there are none, however columns of KXZ and rows/columns of KZZ still correspond to inducing directional derivatives as applicable. No changes to the inference procedure are needed.
We test on a number of UCI datasets [4]: Elevators (D=18 N=16599), Kin40k (D=8, N=40000), Energy (D=8, N=768), Protein (D=9, N=45730), Kegg-Directed (D=20, N=53414). We compare variational GPs without derivatives (SVGP, PPGPR) to variational GPs with p = 1 directional derivative per inducing point (DSVGP1, DPPGPR1). To ensure that methods have equivalent time complexity, we fix the inducing matrix size to 800. Thus, DSVGP1/DPPGPR1 use 400 inducing points, so that 400× (p+ 1) = 800, while SVGP/PPGPR use a full 800 inducing points. We use an 80-20 train-test split for all experiments, and train for 800 epochs. Interestingly, the results in Table 3 show that learning derivative information can improve prediction performance.
Elevators kin40k Energy Protein Kegg-directed RMSE NLL RMSE NLL RMSE NLL RMSE NLL RMSE NLL
SVGP 0.391 0.480 0.163 -0.302 0.231 -0.029 0.680 1.03 0.0870 -1.028 DSVGP1 0.379 0.449 0.134 -0.545 0.087 -1.05 0.654 0.998 0.085 -1.048
PPGPR 0.407 0.404 0.274 -0.874 0.252 -0.899 0.709 0.928 0.0915 -1.57 DPPGPR1 0.391 0.390 0.264 -1.105 0.251 -0.936 0.695 0.892 0.088 -1.643

heading: 6 Discussion

heading: 7 Acknowledgements
text: We acknowledge support from Simons Foundation Collaboration on Hidden Symmetries and Fusion Energy and the National Science Foundation NSF CCF-1934985, and NSF DMS-1645643.

Review summaries:
   The manuscripts proposes a method to accelerate GP inference in situations where derivative observations are available and the number of data points N and the data dimension D tends to be large. Extending previous work on stochastic variational approximations, the central idea of the paper is to use a number p of directional derivatives as variational parameters rather than full derivatives which greatly reduces the computational footprint.     Strengths    + Computational methods to tackle high D, high N GP inference problems with derivative information available are not fully scalable and hence a relevant research subject.    + The proposed approach is clear and builds on established techniques.    + The computational benefits demonstrated in an extensive set of experiments are substantial and convincing.      Weaknesses    - Experimental results are without error bars.

The authors propose a scalable inference method for GP with derivatives. Previous methods could only work for large data in low-dimensions or small data in high-dimensions. They solve the issues by using the inducing points and directional derivatives. They compare with other alternatives on multiple datasets and show that scalable GP with directional derivatives has the best performance.   

In a spirit similar to the SVGPs which sparsify the training data, the authors propose to augment the SVGP framework with derivative information by sparsifying the derivative data with *inducing directional derivatives*. By virtue of these modeling choices, the resulting algorithm scales favorably independent of data size as in SVGPs, and now additionaly also w.r.t. input dimension. Empirically, the authors show much better scalability in both high-N and high-D than existing work.