Scaling Gaussian Processes with Derivative Information Using Variational Inference
Authors: Misha A Padidar
Xinran Zhu
Leo Huang
Jacob R. Gardner
David Bindel
Venue: NeurIPS 2021 Poster (2021)
Keywords: Gaussian Process
Variational Inference
Bayesian Optimization

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(N^3D^3)$ 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 the 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 additionally find that our approach can improve regression performance even in settings where only label data is available.

Scaling Gaussian Processes with Derivative Information Using Variational Inference
Misha Padidar Xinran Zhu' Leo Huang',
Jacob R. Gardner?, David Bindel'
Cornell University, (map454, xz584, ah839, bindel) cornelleed University of Pennsylvania, seccbrg@seas..ueen.eede
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 dominating O(N D2) computational cost when training on N points in imput 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. 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 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 variety of tasks, ranging from high dimensional stellarator fusion regression task to training graph convolutional neural networks on Pubmed using Bayesian optimization. Surprisingly, find that approach can improve regression performance even in settings where only label data available.
1 Introduction
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 or some regression settings found in the physical sciences like estimation of arterial wall stiffness,, derivative information may be available [37 34). In these settings,
Bayesian optimization (BO) 18 30
we have not only noisy function values = fxx + e but also noisy gradients Vy = VxJ(x) + at some set of training points X E XD On paper, GPs are ideal models in these settings, because they allow for training on both labels y and gradients Vy 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 ( D3) time and O(N2 memory. Recent work has addressed this scalability in certain settings, e.g. for many training points low-dimensional space [5] or for few training points in high-dimenn nigh-dimensiona space 3. Despite these advances, training and making predictions for GP with derivatives remains prohibitively expensive in regimes where both and are on the order of hundreds or even thousands.
We introduce 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 (NcurIPS 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 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 set of inducing directional derivatives. The resulting algorithm requires only (M33p² time per iteration of training, where M < N and < D.
We demonstrate the quality of our approximate model by comparing to both exact GPs with derivative information and DSKI from [5] on variety of synthetic functions and surface reconstruction task considered by [S]. We then demonstrate the full scalability of our model on a variety of tasks that are well beyond existing solutions, including training graph convolutional neural network [16] on large scale Stellarator fusion dataset
Pubmed (28] with Bayesian optimization and regression 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 available in the training set.
2 Background
In this section we review the background on Gaussian processes (GP) (Section 2.1). Gaussian processes with derivative observations (Section 2.22, 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 o,y; refer to the jth element of Vy,. In this particular case, this would correspond to the partial derivative observation in dimension for training example x. We also use d,y; to refer to the directional derivative the direction v, i.e. Vy. v.
2.1 Gaussian processes
Gaussian process (GP) is a distribution over functions f GP()(()) x') specified by mean and covariance function μ,, k [23] Given data points X = {x..... XN) and function observations {f(x)).... (xx)) }} placing a GP prior assumes the data is normally distributed with ~ N(ux, Kxx) where μax is the vector of mean values X and Kxx NNN is covariance matrix. Conditioning on noisy observations = f + e where c <~ NN,, I²I) induces a posterior distribution p(("ly) over the value of f points x", which is Gaussian with mean p(x*)= Kx* (XxxXxx o- fx) and covariance k(x*, x") Kx- (Kxx + o²1)-1 Kxx- Thus, standard GP inference takes ((N3) time. Hyperparameters such as a, 0 are generally estimated by Maximum Likelihood. The log marginal likelihood
1
C(X, 8, fly) = (x==x) (Xxx px)
log [xx + 0²1 - log(2x) (1)
2
2
can be optimized with methods like BEGS 19 at a complexity of O(N3) flops per iteration.
2.2 Gaussian processes with derivatives
GPs can leverage derivative information to enhance their predictive capabilities. Notably, as differen- tiation is a linear operator, the derivative of a GP is a GP (20). Derivative observations can then be naturally included in GP by defining multi-output GP over the tuple of function observations and partial derivative observations (y, Vy) 24]. The GP has mean and covariance functions
((x) xx)(x)
A (x,x') = k(x, x')
(V <V= (²x(x,x))
# (x) =
(2)
V x')
While including partial derivative observations can enhance prediction of f. and vice versa, price is paid in the computational cost, as training and inference of GPs with derivatives scale as O D3). This scalability issue has been addressed in the low D setting, and discussed in section
2.3 Stochastic Variational Gaussian Processes
Inducing point methods (2
22 33 10 achieve scalability by introducing a set of inducing points an with associated inducing values, [u,]] Stochastic
"artificial data set" of points Z م

---

Variational Gaussian Processes (SVGP) 9 augment the GP prior p(f X) + p(f u. X. p) Z) and then learn a variational posterior (u) (m, S). Inference for an observation y" x" takes time O(M²):
dl') - t
where or(x)2 = Xxx KxzK2zK3 + K 1 KazKZz5a 7KK22KKKz 72 K 22x is the data-dependent variance. Using Jensen's inequality and the variational ELBO [10 I11, SVGP develops loss that separable in the training data and amenable to stochastic gradient descent (SGD) 25 as the Kullback-Leibler
(KI.) divergence regularization only depends
or oF(xx)² 2 20
ELBOsvGP
KL [a(u)]lz(())
(4)
(x,();,, (F)))))((()))))))))))))))))))( (x;)) are the predicted mean and variance, respectively. The ELBO is maximized over the variational parameters m. and the GP hyperparameters Training with SGD on mini-batches of B data points brings the time per iteration to O(BM² + M²).
In equation 4
While SVGP scales well, its predictive variances are often dominated by the likelihood noise [k 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 or oxxx;) term
ELBOPPGPR
log Nwn (r(x)), of + txx(x))? KL[g(u)]lP())]-
In Section 5 we evaluate our approach as extension to both SVGP and PPGPR, and find that PPGPR gives significant performance gains.
3 Related Work
DSKI and DSKIP (5), derivative extensions of SKI 36 and SKIP 77. are among the first methods to address scaling Gaussian processes with derivative information in low dimensional setting. DSKI and DSKIP approximate derivative kernels by differentiating interpolation kemels Vk(x, x') R E, (w,(x)x(xx.. x') where (x) are interpolation weights used in SKI. Like SKI, DSKI suffers from the curse of dimensionality, and matrix-vector veecorr products cost D6 DSKIP improves the dependence on D, but still costs O(D²(N + M log M + r N log D)) to form the approximate kernel matrices, where & N is the effective rank of the approximation. Thus while these methods exhibit high model fidelity, they are limited to low dimensional settings.
M log M) time.
Recently 3 introduced exact method for training GPs with derivatives in time O(N D (N3)3), which improves on the naive O( when N < D. This method is not applicable as N grows moderately large, while our paper chiefly focuses the high-N and high-D setting.
Bayesian optimization with derivatives was considered 337. Here, the authors consider condition- ing 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, will be equipping each inducing point in 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 large number of inducing points.
4 Methods
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 naive adaptation of stochastic variational Gaussian processes to the setting with derivatives. Second, we argue that this adaptation again scales poorly D. Finally, we show that using additional sparsity on the derivatives gives scalability in both N and D.
3

---

4.1 Variational Gaussian processes with derivatives.
As described in Section 2.3 SVGP creates a dataset of inducing points Z [2z with labels (or inducing values) u = [ j=1- Assume we are given a dataset X = [x:] |=] with labels y = [
and derivative observations Vy [V A natural extension of SVGP to this data is to augment the standard inducing dataset with inducing derivatives, Vu = [Vuy] each of length D, so that each inducing point becomes a triple (zj, Uj, Vuy). This corresponds to new augmented GP prior:
pff, Vf X) pff, Vf u, Vu, X, Z)p(u, Vu Z).
Analogous to SVGP, we introduce a variational posterior:
(u, Vu) N' (m (5) N
((m
( m
S
VS V2s
VS
Here, m and Vm are trainable parameters learned by maximizing the ELBO. We abuse notation and call the second portion of the vector Vm because these variational mean parameters correspond to the M x D inducing derivative values. This also holds for the matrices VS and V
With this augmented variational posterior, the ELBO becomes:
Eg vr) [lo Vy f, Vf) (L(g(u, Vu)llp(u Vu)).
Assuming the typical iid Gaussian noise likelihood for regression and expanding the first term further:
(e ply, Vy Vf)
N
D
Σ BA p(us fJ]+
E. [log 0zf0).
==
Here, we have used linearity of expectation and the conditional independence between yi and 0 given Si to show that the term of the ELBO that depends on training data decomposes as sum over labels ys and partial derivatives 0,yr- Thus, minibatches can contain an arbitrary subset of labels (xa, yt) and partial derivatives (x,, 0,y:), and the minibatch size B remains independent of N and D.
The moments of qf,, Vf) = f p(f, Vf u, Vu)p(u, Vu) du dVu u are similar to those in SVGP, but the kernel matrices have been augmented with derivatives (i.e., using the kernel A (x, x')):
NEVVE Kx (K)z) "mV, Σe,r- = XXX + K xz Z x K (10)
Here, K is a B x B matrix that corresponds to a randomly sampled subset of label and partial derivative information. k is B x M (D++ + 1), and both S and K are M(D+ 1) x M(D- - 1). Similarly, the KL divergence (L(q(u, Vujll [[ Vu)) involves multivariate Gaussians with covariance matrices of size M(D + 1) x MDD + 1). As a result, the running time complexity of an iteration of training under this framework is ((M3 D3) which, grows rapidly with dimension.
4.2 Variational Gaussian processes with directional derivatives.
The procedure above is deceptively expensive despite the asymptotic complexity of a single iteration. Because minibatch of size B contains an arbitrary subset of the N labels and N D partial derivatives rather than simply a subset of the N labels, each epoch in the above procedure must process roughly N+ND minibatches, rather than the usual Additionally, because KY KZ is of size M(D + 1) x MDD + 1), the above procedure is also analogous to SVGP using M M(D+ + 1) inducing points rather than using M. While minibatch training adapts readily N(D + 1) training examples, is rare to use significantly more than 1000 inducing points, which can require specialized numerical tools to make scale even M 10000 (2T). In practice, M( + 1) would rapidly result in matrices KY that make training infeasibly slow.
To make the matrix KYZ not directly scale with the input dimensionality, we replace the induc- ing derivatives from equation (7) with inducing directional derivatives. Rather than the triplet (z,, U4, Vu;) with Vu, having dimension D, each inducing point is now equipped with a set of p distinct directional derivatives (Z;, us @v .i, dv., u:) in the directions Vil...., Vip. We include the inducing directions V = [Vi... VMl Mpx 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 point 2; and some direction and another point and direction 2, and v, the directional-derivative covariance function is:
k(z,, 2j) V: (Zz,)
V k(z, 2j) Tv v.
0e D (2, Bz)
(11)
V 2.2, K (Zz,,
which is of size 2 rather than (D+ + 1) (D 1) with VV (:,-).
Given Mp inducing directions V, per each of the M 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, K Kxz, are:
Kzz V
VzKzz V
Kxz V2 X2 V VxKxz Vxx
K Kzz-
Kxz
(12)
VzZ V
22
the first of which has shape (p + 1) M(p + 1). Constructing K Kx z is inexpensive as we compute them directly from the directional derivative kernel (IT), rather than computing the full gradient kernel A and multiplying by the directions V which would incur cost of O D²).
Variational inference with this model is nearly identical to inference with full inducing gradients. We define variational posterior, this time over the M (p + 1) inducing values and directional derivatives:
g(u, @vu) = / (m,
(13)
where m E R M(--1) and S E R M(++ (pp+))xpp-1)) The model is trained by optimizing the variational ELBO
E (y, [log ply, Vy f. Vf) - (L ayu) [pp(u, uyu)).
(14)
over the variational parameters m. S, the inducing points Z, the inducing directions V. and the hy- perparameters. 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 R M(++1) Inference proceeds by computing qf, Vf) from (IO) by replacing the kernel matrices KX2 and KY with our directional derivative variants Kxz and Kzz.
Learning inducing directions The above algorithm requires the selection of set V 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] () However, for a variational GP with directional derivatives, sharing inducing directions does not improve scalability as the size of the kernel matrices K x 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 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
Derivative modeling with p < D. A key feature of this framework is that it allows for the use of 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 z; has the p = D canonical inducing directions V; = T, through which we recover the full variational GP with derivatives described in section 4.1
5

---

Complexity, For a minibatch size B, when learning directional derivatives per inducing point, the matrices Kxz and Zz become B x [p + 1) and M + 1) x M(p + 1) respectively. As result, the time complexity of variational GP inference with directional derivatives is pp") When using p directions per inducing point, this is computationally equivalent to running SVGP with p times as many inducing points. To counteract the additional matrix size, one may use the whitened formulation of variational inference [7 for GPs when computing equation 10 and use contour integral quadrature in [21].
Experiments
In this section we evaluate the empirical performance of variational GPs with directional derivatives. In sections and 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 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 mishapadidar -eeear/G---DDeiiaaaeeeee were--Devvaaivese-eeeee Variationnlainall Inference
and 5.2 we benchmark the performance of our method on small regression problems
g on single GPU. Code is available at hhttt::/ // //g thub. com/
5.1 Synthetic functions
In this section, we consider low-dimensional regression with derivatives test functions including Branin (2D), SixHumpameel (2D). Styblinksi-Tang (2D) and Hartmann (6D) from [32], a mod-
ified 20D Welch test function I (Welch-m)" and 5D sinusoid f(x) sin(2r (2 | |=)) (Sin-5). We compare variational GPs without derivatives (SVGP, PPGPR) to variational GPs with deriva- tives (GradSVGP, GradPPGPR), exact GPs with derivatives (GradGP), non-variational GPs with derivatives (DSKI), and variational GPs with = 2 directional derivatives per inducing point DDSVGP2DDPPGRR) Exact and variational GPs with full derivatives only tractable in low- dimensional settings due to the scalability issues mentioned in rections[aand 1. Therefore, to apply GradSVGP and GradPPGPR the 20D Welch-m function, onto low dimensional active subspace (2). similar to [5]. To show the limitation of GradS and GradPPGPR, we modified the Welch function have a low-quality low-dimensional active subspace.
first perform dimension reduction
In this mimeeninnsional setting, we find that variational GPs with directional derivatives, DSVGP2 and DPPGPR2, perform comparably to the methods that incorporate full derivatives (DKSI, GradSVGR GradPPGPR, GradGP); see Table [ In Figure [ we compare the negative log likelihood of each method as the inducing matrix size grows 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.
5.2 Implicit Surface Reconstruction
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 inducing directions for 1200 epochs and DSKI on 1 11606 noisy observations of 44818 locations and corresponding noise-free surface normals (gradients of the bunny level sets). DSVGP smoothly reconstructs the bunny and is comparable to DSKI.
5.3 Training Graph Convolutional Neural Networks with Bayesian Optimization
In this section, we demonstrate the full scalability of our approach by training the D = 4035 parameters of two layer graph convolutional neural network (GCN) 16 on the node classification
'The Welch test function has intrinsically 6D active space. We modified it to have a low-quality 6D active subspace and to show the limitation of GGP and GradPPGPR.

---

Branin
Camel
StyTang
Sin-5
Hartmann
Welch-m
RMSE
RMSE
RMSE
RMSE
RMSE
RMSE
NLL
NLL
NLL
NLL
NLL
NLL
(1e-3)
(Ie-3)
(1e-3)
(1e-1)
(1c-1)
(1e-2)
SVGP
1.45
-3.12
5.28
-2.95
3.64
-3.06
6.64
0.99
1.02
-0.69
16.20
-0.39
PPGPR
1.60
-3.21
6.46
-3.10
4.64
-3.17
4.35
0.35
3.02
-1.28
18.08
-0.56
GradGP
15.4
-0.87
25.1
-0.22
44.4
-0.82
2.59
-.23
0.50
-0.74
16.3
-0.38
GradSVGP
0.35
-3.65
2.09
-3.62
1.00
-3.65
4.85
2.31
2.08
0.59
18.94
42.82
GradPPGPR
0.67
-3.32
23.1
-3.14
2.91
-3.30
4.83
0.37
3,95
-1.16
18.92
-0.25
DSVGP2
0.29
-3.10
1.82
-2.50
0.86
-2.97
3.03
1.87
0.92
-0.75
3.74
-0.74
DPPGPR2
0.47
-3.32
8.43
-3.24
1.75
-3.31
4.30
0.05
-1.64
26.08
-0.71
DSKI
0.91
-4.47
3.85
-3.00
1.59
-4.74
N/A
N/A
N/A
N/A
N/A
N/A
Table Regression results on Branin (2D), SixHumpCamel (2D), Styblinksi (2D), Sin-5 (5D),
Hartmann (6D) and Welch-m (20D), each with 10000 training and 10000 testing points. Following
[5]. we train GradGP on 10000/(D + 1) points. The inducing matrix size is 800 for all variational
inducing point methods, while DSKI is trained 800 inducing points per dimension. See the supplementary material for error bars.
PPGPR
GradSVGP
GradPPGPR
DSVGP2
DPPGPR2
2.5
2.0
2.0
1.0
1.5
0.0
1.0
0.5
-1.0
0
-2.0
200
400
800
1200
200
400
800
1600
3200
Inducing matrix size
Inducing matrix size
Figure I: Negative Log Likelihood (NLL) for the various GPs when using different inducing matrix sizes to regress on Sin-5 (Left) and Hartmann (Right). DPPGPR2 often outperforms other methods due to its ability incorporate derivative information while only modestly increasing the inducing matrix size. DSKI is removed because it does not have comparable matrix size. In both figures, shaded regions correspond to standard errors.
Figure 2: Surface reconstruction of the Stanford bunny: (Left) Original surface, (Middle) DSVGP with 800 inducing points and directions, (Right) D-SKI with 303 inducing grid points.
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 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) s 2500 (4035 + 1) > 10 combined function and partial derivative labels. Any one of these datasets would be intractable to existing methods for training GPs with gradient observations.
7

---

For this experiment, we make no effort to modify the Bayesian optimization routine itself to account for the derivative information (e.g., as (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, 1, DPPGPR2 surrogate models. See supplementary material for p = 3 results. We additionally include traditional BO with the Lower Confidence Bound (LCB) [3] acquisition function, Adam [S] and random search. All algorithms were initialized with 400 random evaluations, and the TuRBO variants were with batch size of 20 and retrained over 150 steps. Figure Jummarizes results averaged over trials. We observe that TuRBO with DPPGPR significantly outperform traditional BO and other TuRBO variants. While all Bayesian optimization methods under- peder-perform compared to Adam, conjecture that this performance gap could be narrowed by incorporating the gradient information into the Bayesian optimization algorithm itself.
Random
Adam
BO
TuRBO TURBO-PPGPR
TURBO-DPPGPRI
TURBO-DPPGPNS
10a
1.0
10
Training accuracy
Training loss
0.8
100
0.6
10-1
0.4
10
500
1000
1500
2000
2500
500
1000 1500
2000
2500
Number evaluations
Number evaluations
Figure 3: GCN training on the Pubmed dataset: (Left) training loss and (Right) training accuracy. Averaged over trials for all optimizers. Shaded regions correspond to standard errors.
5.4 Stellarator Regression
In this experiment, we show the capacity of variational GPs with directional derivatives to extend GP regression with derivatives to massive datasets high dimensional settings. We perform regression on N 500000 function and gradient observations gathered from D) = 45 dimensional optimization objective function through the FOCUS code [39]: a code for evaluating the quality of magnetic coils for Stellarator, magnetic confinement based fusion device for generating renewable energy 12]. The dataset is available upon request.
We compare variational GPs with directional derivatives using DPPGPRI. DSVGP2, DPPGPR2) variational GPs without directional derivatives (SVGP, PPGPR). See supplementary material for p = 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 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 variational GPs with directional derivatives significantly enhance regression performance. Even the inclusion 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.
1.2 directions (DSVGPI,
5.5 Rover Trajectory Planning
The rover trajectory planning problem 35 is D 200 dimensional optimization problem with gradients. This experiment validates the use of variational GPs with directional derivatives in setting where no other methods can. We solve 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 myyplane. The rover trajectory is integrated over 100 steps at which forces in and y directions are applied the rover,
Bayesian optimization by leveraging derivative information

---

SVGP
Random
1.4
1100
Objective function value
80
DSVGP2
TuRBO
-1.6
900
TuRBO- -PPBGRL
DPPGPR1
TURBO-OPPGPRG
1.8
700
2.0
500
2.2
300
200
500
800
1000
1200
1400
400
800
1200
1600
2000
Inducing matrix size
Number of evaluations
Figure 4: Stellarator Regression (Left) and Rover (Right). Negative log likelihood of GP variants as the inducing matrix size increases for the 45, N 500000 Stellarator Regression experiment. Rover (Right) shows the value of the objective function over the course of optimization. In both figures, shaded regions correspond to standard errors.
making 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 a algorithms almost immediately.
5.6 Ablation Study Over Number of Inducing Directions
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 the performance of DPPGPR the GCN Bayesian optimization (Section 5.3). stellarator regression (5.4). and rover Bayesian optimization (55) experiments. We find that using p 1 is usually sufficient for capturing much
GCN
Stellarator
Rover
Loss
Accuracy
RMSE (Ie-2)
NLL Objective
(Ie-1)
(%)
PPGPR
4.31
85.28
2,43
-2.22
DPPGPR1
1.59
96.33
2.39
-2.27
258.6
DPPGPR2
0.39
99.67
2.31
-2.31
255.8
DPPGPR3 9.65
98.89
2.32
-2.31
253.0
Table 2: The effect of increasing the number of inducing directions p on the performance of variational GPs with directional derivatives the GCN Bayesian optimization, stellarator regression, and rover trajectory planning tasks. PPGPR was run with 1000 inducing points and DPPGPRp was run with 1000/ (p + 1) inducing points to ensure equivalent computational complexity. Increasing from zero to one shows a significant improvement, but performance improvement diminishes after that.
of the benefit of derivative information, and that increasing p improve performance but at diminishing returns. See supplementary material for more results on the GCN, Stellarator, and Rover tasks.
5.7 UCI Regression
Increasing the number of inducing points for Gaussian process models often results in diminishing returns on final model performance, with 500 This saturation is likely due in part to the ability of sparse Gaussian processes to represent variation
S 2000 often proving sufficient 99 26 13 21).

---

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, variational GP with directional derivatives is initialized with inducing directional derivatives as normal, but during training, minibatches of data only correspond function observations. In other words, rows of the matrix K x never correspond to partial derivative observations because there none, however columns of K x 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-Directee (D=20, N=53414). We compare variational GPs without derivatives (SVGP, PPGPR) to variational GPs with p = 1 directional derivative per inducing point (DSVGP1, DPPGPRI). I). To ensure that methods have equivalent time complexity, we fix the inducing matrix size to 800. Thus, DSVGPI/DPPGPR 1 use 400 inducing points, so that 400 x (p + 1) = 800, while SVGP/PPGPR use 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
Table 3: Variational GPs with no derivative modelling (SVGP, PPGPR) versus with (DSVGP1, DPPGPR1) on UCI benchmark datasets for which no derivative information is available. Bold entries indicate the best method in the column. See the supplementary material for error bars.
Discussion
Augmenting GPs with derivative information can significantly improve their predicitive capabilities but at a significant O(N cost. We introduce a novel method for achieving fully scalable scalable in N and D - GPs with derivative information by leveraging stochastic variational approximations. The resulting model reduces the cost of training GPs with derivatives to O(M3p²) time per iteration of training, where M < N and p < D. A practical limitation of our method is that M, p must be small enough for fast computations, which is not a reasonable assumption in very high dimensional problems. Through a series of synthetic experiments and a surface reconstruction task, we demonstrate the quality of our approximate model in low dimensional settings. Furthermore, we demonstrate the full scalability of our model through training a graph convolutional neural network using Bayesian optimization. in addition to performing regression on large scale Stellarator fusion dataset with derivatives. Lastly, we show that our methods can even have benefit in the regression setting when derivative information is available in the training set, by including new avenue to encode information. While this last result is a surprising benefit of GPs with derivatives, it is not well understood and is thus good direction for future study. While more accessible to practitioners and researchers for calibrating uncertainty estimates, the fundamental assumption that the data is drawn from a GP may be flawed, leading to poor uncertainty estimates and lack of robustness altogether. Researchers and practitioners should take care to understand the reliability of the GP model in their setting rather than relying faithfully on black-box approach.
method may make GPs
7 Acknowledgements
We acknowledge support from Simons Foundation Collaboration on Hidden Symmetries and Fusion Energy and the National Science Foundation NSF CCF- 193495,,, and NSF 1645643.
10

---

References
[1] Einat Neumann Ben-Ari and David M Steinberg. Modeling data from computer experiments: An empirical comparison of kriging with MARS and projection pursuit regression. Quality Engineering, 1944)3327-338, 2007.
[2] Paul G Constantine. Active subspaces: Emerging ideas for dimension reduction in parameter studies. SIAM, 2015.
[3] Filip de Roos, Alexandra Gessner, and Philipp Hennig. High- High-dimensional Gaussian process inference with derivatives. 38th International Conference on Machine Learning. 2021.
[4] Dheeru Dua and Cascy Graff. UCI Machine Learning Repository. Irvine, CA: University of California, School of Information and Computer Science, 2017.
[5] David Eriksson, Kun Dong, Fric Hans Lee, David Bindel, and Andrew Gordon Wilson, Scaling Gaussian process regression with derivatives. In Advances Neural Information Processing Systems, volume 31, 2018.
[6] David Eriksson, Michael Pearce, Jacob R Gardner, Ryan Tumner, and Matthias Poloczek. Scalable global optimization local Bayesian optimization, Advances in Neural Information Processing Systems, volume 32, 2019.
Jacob Gardner, Gcoff Pleiss, Ruihan Wu, Kilian Weinberger, and Andrew Wilson, Product kernel interpolation for scalable Gaussian processes. International Conference on Artificial Intelligence and Statistics, pages 1407-1416. PMLR, 2018.
[8] Jacob Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, and Andrew Gordon Wilson. Gpy- torch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. In Advances Neurul Information Processing Systems, volume 31, 2018.
[9] James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data. In Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence, page 282-290, Arlington, Virginia, USA, 2013. AUAT Press.
[10] James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification. In Artificial Intelligence and Statistics, pages 351-360. PMLR, 2015.
[11] Matthew D Hoffman, David M Blei, Chong Wang. and John Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14(5), 2013.
[12] Lise-Marie Ambert-Gerard, Elizabeth J Paul, and Adelle Wright. An introduction to stellarators: From magnetic fields to symmetries and optimization, arXiv preprint arXiv: 1908,05360 2019.
[13] Martin Jankowiak, Geoff Pleiss, and Jacob Gardner. Parametric Gaussian process regressors. Interna- tional Conference on Machine Learning, pages 4702-4712. PMLR, 2020.
[14] David Jones, David Stenning. Eric Ford, Robert L Wolpert, Thomas Loredo, Christian Gilbertson. and Xavier Dumusque. Improving exoplanct detection power: Multivariate gaussian process models for stellar activity. arXiv preprint 171100 2017.
[15] Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. Proceedings of the 3rd International Conference on Learning Representations, 2014.
[16] Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In Proceedings of the 5th International Conference on Learning Representations, 2017.
[17] Alexander Graeme de Garis Matthews. s. Scalable Gaussian process inference using variational methods. PhD thesis, University of Cambridge, 2017.
[18] Jonas Mockus. Bayesian approach to global optinization: theory and applications, volume 37. Springer Science & Business Media, 2012.
[19] Jorge Nocedal and Stephen Wright. Numerical optimization. Springer Science & Business Media, 2006.
[20] Emanuel Parzen. Stochastic processes. SIAM, 1999.
[21] Geoff Pleiss, Martin Jankowiak, David Eriksson, Anil Damle, and Jacob R Gardner. Fast matrix square roots with applications to Gaussian processes and Bayesian optimization. In Advances in Neural Information Processing Systems, volume 33, 2020.

---

[22] Joaquin Quinonero- Candell and Carl Edward Rasmussen. A unifying view of sparse approximate Gaussian process regression. The Journal of Machine Learning Research, 6:1939-1959, 2005.
[23] Carl Edward Rasmussen, Gaussian processes in machine learning. In Summer school on machine learning. pages 63 71. Springer, 2003.
[24] CE. Rasmussen and CKI. Williams. Gaussian Processes for Machine Learning. MIT Press, Cambridge, MA, USA, 2006.
[25] Herbert Robbins and Sutton Monro. A stochastic approximation method. The Annals of Mathematical Statistics, pages 400-407, 1951.
[26] Hugh Salimbeni, Ching-An Cheng, Byron Boots, and Marc Deisenroth. Orthogonally decoupled variational Gaussian processes. In Advances in Neural Information Processing Systems, volume 31, 2018.
[27] Peter Schulam and Suchi Saria. Reliable decision support using counterfactual models. In Advances in Neural Information Processing Systems, volume 30, 2017.
[28] Prithviraj Sen, Galileo Namata, Mustafa Bilgic, Lise Getoor, Brian Galligher, and Tina Fliassi-Rad. Collective classification in network data. Al Magazine, 29(3):93-93, 2008.
[29] Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pscudo-inputs. In Advances in Neural Information Processing Systems, volume 18, 2006.
[30] Jasper Snoek, Hugo I Larochelle, and Ryan P Adams. Practical Bayesian optimization of machine learning algorithms. In Advances in Neural Information Processing Systems, volume 25, 2012.
[31] Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias W Seeger. infereitionthhooreticc regret bounds for Gaussian process optimization the bandit setting. IEEE Transactions Information Theory, 58(5):325033265, 2012.
[32] Sonja Surjanovic and Derek Bingham. Virtual library of simulation experiments: Test functions and datasets. Retrieved May 27, 2021, ftom/att::///ww.cc. sfu. ca/ "ssurjano" 2013.
[33] Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Artificial intelligence and statistics, pages 567-574. PMI.R, 2009.
[34] Selvakumar Ulaganathan, Couckuyt, Tom Dhaene, Joris Degroote, and Eric Laermans, Performance study of gradient-enhanced kriging. Engineering with computers, 32(1):15-34, 2016.
[35] Zi Wang, Clement Gehring, Pushmect Kohli, and Stefanic Jegelka. Batched large-scale Bayesian optimiza- tion in high-dimensional spaces. In International Conference on Artificial Intelligence and Statistics, pages 745-754. PMLR, 2018.
[36] Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured Gaussian processes (KISS-GP). In International Conference on Machine Learning, pages 1775-1784. PMIR, 2015.
[37] Jian Wu, Matthias Poloczek, Andrew Gordon Wilson, and Peter Frazier. Bayesian optimization with gradients. In Advances in Neural Information Processing Systems, volume 30, 2017.
[38] Zhengkun Yi, Roberto Calandra, Filipe Vciga, Herke van Hoof, Tucker Hermans, Yilei Zhang, and Jan Peters. Active tactile object exploration with gaussian processes. In 2016 JEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 4925-4930. IEEE, 2016.
[39] Caoxiang Zhu, Stuart R Hudson, Yuntao Song, and Yuanxi Wan. New method to design stellarator coils without the winding surface. Nuclear Fusion, 580 ))016008 2017.
12


Reviews:
[{"review_id": "yyBSBXu3cSG", "reviewer": "Reviewer_GftX", "summary": "   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.\n\n   Strengths\n   + 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.\n   + The proposed approach is clear and builds on established techniques.\n   + The computational benefits demonstrated in an extensive set of experiments are substantial and convincing.\n\n\n   Weaknesses\n   - Experimental results are without error bars.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "   1) Clarity\n\n   The paper is well written. The notation is concise and the technical content is easily accessible. The experiments are well described.\n\n   2) Originality\n\n   The proposal to use directional derivatives as variational parameters to render GP inference with derivatives more scalable seems novel.\n\n   3) Significance\n\n   Although the setting of GP regression with high N, high D and derivative information available sound like a niche application, the presented experiments show that these settings do exist and could not be tackled before. Beyond this, the scope of the method seems rather limited.\n\n   4) Reproducibility\n\n   The data and the code to run the experiments will be released upon acceptance. Although not available during the review process, it should be possible to reproduce the results.\n\n   5) Empirical analysis\n   \n   The manuscript contains a wide range of experiments illustrating various aspects of the method on a number of data sets along with comparative results relative to the relevant baseline methods.\n\n   6) Minor Issues and Typos\n\n   line 197: DSKIP rather than DKIP\n   line 206: learning tasks, we perform\n   Table 1, caption: Hartmann rather than Hatrmann\n   Figure 1, caption: ability to incorporates\n   line 256: In this experiment, we show\n   line 283: in Figure 4\n   line 312: reasonable assumption in very high", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"ethical_concerns": "No concerns.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "uquZxEU-kC", "reviewer": "Reviewer_TvZh", "summary": "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.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall, the paper checks all the boxes in terms of writing and narrative framing, presenting a clean linear story. The contribution is important, given that it brings a general approach to augment variational GP models with derivative information via PPGPR to allow heteroscedasticity. The choice of using directional derivatives is natural in hindsight, making this contribution even more important.\n\nOne setting which seems to be missing is the choice of using a fixed set of inducing directions for all inducing inducing points, instead of picking unique directions for every inducing point. Intuitively, is this choice the natural first step? Would the flexibility of the model suffer if we started with directions shared across all inducing points?\n\nSection 5.6 on using DPPGPR even when no derivative information is available is interesting, but not entirely unclear. How does the model in Eq. $(6)$ change if we do not have derivative information. This does not appear to have been spelled out clearly, or isn't obvious.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "KQe0qHlb6ww", "reviewer": "Reviewer_VwPg", "summary": "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. \n\n", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The directional derivative idea to map high-dimensional derivative kernels down to a scalar is novel and seems to work practically well. However, I'm still a bit skeptical about the choice of directions which could have a potential effect I assume. But this is not analyzed in the paper at all. The paper needs to provide better intuitions, especially in the experiment part to convince readers about that. The paper needs more edits to make it clear and easy to follow.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "1. mu_f(x_i) is not defined before or after the use in eq. 4.\n\n2. How to construct kernel matrices in eq. 12 from eq. 11?\n\n3. How to select the directions for v? Especially with large D. I assume mapping high-dimensional derivatives down to a scalar with v and p=1,2 would be quite tricky. \n\n4. I don't understand this sentence very much \"as this would confound the performance improvements achieved by higher fidelity modeling by incorporating derivative information.\" Could you explain it? I'm confused that if the derivative information is not considered in BO, what is \\nabla f in this example? \n\n5. Why DSKI have N/A in table 1? \n\n6. In sec 5.2, how many inducing points, i.e., M, are used with DSVGP? DSKI on 11606 observations of 34818 locations. What are 11606 and 348181?\n\n7. Since DSKI claims that it could scale to high dimensions, it would still be helpful to include DSKI for applications with large D, like sec 5.3 and 5.4.\n\n8. Flip the order of the two figures in figure 4 since stellarator is introduced before rover.\n\n9. It's a bit confusing to arrange the experiments with the current order: 5.1 regression, 5.2 surface reconstruction, 5.3 BO, 5.4 regression, 5.5 BO, 5.6 regression without derivative information. It feels like the comparable methods switch between different groups. This makes the logic hard to follow. It might be helpful to have a table summarizing what application the experiment is and whether it's large N or large D or both.\n\n\n\n\n\n\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2 hours"}, "extra_text": {}}, {"review_id": "Ei0jrgeQ_sF", "reviewer": "Reviewer_g1G4", "summary": "This paper presents a new method to enable Gaussian processes with derivatives to deal with large and high-dimensional datasets. Several synthetic and real-world experiments are conducted in evaluating the proposed method. ", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality and Significance: To be honest, I am a little bit disappointed when reading the technical part and experiments. Since two key ingredients discussed in the introduction:\n1. heteroscedastic noise\n2. principal to select $p$\n\nI am expecting the authors can address 1. what is the advantages of heteroscedastic noise rather than i.i.d noise? 2. how will the model perform with different $p$. \n\nIt is natural that selecting a small number of $p$ can reduce the time complexity, however, the authors did not really mention how to determine $p$ for different tasks? PCA, t-sne, or even gplvm can do this dimensionality reduction right? I am expecting more discussion relating to derivative and dimensionality reduction. This should be the golden apple of this paper.\n\nMoreover, Sec 4.1 seems to be simple adaption of VSGP to VSGP with derivatives, is there any difficulty in adopting VI? Any other benifits in adopting VI other than reducing time complexity?\n\nThe two aspects above are not extensively discussed and evaluated in the paper which makes the contribution incremental. \n\nQuality & Clarity: The writing of this paper is really smooth and easy to follow. I really appreciate the efforts the authors put in polishing this paper. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "refer to the detailed comments above.", "review_text": "", "extra_scores": {"ethical_concerns": "N.A.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper addresses the problem of scalable inference in Gaussian process (GP) regression with derivative information for the general case where both the number of observations (N) and the dimensionality (D) are large. While previous work has addressed the large-N-low-D and low-N-large-D regimes independently, the paper proposes the use of inducing directional derivatives (which play a similar role to that of the inducing variables in standard variationally sparse GP models) rather than full derivatives to derive a scalable variational inference algorithm. A comprehensive set of experiments (on synthetic data, implicit surface reconstruction, training graph convolutional networks with Bayesian optimization, large-scale regression with derivative information, Rover trajectory planning and standard GP regression on UCI datasets) show the benefits of the proposed approach with respect to previous methods. 

I believe the contribution of this paper to be significant to the NeurIPS/GP community as it develops a new method for settings where previous work cannot be applied. The reviewers raised several concerns regarding the choice of inducing directions and the use of different directions per inducing point. I believe the authors have addressed these successfully. A limitation of this work is the very small number of directional derivatives (p) used in the experiments (p=1, 2). As the complexity of the proposed algorithm scales as (M^3 p^3), exploring larger p values will necessarily lead to fewer inducing variables (M). 

Decision: Accept (Poster)