Dual Parameterization of Sparse Variational Gaussian Processes
Authors: Vincent ADAM
Paul Edmund Chang
Mohammad Emtiyaz Khan
Arno Solin
Venue: NeurIPS 2021 Poster (2021)
Keywords: Gaussian processes
sparse variational inference
natural gradients

Abstract: Sparse variational Gaussian process (SVGP) methods are a common choice for non-conjugate Gaussian process inference because of their computational benefits. In this paper, we improve their computational efficiency by using a dual parameterization where each data example is assigned dual parameters, similarly to site parameters used in expectation propagation. Our dual parameterization speeds-up inference using natural gradient descent, and provides a tighter evidence lower bound for hyperparameter learning. The approach has the same memory cost as the current SVGP methods, but it is faster and more accurate.

Dual Parameterization of Sparse Variational Gaussian Processes
Vincent Adam" Aalto University Seconndiind.ai Espoo, Finland Cambridge, UK vincent adam@aalto fi
Paul E. Chang" Aalto University Espoo, Finland paul. changgaalto. fi
Mohammad Emtiyaz Khan RIKEN Center for AI Project
Arno Solin
Aalto University
Tokyo, Japan emtiyaz. wwankriken.n
Espoo, Finland arno ceno.ooli@aalto
Abstract
Sparse variational Gaussian process (SVGP) methods are common choice for non-conjugate Gaussian process inference because of their computational benefits. In this paper, we improve their computational efficiency by using dual parame- terization where each data example is assigned dual parameters, similarly to site parameters used in expectation propagation. Our dual parameterization speeds-up inference using natural gradient descent, and provides tighter evidence lower bound for hyperparameter learning. The approach has the same memory cost the current SVGP methods, but it is faster and more accurate.
Introduction
Gaussian processes (GPs, [31]) have become ubiquitous models in the probabilistic machine learning toolbox, but their application is challenging due to two issues: poor O(n²) scaling in the number of data points, n, and challenging approximate inference in non-conjugate noon-Gaussian) models. In recent years, variational inference has become the go-to solution to overcome these problems, where sparse variational GP methods [35] tackle both the non-conjugacy and high computation cost. The computation is reduced to (nnz² by using small number of < inducing- locations. For large problems, SVGP framework a popular choice as enables fast stochastic training, and reduces the cost to O(m + mmm??) per step, where nı is the batch size [10, 11].
is a common practice SVGP to utilize the standard mean-covariance parameterization which requires O(m²) memory. Inference is carried out by optimizing an objective C(g) that uses Gaussian distribution & parameterized by parameters E = (m, L) where is the mean and L is the Cholesky factor of the covariance matrix. We will refer to the SVGP methods using such parameterization as the VVPV methods. The advantage of this formulation is that, for log-concave likelihoods, the objective is convex and cradient-based optimization works well [2]. The optimization can further be improved by using natural gradient descent (NGD) which is shown to be less sensitive to learning rates [33]. The NGD algorithm with SVPV parameterization is currently the state-of-the-an and available in the existing software implementations such as GPflow [26] and GPyTorch [9].
An alternate parameterization to the SVP parameterization is the one where every likelihood is assigned two of parameters which require O(n) memory. Existence of such parameterizations was initially shown by Csató and Opper [8] for general GP models, and later on extended to variational
Both authors contributed equally.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

f-SVGP (ours)
1-SVG (ours) gVVV
g-SVGP
ELBO, Lels" (0.r)). .)
Q-SVGP (whitened)
8/ (magnituds)
-63.4
Tighter hound
gradient match
ELBO and
10
108 66
0
0.2
0.4 0.6 0.8
0.8
0.9
1
01 (length-scale)
Figure 1: The dual parameterization gives a tighter bound compared to standard SVGP (whitened/ unwhitened) as shown on the left for GP classification tasks and varying kernel magniude 0. The tighter bound helps take longer steps, speeding up convergence of hyperparameters, as shown on the right for the banana classification task with coordinate ascent w.r.t. @ and variational parameters &.
objectives for GPs [28, 29], and also to latent Gaussian models by using Lagrangian duality [20, 16]. Due to this later connection, we refer to this parameterization as the dual parameterization where the parameters are the Lagrange multipliers, ensuring that the marginal mean and variance of each latent function is consistent to the marginals obtained by the full GP [16]. Expectation propagation (EP, (27]) too naturally employ such parameterizations, but by using site parameters. Although unrelated to duality, such methods are popular for GP inference, and due to this connection, we will refer to the methods using dual parameterization as 1-SVGP (the letter "t" refers the sites). To the best of our knowledge, the dual parameterization for SVGP has only been used to speed up computation and inference for the specific case of Markovian GPs [4, 38].
The main contribution of this work is to introduce the dual parameterization for SVGP and show that it speeds up both the learning and inference. For inference, automatically obtained through a different formulation of NGD, written in terms of the expectation parameters [17, 15. 18]. The formulation is fast since it avoids the use of sluggish automatic differentiation to compute the natural gradients. We also match the typical O(m²) memory complexity in other SVGP methods by introducing a fied parametrization. For learning, we show that the dual parameterization results in a tighter lower bound to the marginal likelihood, which speed-up the hyperparameter optimization (see Fig. 1). We provide extensive evaluation on benchmark data sets, which confirms our findings. Our work attempts to revive the dual parameterization, which was popular in the early 2000s, but was somehow forgotten and not used in the recent SVGP algorithms.
show that the dual parameters
Background: Variational Inference for Gaussian Processes Models
Gaussian processes (GPs, [31]) are distributions over functions, commonly used in machine learning endow latent functions generative models with rich and interpretable priors. These priors can provide strong inductive biases for regression tasks in the small data regime. GP-based models are the ones that employ GP prior over the (latent) functions f(-) GP(p(-), s(,, ')), where the prior is completely characterized by the mean function μu(-) and covariance function -") We denote f(+) as a function but occasionally simplify the notation to just f in the interest of reducing clutter. Given a data set D = y) = {(x, 3/. F of input-oulput pairs, we denote by f the vector of function evaluations
The function evaluation at x; are passed through likelihood functions to model the outputs y R, e.g. by specifying p(y f) == [ P(v: f:). Prediction at a new test input x, is obtained by computing the distribution p(f(x.) D,x.), For Gaussian likelihoods N(v: f., a²), the predictive distribution is available in closed form as Gaussian distribution (x.(x.) MGGPR (x.). UGPR (xPR(x.)) with mean and variance defined as follows:
the inputs {f(x;)
MGPR (x.) == kT, (K + o II,) -ly, and UGprR (x.R((x)) ). == K.n kT, (K + +o
(1)
where kr. is a vector of s(x,, x;) as the 7th element for x; E X, KR is an n x a matrix with x(x,, Xy) as the 2h entry, and Rew x()..

---

2.1 Variational Expectation-Maximizaale for GPs with Non-Conjugate Likelihoods
For non- non-Gaussian likelihoods, the posterior and predictive distributions are no longer Gaussian, and we need to resort to approximate inference methods. Variational inference is popular choice because it allows for fast posterior approximation and hyperparameter learning via stochastic training [35, 10]. Denoting kernel hyperparameters by 0 and the corresponding GP prior by (f( the posterior distribution can be written as (f yy) po((f) p(y f)/po(y), where po(y) marginal likelihood of the observations. We seek to approximate po y) = qr (f) by Gaussian distribution whose parameters can be obtained by optimizing the following evidence lower bound (ELBO) to the log- marginal likelihood,
po (f, y)df the
log (y) 2 (gr, 0) = -. (F.) [logp(s [A) - DKLlgr le
(2)
For the variational approximation, is a standard practice to choose the mcan-covariance parameteri- zation, denoted by (m, S). is also common uniquely determined for a covariance matrix. The multivariate normal distribution is part of the expo- nential family [37]. i.e. it's probability density function take the form p(x) = cxp(1) T(x) a(n)). with natural parameters 77 - SS11/2
the Cholesky factor L instead [3] since it is
and sufficient statistics T(x) [x, xx This
natural parameterization is also a common choice, along with the associated expectation parameteriza- tion μ E3, [x) (m, S + mm ). A final choice of the parameterization, called the whitened parameterization [36], uses a variable v ~ N(v; my, Sv) along with the transformation f Lv, to LS.LT.
parameterize m Lm, and
In the following, we will consider the ELBO with several parameterizations and, to make the notation clearer, we will indicate the parameterization used with subscript with L. For example, we may have (g 0) C, (77. 0) C (4, 0), which are all clearly equal due to unique mapping between the parameterizations; see [25] for more details on the maps.
The ELBO can be optimized by using a variational expectation expectation-maxiimzzttoon (VEM) procedure, where we alternate between optimizing variational parameters, say &, and hyperparameters O,
E-step: & + arg maxg Cg(E. 0.),
M-step: 0,+1 < arg max Le 0),
(3)
where denotes the iterations, and we have explicitly written the dependence of optimal parameter (0.) as a function of the old parameter 0,. Both and and M-steps can be carried out with gradient
descent, for example, using an iteration of the form ( + ((h) + P Te (k) 0,) for E-step,
which would ultimately converge 37. similar iterative method can be used for the M-step.
2.2 Inference via Natural-Gradie Gradient Descent (NGD)
A popular strategy for the E-step is to use natural gradient descent where we replace the gradient by the one preconditioned using the Fisher information matrix F(E) of gr()). We denote natural gradients by VeLe 0) F(E) V Ce 0), to get the following update,
(((+1) + {! + Pk 5 (K) 0.).
Such updates can converge faster than gradient descent [33, 17, 18], and at times are less sensitive to the choice of the learning rate Ps due the scaling with the Fisher information matrix. The implementation simplifies greatly when using natural parameterizations,
(k+1) + Ti (M) T
+ Px VuC (k) 0.).
(5) because = Vu, 0). that is, the natural gradients with respect n are in fact the gradients with respect to the expectation parameter μs [18]. In the remainder of the paper we will frequeetly use this property, and refer to the natural gradient with respect to by the gradients with respect to expectation parameterization μ. The VEM procedure with NGD has recently become a popular choice for sparse variants of GPs, which explain next.
2.3 Sparse Variational GP Methods and Their Challenges
Inference GP models, whether conjugate or non-conjugate, suffers from an Onn)) computational bottleneck required to invert the posterior covariance matrix. A common approach to reduce the
3

---

computational complexity is to use sparse approximation relying on small number m < representative imputs, also called inducing inputs, denoted by Z := (Z)
Zz..... Zm) [34, 7, 30, 39].
Sparse variational GP methods [35, 10, 11, 5, 6, 32] rely on Gaussian approximation (u) over
the functions u = (z1), (z2), (()).....,,(.)m) f(mm) to approximate the posterior over arbitrary locations,
qu (u.0((())) - f (-(((-) |u) (u)
(6)
where po (f(-) [u) is a conditional of the GP prior. For example, for a Gaussian (u(u) = Nuu; mu Su), the posterior marginal of fi = (xx) takes the following form,
qu.o(fi) = N (ft a mu, Re a. (Kuu - Su)a:).
where a, = Kuukui with K the prior covariance evaluated Z. and kui m-length vector of n(zz, x,), Vj. The parameters Eu = (m, Su) can be learned via an ELBO similar to Eq. (2),
= Ei= 11 E. (f.) [log If)] - Dxal9. (u) || (u)).
(8)
The variational objective for such sparse GP posteriors can be evaluated at a cost O(nm² + m²). and optimization can be performed in O(m²+n,mm² per iteration via stochastic natural- -auuurlnggrrdnttttt methods with mini-batch size nh [10]. This formulation also works for general likelihood functions [1 I]. This and the low computational complexity has lead to wide adoption of the SVGP algorithm. It currently the thr for sparse variants of GP and is available in the existing software implementations such as GPflow [26] and GPyTorch [9].
Despite their popularity, the current implementations are cumbersome and there is plenty of room for improvements. For example, the methods discussed in Salimbeni et al. [33] (see App. D for a summary) rely on the mean-covariance parameterization and NGD is performed in 17-space. However, the natural gradients are implemented via chain rule: Ve CE 0) utilizing the gradients in the E-space, which requires computation of additional Jacobians and multiplication operations. Since other operations are done via covaricove parameterization, we need to go back and forth between 1) and g, which further increases the cost. In addition, many existing implementations currently compute the natural gradient of the whole ELBO, including the KL term which is not required; see Khan and Rue [19, Sec. 2.2]. Finally, the M-step is dependent on the choice of the parameterization used in E-step and can affect the convergence speed. To the best of our knowledge, this has not been investigated in the literature.
In what follows, we argue to use a dual parameterization instead of the usual mean-covariance parameterization, and show that this not only simplifies computations of natural-gradients but also gives rise to tighter bound for hyperparameter learning and speed-up the whole VEM procedure.
3 The G-P Method: Dual-Parameter Based Learning for GPs
We start with a property of the optimal 97 of Eq. (2). Khan and Nielsen [18, Eq. (18)) show, that o7 can be parameterized by 2D vectors A (Ai i X5.) used in site functions (fs).
of (f) x [[] (A)T(()))
where X = VV (n log f().
(9)
The vectors X are equal to the natural gradient of the expected log-likelihood, where the expectation taken with respect to the posterior marginal (f:) = N(Jmm ( with S% as the i'th diagonal element of S*, and the gradient taken with respect to its expectation parameter μ4; (m, m + S and evaluated at μ.a. The vector T(f:) (fi, f?) are the sufficient statistics of a Gaussian. The oi uses local unnormalized Gaussian sites & (ff), similarly to those used the Expectation Propagation (EP) algorithm [27]. The difference here is that the site parameters X; are equal to natural gradients of the expected likelikel which are easy to compute using the gradient with respect to p4s
The parameters x can be seen as the optimal dual parameters of a Lagrangian function with moment-matching constraints. We can show this in two steps:
1. For each p(u: [i), we introduce a local Gaussian T6 (fz; (i:) with expectation parameters ñ. 2. Then, we aim to match मú with the marginal moments μi of the global Gaussian gr (f: 11).

---

This is written below as a Lagrangian where the middle term (shown in red) 'decouples' the terms using the local Gaussians from those using the global Gaussian,
Laagrange (1., pr. A) > (f.i.) Dog p(ys) f(]l -i=101 pb) Dxulgr(f; (1) pe (f)
(10)
The parameter X; is the Lagrange multiplier of the momentt-mtchinng constraint fı; 14, (here, X and j denote the sets containing all X; and f4:). The optimal X; is equal to X shown in Eq. (9).
We can show this by, first setting the derivative with respect to X; to 0, finding that the constraints j μ are satisfied. Using this and by setting the derivative with respect to fa; to O, we see that the optimal X is in fact equal the natural gradients, as depicted in Eq. (9).
The optimal natural parameter of of (f), denoted by 17" has an 'additive' structure that, as we will show, can be exploited to speed-up learning. The structure follows by setting derivatives w.r.t μ to O,
DKL par (f:; pz" pe (f)] = X =
'¹' = (0) X"
(11)
The second equality is obtained by using the result that Va.DKL μ1) po( = (0) where 7 and (0) are natural parameters of qf;; μ4) and (f) respectively [19, Sec. 2.2]. The final result
follows by noting that the right-hand side the natural parameter of (f) from Eq. (9). This implies that the global (f::: μ²) of (f) is also equal to the optimal approximation, as desired.
The Lagrangian formulation is closely related the maximum entropy principle [13] which forms the foundations of Bayesian inference [14]. Through moment matching, the prior is modified to obtain posterior approximations that explain the data well. Since X are the optimal Lagrange multipliers, they measure the sensitivity of the optimal of (f) to the perturbation in the constraints, and reveal the relative importance of data examples. Therefore, the structure of the solution 17 shown in Eq. (11) is useful for estimating 0. The additive structure can be used to measure the relative importance of the prior to the dual parameters X Our main idea is to use the structure to speed-up learning for SVGPs.
Eq. (1 can be rewritten in terms of the rrea-covarianee parameterization, to gain further insight about the structure. To do so, we use Bonnet and Price's theorem, and rearrange to get the following (see Eqs. 10 and 11 in [19] for a similar derivation),
m = =Kra" (S*)- K + diag(8"), where 3 is vector of B = ( [- p(1: fa]] (13)
where a is a vector of a: = (x) V log p(u: [f)).
(12)
The variables of and 33 can be easily obtained by using the gradient and Hessian of the lokelikeldhood, and using those we can get X (8)m, + af, (.))
Several other works have discussed such parameterizations, although our work is the first to connect it to natural gradients as the optimal Lagrange multiplier. The representation theorem by Kimeldorf and Wahba [22] is perhaps the most general result, but Csató and Opper [8] were the first to derive such parameterization for GPs; see Lemma in their paper. Their result is for exact posteriors which minor difference there is that their parameterizations use the integrals of likelihoods (instead of log-likelihoods) with respect to the GP prior (instead of the posterior), but we can also express them as Eq. (9) where of (f) is replaced by the true posterior (f y
intractable while ours is for Gaussian approximations and easy to compute.
Parameterization of the variational posterior similar to ours are discussed in [28, 29], but the one by Khan et al. [20] is the most similar. They establish the first connection to duality for cases where ELBO convex with respect to the mean-covariance parameterization. Khan [16] extends this to non- convex ELBO using the Lagrangian function similar to ours, but written with the covariccee parameterization to get the solutions shown in Eqs. (12) and (13). As shown earlier, their (a, (3) parameterization is just reparameterization of our X parameterization. Here, we argue in favour of our formulation which enables the reformulation in terms of site functions in Eq. (9) and also allows us to exploit the 'additive" structure in Eq. (11) to speed up hyperparameter learning. The mean-covariance parameterization does not have these features.
3.1 Improved Objective for Hyperparameter Learning
We will now discuss a method to speed-up VEM by using the dual parameterization. The key idea is to exploit the form given in Eq. (11) to propose a better objective for the M-step.
Standard VEM procedures, such as those shown in Eq. (3), iterate pairs of E and steps which we here describe in the context of the dual parameterization. In the E-step, starting from hyperparameter
5

---

9,, the optimal variational distribution (f) maximizing the ELBO in Eq. (2) is computed. For the dual parameterization, we get the optimal variational parameters ni (e.) + x; Here, the subscripts I. in ny and X indicate the dependence of the E-step iterations on θ., while 17 (0.) indicates a direct dependence of the prior natural parameter over 0,. The standard M-step would then be to use ni in the ELBO in Eq. (3) as shown below, while we propose an alternate procedure where the prior 770 (0) is left free (shown in red):
Standard M-step: 0,+1 = arg ming (7o + x; 0)
(14)
Proposed M-step: 0N+1 arg ming (, (0) + Xi 0)
(15)
This proposed objective is still lower bound the marginal likelihood and it corresponds to the ELBO in Eq. (2) with a distribution whose natural parameter is r): (0) (7))) + X and thus depends on O. We denote this distribution by gr(f: (f: n (0)) The ELBO is different from the one the distribution obtained after the E-step, with natural parameter fi which is independent of 0. We denote this distribution by qr(f; Clearly, at 0- O. both objectives match, and do their gradient with respect to 0, but they generally differ otherwise. We argue that the proposed M-step could lead to a tighter lower bound; see Fig. 1 for illustration.
In the standard M-step, the dependency of the bound on 0 is only via the KL divergence in Eq, (2). In the M-step we propose, this dependency is more intricate because the expected log- lieliikl also depend on O. Yet, as we show now, it remains simple to implement. The lower bound in the proposed M-step takes form where an existing implementation of GP regression case can be reused.
I p(ys fi)potf
(7) + Ni (Fv)
log
log ())
+ cc0), (16)
IP 4 (fi)pott)
2,
where c(0) = -i=1 (f log plu: (A) and log (0) is the log-partition of gr (f: (0)).
(7,)
log (0) = - l0g( 1 log| diag(B,)) (+) 1y [diag((;)) + [+) y.
(17)
Here, y a vector of To = - A2i and we have explicitly written (0) to show its direct dependence on the hyperparameter O. The gradients of (0) can be obtained using GP regresssion code, while the gradient of c(O) can be obtained using standard Monte-Carlo methods. A similar lower bound was originally used in the implementation" provided by Khan and Lin [17], but they did
not use it for hyperparameter learning.
For GP regression, we recover the exact uog-magginal likelihood log pe (y ()) for all values of 0. Indeed X; = /v= - c(0)= O. This also gives us T7+ y/s and 3: 1/o², and we get log (0) log po (y DD)
1/(20²)), which means that the sites exactly match the likelihood terms so
For non- conjugae problems, we found it to be tighter bound than the standard ELBO (Eq. (14))) which could speed-up the procedure. This is illustrated in Fig. (top row) where the proposed ELBO is compared to two other parameterizations (mean-covariance and whitened) for many values of 0, Oaa- We see that the maximum value (shown with a dot) remains rather stable for the proposed method compared to the other two. This is as expected due to Eq. (15) where we expect the solutions to become less sensitive to 02 because we have replaced no (0.) by (0)). The bottom row in Fig. shows the iterative steps (0.+1, 0z) for a few iterations, where we see that, due to the stable solutions of the new ELBO, the iterations quickly converge to the optimum. Exact theoretical reasons behind the speed-ups are currently unknown to us. We believe that the conditioning of the ELBO is improved under the new parameterization. We provide some conditions App. A under which the new ELBO would provably be tighter.
3.2 Faster Natural Gradients for Inference Using the Dual Paramterization
So far, we have assumed that the both E and M steps are run until convergence, but it is more practical use stochastic procedure with partial E and M steps, for example, such those used in [10, 12]. Fortunately, with natural-gradient descent, we can ensure that the iterations also follow the same structure as that of the solution shown in Eq. (9) and Eq. (11). Specifically, we use the method of
"See

---

Khan and Lin [17], expressed in terms of the dual parameters À and natural gradients of the expected
log-likelihoods g: (k) Vv, of (1. Loggloopuu [fi] (see also [19, Sec. 5.41),
(K+1) Of (f) x (f( IT- e (2 (R (((A) where A x(++1) = (1 - -rr) + rxg. (R)
(18)
(A)
The convergence of these iterations guaranteed under mild conditions discussed in [21]. The natural parameter of Of (f) at iteration k can be written in terms of 0 as follows,
[8)
(e) + X (A)
(19)
and the expectation parameters μi (R) required to compute the natural gradients of the expected log-likelihood, can be obtained by using map from the natural parameter 17(
The updates hold for any 0 and can be conveniently used as 0 = O: at the E-step of the 4h EM iteration. We name I-VGP the EM-like algorithm with 1) an E-step consisting of the natural gradient updates of Eq. (18), and, 2) the proposed M-step introduced Eq. (15).
4 The t-SVGP Method:
Based Inference for SVGP
We now extend the dual parameterization based stochastic VEM procedure to the SVGP case and refer the resulting algorithm (-SVGP. The optimality property shown in Eq. (9) shared by the ELBO given in Eq. (8). That is, we can express the optimal qu, (u) in terms of n 2D parameters X;
La (u) x po(u) (aaT u)) where X; = - (r p(ue f:)]-
(20)
(2(m)
The difference here is that the site parameters use the sufficient statisties T(a, u), defined via the projections a, = Kuiku Kui- The natural gradients of the expected log-likelihood are computed by using the marginal qu (f:) defined in Eq. (6) by using & = (m St.) evaluated the expectation parameters M Note that both a, and 9 (f:) depend on 8, but we have suppressed the subscript for notation simplicity.
Similarly to the VGP case, the X are the optimal dual parameters that measure the sensitivity of the solution to the perturbation in the moments of the posterior marginal (fi). This suggests that can design a similar VEM procedure that exploits the structure of solution in Eq. (20). The structure is shown below in terms of the natural parameterization of du (u) for sufficient statistics T(u).
(S;) m
K
and (S:) -1 - K mu + K. (E
Kmn KT. K. (21)
=A;
the quantities Kin and Kui directly depend on O and we can express the ELBO as the partition function of a Gaussian distribution, similarly to Eq. (17) (exact expression in App. B). For large data sets, storing all the [x;]% might be problematic, and we can instead store only Xi m-length vector, and As. a m X m matrix. This tied parameterization is motivated from the site-tying setting in sparse EP [1. 24] where the goal is to reduce the storage. The parameterization ignores the dependency of k. over f and may reduce the coupling between A and Zu but it is suitable for large data sets. An alternative tying method consists in storing the sums and the flanking K-1 terms.
We detail the final stochastic variational procedure which we refer to as t-SVGP Given a parameter 0., we run few iterations of the E-step. At each iteration & of the E-step, given m! and ((1), sample a minibatch M and compute the natural gradients by first computing
(k)
(k)
9u (1.) v pog fi)l
and 2((A
(e) (f.) -- logp(w fi)
using the marginals (k) (fi) from Eq, (6). The natural gradients of the expected log-likelihood for
the jh site is then equal to g. (k) (a m (R) (k). aa; 2(4)).
Using these natural gradients we can use an iterative procedure similar to Eq. (18) but now on the tied parameters.
+
X (k+1) <- (1-rx)ā (R) + Fk Σiem knigl.) (k)
(22)
A (++1) t- (1 (1-ra) (1-rx)Ag (K) + Tk Σiem TT (K)
(23)
7

---

Le(E (Bota).0) 9)
(E (Bora), 8)
Le(E (8ota), 8)
2
2
Rold
1.5
Oaid
1.5
0.dd
1.5
1.5
2.5
2,5
2
2
2
Pold
1.5
0oid
1.5
Roid
1.5
1
0.5
0.5
0.5
1.5
2.5
0.5
1.5
2.5
0.5
1.5
2.5
0
(a) 1-SVGP (ours)
(b) q-SVGP (whitened)
(c) q-SVGP
Figure EM iterations for hyperparameter learning on classification task on a a toy dataset. Top row: for the M-step, the optima of the t-SVGP objective (left) less sensitive to the initial hyperparameter value Oold. compared to q-SVG (middle and right). Bottom row: the EM iterations shown in white top of the EM objectives for range of starting values for Oold- The starting point 0 = 2.5 is marked with a white dot and the optimum with a star. t-SVGP converges much faster (2 iterations in leftmost plot) compared to 9-SVGP (middle and rightmost plots which take >5 iterations).
The natural parameter required can be obtained using Eq. (21),
S
K_1 m K- A K and m +- ((K) (
(24)
After few E-steps, we update the parameters with gradient descent step using the gradient of the log-pan partition function of Of (k) (u) with respect O. In App. we detail how to efficiently make predictions and compute the ELBO under parameterization Eq. (21). The full algorithm is given in
App. The convergence the sequence of stochastic updates the E-step is guaranteed under mild conditions discussed in [21] for the untied setting. Site-tying introduces bias but does not seem to affect convergence practice.
Empirical Evaluation
We conduct experiments to highlight the advantages of using the dual parameterization. Firstly, we study the effects of the improved objective for hyperparameter learning of t-SVGP versus Q-SVGP We study the objective being optimized for a single M-step, after E-step until convergence. We then show a full sequence of EM iterations on small data sets. For large-scale data, where running steps to convergence is expensive, bound and faster natural gradient computations show benefits in both settings. It is worth noting that the E-step for both VVP with natural gradients and t-SVGP are identical up to machine precision, and any differences in performance are to be attributed to the different parameterization.
partial E and M-steps and mini-batching. Our improved
The Role of the Learning Objective In Fig. 1, we learn the kernel hyperparameters 0 in GP classification task via coordinate ascent of the lower bound (E 0), where are the variational parameters, iee via EM. Starting at hyperparameter Oold. we denote by 8* (0.ia) the associated optimal

---

-30 200 200
TUN
CONCRETE
IUNOSPRERE
-001
DIABETES
-30
2200
500
1000 SBC
400
5-
101
15
20
10
10
16
12 IAE
NL.PO
211
Stops
Stepe
Steys
Stegs
Steps
9-SVGP (whitened)
v-SVGP (whitened, with natgrads) 1-VGP (with natgrads)
Figure 3: Comparison of convergence in terms of ELBO and negative log- seg-predictive density (NLPD) as averages over 5-fold cross- cross- validation runs, on UCI regression and classification tasks. All methods were trained (incl. hyperparameters and 50 inducing input locations) with matching learning rate. Natural gradient-based training is superior, and the I-SVGP parameterization improves stability.
variational parameters. Updating 0 consists in optimizing Ce (O 0) which we show on the left panel for the dual parameterization (blue), the standard whitened (orange) and unwhitened (green) SVGP parameterizations. The dual parameterization leads to a tighter bound and thus to bigger steps and faster overall convergence as shown on the right for the illustrative toy classification task, starting (0,, 02) (1, 1). in the extreme case of taking both the and M step to convergence. For the toy data set we use m 10 inducing points (see details in App. F).
In Fig. 2, also use the toy data and parameters in Fig. 1, but we show how the learning objective changes over iterations. The blue contours show, for all initial Ood.. the objective maximized in the M-step, i.e. (Gtaa( 0). The orange lines show, for all initial Oold. the outcome of an E-step followed by and M-step, i.e. 0' (f = arg max, (E (Oola, 0). The EM iterations converge to the fixed points of 0", i.e. its intersection with the diagonal line of the identity function. flatter line around the optimal value is more desirable as it means the iterations converge faster to the optimum value which this experiment is just below one, while line close to the diagonal leads to slow convergence. Here t-SGG has the fastest convergence, 9-SV performs poorly, although whitening clearly helps the optimisation problem. The dark dashed lines show how optimising 0 would look starting from Og = 2.5 and running 8 iterations for the different models.
Evaluation on UCI Classification and Regression Tasks We use common small and mid-sized UCI data sets to test the performance of our method against q-SVGP with natural gradient optimisation and normal g- SVP trained with Adam optimizer for the variational parameters. All methods use Adam for the hyperparameters. The exact details of the data sets can be found in App. F. Here we again take the approach that the optimal way to optimize the ELBO computational budget allows is to alternate between performing E and M-steps till convergence. We plot how the different inference schemes perform for ELBO and NLPD on a hold test set. We perform 5-fold cross validation with the results in Fig. 3 showing the mean of the folds for ELBO and NLPD. Natural gradient variants of q-S VVPG clearly perform better than non natural gradient q-S Our method 1-SVGP seems more stable specifically in NL.DD for most data sets if not equal to SVPG For gVSVG we have used the whitened version which as we noted helps with hyperparameter optimisation.
Improved Efficiency in Large-scale Inference To highlight practical benefits, we show the performance of our stochastic and sparse t-SVGP framework on the MNIST ([23), available under task (for details see App. F). Given the data set is n = 70,000,
CC BY-SAA 3.0) minibatching is needed and we adopt the parameterization of Eq. (21) meaning we match SVVG for parameter storage complexity. We compare against the natural gradient SVGP implementation but not the nan-natural gradient version since it produces considerably worse performance. In large scale minibatching experiments performing full E and M steps may not be efficient. Instead we perform partial steps for both. A single and M step can be thought of as the approach outlined in [33]. All experiments are performed with batch size of ni 200 and m = 100 inducing points and the

---

-10
T.SYPP 9-SVGP
-0.5
ELBO
-1.5
SVSP
0.5
SVGP
40
80
20
60
Walk-clock time (4)
Wall-clock time (s)
Figure 4: Comparison of practical inference and learning on MNIST. We compare training time laptop between t-SVGP to the q-SVGP model in GPilow in terms of wall-clock time of training for 150 steps, where both methods use natural gradient updates and share the same learning rates.
optimization is ran until convergence using the Adam optimizer for the hyperparameters (M-step). Table shows different variations of learning rates and iterations of E- and M-steps. The results suggest some benefits in running partial EM steps. The t-SVGP formulation performs equally not better than q-SVGP under all settings.
In Fig. 4, show the speed advantage of t-
Table I: NLPD on MNIST benchmarks for dif-
SVGP over q-SVGP due to cheaper natural gra- dient updates. We compare against the state-of- the-art implementation of SVGP in GPflow ([26). 12.2.1) and closely matched implementation of our method in GPflow. We compare wall-clock time compute 150 steps of the algorithm for both methods in terms of NLPD and ELBO taking single E and M-steps (MacBook pro, 2 GHz CPU, 16 GB RAM). Our implementation avoids the use of sluggish automatic differentiation to compute
ferent learning rates and E and M steps.
NLPD
LR
9-SVGP
1-SVGP
0.304:10.01
0.304 006
0.040
0.05
0.289:10.010
0.283:1.000
0.035
0.10
0.293-20.00
0.2814:0 010
0.030
0.259
0.255+0
0.025
0.03
0.282:10.00
0.283:10
0.243310.003
0.23010.000
0.030
the natural gradients, and, even if our implementation is not as optimized as SVGP in GPflow, is roughly 5 times faster on this standard benchmark.
6 Discussion and Conclusion
Sparse variational GP (SVGP) methods are the current de facto approach to allow GPs to scale to large problems. In this paper, we introduced an alternative parameterization to variational GPs that leads improved loss landscape for learning (cf., Fig. 1). This improvement hinges writing the variational problem in terms of its aual-simille to the conjugate-computatio variational inference (CVD) approach by Khan and Lin prr=-prrmmtterization approximate posterior decomposes into a prior contribution and a Gaussian approximate likelihood contribution. Variational inference under this model can conveniently be implemented by mirror descent and corresponds to natural gradient based learning. thus improving convergence in variational parameter optimization (the 'E-step', ), at the same time as improving hyperparameter optimization (the 'M-step'), due to the tighter evidence lower bound.
capture sites: assume the
We further show that we can derive the sparse equivalent of this method, which also allows for stochastic training through mini-batching, reducing the computational complexity to O(m² + mgm?? per step. Our method matches the asymptotic computational cost of other SVGP methods, while marginally reducing compute due to simpler expressions to paok-propagate through (see discussion in Sec. 5). Our empirical validation across wide variety of regression and classification tasks confirms the benefits suggested by our theory: The proposed strategy typically allows for improved stability gold-standard SVGP methods even when the learning rates remain the same. It allows for higher learning rates, and reduces computational cost leading to improved learning both in terms of reduced steps as well expected wall-clock time.
provide reference implementation of our method under the GPflow framework at https: // //gi thub.
10

---

Acknowledgments and Disclosure of Funding
AS acknowledges funding from the Academy of Finland (grant numbers 324345 and 339730). We acknowledge the computational resources provided by the Aalto Science-IT project. We thank Stefanos Eleftheriadis, Richard E. Turner, and Hugh Salimbeni for comments on the manuseript.
References
[1] T. D. Bui, J. Yan, and R. E. Turner. A unifying framework for Gaussian process pseudo-pos poouudto-ootttooo approxi- mations using power expectation propagation. Journal of Machine Learning Research, 18( 18(1):6449-3720, 2017.
[2] E. Challis and D. Barber. Concave Gaussian variational approximations for inference in large-scale Bayesian linear models. In Proceedings of the Fourteenth International Conference Artificial Intelligence and Statistics (AISTATS), volume Proceedings of Machine Learning Research, pages 199-207 PMLR. 2011.
[3] E. Challis and D. Barber. Gaussian Kullback-Leibler approximate inference. Journal of Muchine Learning Research, (4(32)2239-22335 2013.
[4] P. F. Chang, W. Wilkinson, M. F. Khan, and A. Solin. Fast variational learning state-space Gaussian process models. In 2020 JEEE 30th International Workshop on Machine Learning for Signal Processing (MLSP), pages 1-6. IEEE, 2020.
[5] C.-A. Cheng and B. Boots. Incremental variational sparse Gaussian process regression. In Advances in Neural Information Processing Systems 29 (NIPS), pages 4410-4418. Curran Associates, Inc., 2016.
[6] C.-A. Cheng and B. Boots. Variational inference for Gaussian process models with linear complexity, In Advances in Neural Information Processing Systems 30 (NIPS), pages 5184-5194, Curran Associates, Inc., 2017.
[7] L. Csató. Gaussian Processes: Iterative Sparse Appreximations. PhD thesis, Aston University, Birmingham, UK, 2002.
L. Csató and M. Opper. Sparse on-line Gaussian processes. Neural Computation, 14(3):641-663 2002.
[9] J. Gardner, G. Pleiss, K. Weinberger. D. Bindel, and A. G. Wilson. GPyTorch: Blackbox matrix matrix Gaussian process inference with GPU acceleration. In Advances Neural Information Processing Systems 3/ (NeurIPS), pages 7576-7586. Curran Associates, Inc., 2018.
[10] Hensman, Fusi, and N. D. Lawrence. Gaussian processes big data. Proceedings of the 29th Conference on Uncertainty in Artificial Intelligence (UAI), pages 282 290. AUAI Press, 2013.
[11] Hensman, A. Matthews, and Z. Ghahramani. Scalable variational Gaussian process classification. In Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics (AISTATS). volume 38 of Proceedings of Machine Learning Reseurch, pages 351-360. PMLR, 2015.
[12] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14(4): 1 303-1347, 2013.
[13] E. T. Jaynes
Information theory and statistical mechanics. Physical Review, 106:620-6333 630, 1957.
[14] E. T. Jaynes. On the rationale of maximum-entropy methods. Proceedings of the IEEF, 70K9):93-955,,
1982.
[15] M. Khan, D. Nielsen, V. Tangkaratt, W. Lin, Y. Gal, and A. Srivastava. Fast and scalable Bayesian deep learning by weight-perturbalte in Adam. In Proceedings of the 35th International Conference Machine Learning (ICML), volume 80 of Proceedings of Machine Learning Research, pages 2611-2620 PMLR, 2018
[16] M. E. Khan. Decoupled variational Gaussian inference. In Advances in Neural Information Processing Systems 27 (NIPS), pages 1547-1555. Curran Associates, Inc., 2014.
[17] E. Khan and W. Lin. in non non-conje models to inferences in conjugate models. In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics (AISTATS), volume 54 of Proceedings of Machine Learning Research, pages 878 887. PMLR, 2017.
variational inference: Converting variational inference

---

[18] M. E. Khan and D. Nielsen. Fast yet simple natural-gradien descent for variational inference in complex models. In 2018 International Symposium on Information Theory and Its Applications (ISITA), pages 31-35. IEEE, 2018.
[19] E. Khan and Rue. Learning- algoritths from Bayesian principles. arXiv preprint arXiv, 2021.
[20] M. E. Khan, A. Aravkin, M. Friedlander, and M. Seeger. Fast dual variational inference for non-conjugate latent Gaussian models. In Proceedings of the 30th International Conference on Machine Learning (ICML). volume 28 of Proceedings of Machine Learning Research, pages 951-959. PMLR, 2013.
[21] M. E. Khan, R. Babanezhad, W. Lin, M. Schmidt, and M. Sugiyama. Faster stochastic variational inference using mexxiallggadienn methods with general divergence functions. In Proceedings of the Thirty- Seccoy-Second Conference on Uncertainty in Artificial Intelligence (UAT), pages 319-328. AUAI Press, 2016.
[22] G. Kimeldorf and G. Wahba. Some results on Tchebycheffian spline functions. Journal of Mathematical Analysis and Applications, 33(1):82-95, 1971.
[23] Y. LeCun, C. Cortes, and C. J. Burges. The MNIST database of handwritten digits, 1998. URL http: //yann. .lecun. com/exdb/mniste
[24] Y. Li, J. M. Hernández-Lobato, and E. Turner. Stochastic expectation propagation. In Advances in Neural Information Processing Systems 28 (NIPS), pages 2323-2331. Curran Associates, Inc., 2015.
[25] L.. Malagò and G. Pistone. Information geometry of the Gaussian distribution in view stochastic optimization. In Proceedings of the 2015 ACM Conference on Foundations of Genetic Algorithms XIII. pages 150-162, 2015.
[26] A. G. d. G. Matthews, M. van der Wilk, T. Nickson, K. Fujii, A. Boukouvalas, P. León- Villagrá, Z. Ghahra- mani, and J. Hensman. GPfoow: A Gaussian process library using TensorFlow. Journal of Machine Learning Research, 18(40):1-6, 2017.
[27] T. P. Minka. Expectation propagation for approximate Bayesian inference. In Proceedingy of the Seven- teenth Conference on Uncertainty Artificial Intelligence (UAI), volume 17, pages 362-369. AUAI Press. 2001.
[28] H. Nickisch and C. E. Rasmussen. Approximations for binary Gaussian process classification, Journal of Machine Learning Research, 90ct1):035-20778 2008.
[29] M. Opper and C. Archambeau. The variational Gaussian approximation revisited. Neural Computation, 21 (3):786-792, 2009.
[30] J. Quilooneoo-anneee Candela and C. E. Rasmussen. A unifying view of sparse approximate Gaussian process regression, Journal of Machine Learning Research, 6(Dec):1939 19599,999.01999 2005.
[31] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning MIT Press, Cambridge, MA, 2006.
[32] H. Salimheni, C.-A. Cheng, B. Boots, and M. Deisenroth. Orthogonally decoupled variational Gaussian processes, Advances in Neural Information Processing Systems 31 (NeurIPS), pages 8711-8720. Curran Associates, Inc., 2018.
[33] H. Salimbeni, S. Eleftheriadis, and J. Hensman. Natural gradients in practice: Non-conjugate variational inference in Gaussian process models. In Proceedings of the Twenty-Fine First International Conference on Artificial Intelligence and Statistics (AISTATS), volume 84 of Proceedings of Machine Learning Research, pages 689-697. PMLR, 2018
[34] M. Seeger. Bayesian Gaussian Process Models: PAC-Bayesian Generalisation Error Bounds and Sparse Approximations. PhD thesis, University of Edinburgh, Edinburgh, UK, 2003.
[35] Titsias. Variational learning of inducing variables sparse Gaussian processes. In Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics (AISTATS), volume of Proceedings of Machine Learning Research, pages 567-574. PMIR, 2009.
[36] M. van der Wilk, V. Dutordoir, S. John, A. Artemev, V. Adam, and J. Hensman. A framework for interdomain and multioutput Gaussian processes. arXiv preprint 2003:2003300111 2020.
[37] M. J. Wainwright, M. I. Jordan, et al. Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning. 1(1-2):1-305, 2008.
12

---

[38] W. Wilkinson, A. Solin, and V. Adam. Sparse algorithms for Markovian Gaussian processes. In Proceedings of The 24th International Conference Artificial Intelligence and Statistics (AISTATS), volume 130 of Proceedings of Machine Learning Research, pages 1747-1755. PMLR, 2021.
[39] C. K. Williams and M. Seeger. Using the Nyström method to speed up kernel machines. In Advances in Neural Information Processing Systems 13 (NIPS). pages 682-688. MIT Press, 2001.
13


Reviews:
[{"review_id": "xzbSJFmt-MN", "reviewer": "Reviewer_zj8L", "summary": "This paper proposes to speed-up approximate inference in the sparse variational GP model using a dual parameterisation (similar to an EP approach) where each data point is assigned an additional dual variable that may be updated efficiently using natural gradient descent with reduced complexity on the natural gradient updates that also lead to a tighter bound on the evidence. The paper illustrates the efficacy of the approach over standard SVGP variational inference on a number of standard benchmark datasets and demonstrates good empirical performance whilst reducing training time.", "questions": "", "limitations": "", "rating": 8, "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": "I found the paper to be clear to read with a good coverage of the background material and justification for the proposed methodology. The authors highlight that the proposal brings together a number of approaches that are already established (e.g. SVGP, natural gradient descent and EP) but, in my opinion, this is not a limitation as the paper treats the material well, does not overclaim the novelty but instead makes the point that bringing these things together in this way yields a very satisfactory approach that is neither well-known nor well used (to the best of my knowledge). That it performs competitively in established code-bases (such as GPflow) is a good indication of importance of the result and its general applicability. In particular, the recognition that improving the calculations of the natural gradient updates is probably the most satisfying part of the paper and the case for this is well made.\n\nOverall, I believe that the paper provides sufficient empirical evidence to support the theoretical ideas behind the proposed method. I can see that there might be objections to datasets such as MNIST but I agree that it is the scaling we are interested in here and the separation of the objective wrt hyperparameters and variational params provides a strong motivation and good result for the paper. \n\nTo the best of my knowledge the approach described has not been explicitly presented before and I am in favour of accepting the paper.\n\nOther comments:\n\n- Thinking about fair comparison for the best q-SVGP - a number of papers have questioned whether the alternation approach of EM does always lead to improved convergence (i.e. is coordinate ascent always better than something like a Gauss-Newton approach that approximates the off-diagonal blocks of the Hessian somehow). I wonder if the authors could comment upon whether or not they have tried these sorts of approaches? I can well believe that the natural gradients lead to better conditioning (effectively the appropriate blocks of the Hessian for approx second-order optimisation) but, if I understand correctly, they will assume “zero blocks” in the Hessian in the cross terms between, say, the natural gradient parameters and the hyper parameters. My understanding is that these are usually dealt with first order approaches with momentum, etc.. but I don’t know if there has been any more in-depth investigation for this approach? E.g. something like L-BFGS or similar?\n\n- Could the authors avoid using “e” as a parameter (e.g. (11)) since its also used as the exponential (e.g. (12))?\n\n-  Typos: lines 81 (-them by-), 91 (_integrals_), 110 (Gradient Descent?)", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "No ethical concerns.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "v-b5LEKC81A", "reviewer": "Reviewer_SLDR", "summary": "This paper proposed a new dual parametrization method to improve the optimization efficiency for variational sparse Gaussian processes. The method speeds up inference than natural gradient descent and achieves a tighter lower bound for hyperparameter learning. ", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "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: This paper introduces a new method for optimizing the variational sparse Gaussian processes.  \n\nQuality: The authors should focus more on the technical writing to make them error-free:\nLine 68: makes it a drawback .....\nLine 130: can be approximated....\nLine 519 (appendix)\n......\nClarity: As mentioned above, the technical writing is not very sound and fluent which is misleading sometimes.\n\nSignificance: The experimental results generally demonstrated the technical part in the paper. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Focus more on the writing", "review_text": "", "extra_scores": {"ethical_concerns": "N.A. ", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "tO9Ml696Gmh", "reviewer": "Reviewer_aCDy", "summary": "The paper revisits a known, but seemingly forgotten dual parametrization for multivariate Gaussian distributions that naturally lends itself to variational inference in Gaussian process models. This dual parametrization was introduced in previous works for variational Gaussian inference using a Lagrangian relaxation of the ELBO loss. Similarly to commonly used parametrizations, the dual parametrization can be combined with natural gradients based optimization that allows for faster training of the variational distribution by exploiting the information geometry of distribution space. However, the required natural gradients are considerably cheaper to compute in this setting due to only requiring first-order gradients. This parametrization is then incorporated into a sparse variational inference scheme for GPs. One of the computational difficulties when training the the SVGP lies with the coupled optimization of the variational distribution and the prior hyperparameters, for which several heuristics have been introduced in the community. The experiments compare the proposed parametrization for SVGP with the commonly used whitened one. They demonstrate that for small and mid-sized datasets, gradient descent along the trajectory of optimal variational distributions results in more stable learning of the hyperparameters. For larger datasets, the coupled stochastic optimization of the two (with only partial optimization of the variational distribution) leads the dual formulation to equivalent or slightly better optima measured in terms of predictive density, but it is considerably faster than the baseline due to cheaper natural gradients. ", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "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": "This is a well-written paper that recycles a previously known parametrization for variational Gaussian inference that has not been employed in the SVGP context before. The main idea is therefore relatively simple, which in concept follows from simply putting these together. However, in my view, there is considerable value in this due to the nontrivial insight that is required to carry out the required calculations, put them into an implementation and evaluate them experimentally. Especially so, because the main problem that the paper is trying to address (optimization in SVGP) is of considerable significance, and the results are of interest. Judging purely by the results, the approach does seem promising and I think this paper can be a valuable contribution to the GP community.\n\nOne obstacle I see to these ideas becoming adopted and their use more wide-spread is that the description of the proposed approach is not very clear at the moment. In particular, Section 3 is hard to follow because no intuition is provided about the dual parametrization, and only a brief survey is given (although with several references), which is not an adequate initiation to the approach for those not already familiar with it. Given that the main tenet of the paper is that this parametrization seems to be unexplored in the SVGP context, I would have appreciated if a more detailed description was given. In contrast, Section 2 spends about two pages on recapping well-known ideas. Some of this space could probably be reallocated to the description in Section 3. This also makes Section 4 somewhat hard to follow, where the idea is combined with the SVGP. Section 5 on the other hand is very clear, the experimental evaluation and the results are interesting.\n\nOverall, although neither of the core ideas in the paper are novel, I believe that a novel combination of known ideas can also be a valuable contribution. However, my main criticism of the paper is regarding the presentation of the idea that they are trying to introduce, which in its current form might limit their use due to lack of clarity. If the authors can improve on the presentation in these sections, I would be inclined to slightly raise my rating.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "At the moment there is no discussion of limitations, which may be useful to include with the experimental evaluation. In particular, if there are any drawbacks of the proposed parametrization compared to the standard one would be useful to know. It is also not clear to me if the discussed dual parametrization limits the space of approximating distributions, or if any multivariate Gaussian can be represented in this form as well. Also, did the authors face any challenges during the implementation of the proposed model, or during its optimization?", "review_text": "", "extra_scores": {"ethical_concerns": "No concerns.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "p0zaV3m9rL", "reviewer": "Reviewer_r5PB", "summary": "The authors argue for a different (existing) parametrization for training (non-conjugate) Gaussian Processes.\nThey show that this parametrization can be interpreted as the optimal result of a dual problem, and that one can compute easily natural gradient updates.\nAdditionally this parametrization \"may\" lead to a tighter ELBO because of a stronger dependency on the hyperparameters (kernel parameters).\nExperiments are run to show the improvement over the standard \"q-SVGP\".", "questions": "", "limitations": "", "rating": 6, "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 paper is well written and exposes clearly all the different elements to understand the theory and method introduced.\nThe main contributions are to bring a different view to the parametrization introduced in Csato 02' as an optimization problem that is mathematically interesting and derive natural gradient updates for this new parametrization.\nSince the parametrization is already well known, the contribution is only limited and some of the mentioned advantages have flaws.\n\n- There is a large flaw in the reasoning and the comparison with Salimbeni 18'. It is argued that  \"[t-SVGP] is much cheaper than previous calculations of the natural gradient\", because one supposedly does not have so many reparametrizations to work with. However $\\nabla_\\mu E_{q(f)}[\\log p(y|f)]$ requires the same parametrization as Salimbeni 18' since gradients given the second expected parameter $(mm^\\top + S)$ cannot be found directly (and therefore require a reparametrization as well)\n-  The mention \"It is a first order update\" (L175)  sounds wrong since gradient derivatives given $S$ imply second order derivatives.\n- Using natural gradients in Salimbeni 18' requires extra-care with the learning rate, especially in the first iterations. Some additional comments in this regard for t-SVGP would be good.\n- It is not enough in the experiments to compare only to the standard q-SVGP, there should be at least a direct comparison with Opper 09', so one can see what the natural gradient approach actually brings.\n- There are missing references for relevant work on this area, notably via the use of latent variable augmentation, see e.g. Wenzel 19' or more generally Galy-Fajou 20'\n- Minor comments:\n\t- L100: \"leads to a convex objective\" Since the objective is defined to be the ELBO, it should be concave.\n\t- L192: \"after one E-step\" Repetition with the beginning of the sentence\n- References\n\t- Csato 02' : Sparse on-line Gaussian processes, 2002, Csato and Opper\n\t- Opper 09' : The Variational Gaussian Approximation Revisited, 2009, Opper and Archambeau\n\t- Salimbeni 18' : Natural Gradients in Practice: Non-Conjugate Variational Inference in Gaussian Process Models, 2018, \t\tSalimbeni, Eleftheriom and Hensman\n\t- Wenzel 19' : Efficient Gaussian Process Classification Using Polya-Gamma Data Augmentation, 2019, Florian Wenzel, Theo Galy-Fajou, Christan Donner, Marius Kloft and Manfred Opper\n\t- Galy-Fajou 20' : Automated Augmented Conjugate Inference for Non-conjugate Gaussian Process Models, 2020, Théo Galy-Fajou, Florian Wenzel, Manfred Opper", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Authors mention limitations in Section 6, but there is no mention on how the method could fail.", "review_text": "", "extra_scores": {"ethical_concerns": "Not applicable", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "HGgg5I-NFmK", "reviewer": "Reviewer_XG5T", "summary": "The authors introduce an efficient natural gradient based approximate inference scheme for sparse variational Gaussian Process models. The approach leverages a \"dual parameterization\" that exploits Lagrange multipliers along the lines of \"Decoupled variational Gaussian inference\" (Khan 2014). These dual parameters are conceptually similar to site parameters in expectation propagation. In experiments it is shown that this dual parameterization scheme results in faster optimization and tighter bounds (although solid theory is lacking for the latter claim, at least for the non-conjugate case).", "questions": "", "limitations": "", "rating": 6, "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": "Gaussian process (GP) models are popular due to their good uncertainty quantification, interpretability, and ability to include prior information through kernel design. Thus improvements to GP inference are likely to be of interest to an active if somewhat moderately sized NeurIPS subcommunity. This submission has pretty clear exposition and does a reasonably good job of positioning itself w.r.t. prior work. The experiments are moderately comprehensive and give some amount of confidence in the general claim that this method offers non-negligible inference speed-ups. The work is somewhat incremental, as it directly builds on other references, including e.g. various Khan et. al. and Hensman et. al. papers cited in the submission. Nevertheless, I think this work could be of interest to the NeurIPS community, since it offers concrete speed-ups over natural gradient approaches to SVGP as well as an overview/application of approaches like \"Decoupled variational Gaussian inference\" that, as the authors say, can be described as a \"known, yet rarely used parameterization for Gaussian distributions.\"\n\nGiven the somewhat incremental nature of this work, I think the work could be considerably improved if the empirical evidence for inference speed-ups could be expanded. In particular as far as I can tell all the experiments use a pretty modest number of inducing points (100). Does the quoted 5x speed-up (line 322) persist for 500 or 1000 inducing points? Does the performance improvement disappear or largely disappear when implemented on GPU? As far as I can tell the experiments are all performed on CPU.\n\nAdditional comments/questions:\n- “The formulation is fast since it avoids the use of sluggish automatic differentiation to compute the \nnatural gradients. “ (line 42) Can you clarify why this is exactly? As I understand it we are talking about a multiplicative constant and not an asymptotic speed-up. Has this more to do with the limitations/flexibility of current AD systems or is this something more fundamental?\n- Typos: mirrordescent (x2), optimial\n\n## After author response\n\nI thank the authors for their response. Given the (somewhat) incremental nature of the work I maintain my original score. Nevertheless I would be in favor of seeing this submission accepted for publication, as I believe it could be of interest to the community. I suggest the authors include experiments with a larger number of inducing points. In most practical applications (unless one is particularly compute constrained) M ~250-500+ inducing points often deliver tangible performance improvements over much smaller numbers of inducing points. As such this represents a common and important regime.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors do not explicitly discuss potential negative societal impact of their work.\n\nAs suggested above, I think it would be valuable if the authors further delineated regimes (and verified empirically) in which they do and do not expect substantial speed-ups. In particular CPU versus GPU and few inducing points versus many inducing points.\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1.2"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper revisits a known alternative, so-called dual, parameterization for models with Gaussian priors and iid likelihoods and applies it to stochastic variational inference in Gaussian process models. It shows _empirically_ that this parameterization results in faster optimization and tighter bounds that improve hyper-parameter learning. Results are shown on small UCI datasets and on MINST (n=70,000 datapoints).

Overall, the reviewers believe that the proposed approach is incremental as it takes a known idea (that of the dual parameterization) and applies to the SVGP context. However, all the reviewers agree that the paper provides a significant contribution to the NeurIPS community and that NeurIPS will benefit from knowing the details and findings in this paper. 

Besides its incremental nature, one of the major criticisms of this paper was with respect to its clarity (Reviewer aCDy), which the authors seemed to have addressed satisfactorily. Another important point (raised by reviewer XG5T) concerns whether the benefits of the approach would extend to larger problems for which one would need a higher number of inducing variables ($m>100$). The authors have clarified that the bigger the $m$ the smaller the marginal gain obtained by the proposed method. However, the proposed approach still provides benefits for settings of practical interest. 

Decision: Accept (Poster)