Spatio-Temporal Variational Gaussian Processes
Authors: Oliver Hamelijnck
William J. Wilkinson
Niki Andreas Loppi
Arno Solin
Theo Damoulas
Venue: NeurIPS 2021 Poster (2021)
Keywords: Gaussian processes
Variational Inference
Spatio-Temporal Analysis

Abstract:   We introduce a scalable approach to Gaussian process inference that combines spatio-temporal filtering with natural gradient variational inference, resulting in a non-conjugate GP method for multivariate data that scales linearly with respect to time. Our natural gradient approach enables application of parallel filtering and smoothing, further reducing the temporal span complexity to be logarithmic in the number of time steps. We derive a sparse approximation that constructs a state-space model over a reduced set of spatial inducing points, and show that for separable Markov kernels the full and sparse cases exactly recover the standard variational GP, whilst exhibiting favourable computational properties. To further improve the spatial scaling we propose a mean-field assumption of independence between spatial locations which, when coupled with sparsity and parallelisation, leads to an efficient and accurate method for large spatio-temporal problems.

Spatio-Temporal Variational Gaussian Processes
Oliver Hamelijnck'
William J. Wilkinson" Aalto University william. wi wilkiinson@alllo....
Niki A. Loppi NVIDIA Aloppi@nvidia.com
The Alan Turing Institute University of Warwick phameli.jcckkturing ac.uk uk
Arno Solin
Theodoros Damoulas
Aalto University
The Alan Turing Institute / University of Warwick
arno. arno.soling iollinnaalt....
tdamoulas@turing, ac. uk
Abstract
We introduce scalable approach to Gaussian process inference that combines spatio-temporal filtering with natural gradient variational inference, resulting in a non-conjugate GP method for multivariate data that scales linearly with respect to time. Our natural gradient approach enables application of parallel filtering and smoothing, further reducing the temporal span complexity be logarithmic in the number of time steps. We derive a sparse approximation that constructs a state-space space model over a reduced set spatial inducing points, and show that for separable Markov kernels the full and sparse cases exactly recover the standard variational GP, whilst exhibiting favourable computational properties. To further improve the spatial scaling propose meennfiield assumption of independence between spatial locations which, when coupled with sparsity and parallelisation, leads to an efficient and accurate method for large spatio-temporal problems.
Introduction
Most real-world processes occur across space and time, exhibit complex dependencies, and are observed through noisy irregular samples. Take, for example, the task of modelling air pollution across city. Such a task involves large amounts of noisy, partially-observe effects governed by weather, traffic, human movement, etc. This setting motivates probabilistic perspective, allowing for the incorporation of prior knowledge and the quantification of uncertainty.
data with strong seasonal
Gaussian processes (GPs, [38]) provide such probabilistic modelling paradigm, but their inherent cubic computational scaling in the number of data, N, limits their applicability to spatio-temporal tasks. Arguably the most successful methods address this issue are sparse GPs [37], which summarise the true GP posterior through reduced set of M inducing points and have dominant computational scaling M²), and spatio-temporal GPs [43], which rewrite the GP prior state-space model and use filtering to perform inference in O (Nd?), where d is the dimensionality of the state-space. Sparse GPs and apatio-tempone GPs have been combined by constructing Markovian system in which a set of spatial inducing points are tracked over time [24, 51].
the prior conditional model [24] or do not exploit natural gradients [45], meaning they do not provide the same inference and learning results state-of-the-art variational GPs [26] in the presence of non-conjugate like- lihoods or sparsity, which has hindered their widespread adoption. We introduce spatio-temporal
However, existing methods for spatio-temporal GPs make approximations
"quul contribution
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

SVGP
ST-SVGP recovers same spatio-temporal model as SVGP if the inducing points repeat over time
ST-SVGP
ST-SVGP's effective inducing point count scales with the data, but its computationnl scaling lincar in the number of time stops
Laducing points are sharee over the tempural
20101000
Figure 1: A demonstration of the spatio-temporal sparse variational GP (ST-SVGP) applied to crime count data in New York. ST-SVGI tracks spatial points over time via spatio-temporal filtering. The colourmap is the posterior mean, and the red dots are spatial inducing points. The diagram shows the difference between how inducing points are treated in ST-SVGP and SVGP.
variational GPs (ST-VGI) which provide the exact same results as standard variational GPs, whilst reducing the computational scaling in the temporal dimension from cubic to linear. ST-VGP is derived using a natural gradient variational inference approach based on filtering and smoothing. We also derive this method's sparse variant, and demonstrate how it enables the use of significantly more inducing points than the standard approach, leading to improved predictive performance.
We then show how the spatio-temporal structure can be exploited even further to improve both the temporal and spatial scaling. We demonstrate for the first time how to apply parallel filtering and smoothing [41] to non-conjugate GPs to reduce the temporal (span) complexity to be logarithmic. We then reformulate the model enable an efficient mean-field approximation across space, improving the complexity with respect to the number of spatial points. We analyse the practical performance and scalability our proposed methods, demonstrating how they make possible to apply GPs to large-scale spatio-temporal scenarios without sacrificing inference quality.
1.1 Related Work
modelling in both machine learning and spatial statisties [38, 22, 14, 6]. Many approaches to overcome their computational burden have been proposed, from nearest neighbours [17] to parallel algorithms on GPUs [48]. Within machine learning, the sparse GP approach is perhaps the most popular [37, 46), and is typically combined with mini-batching to allow training on massive datasets [26). However, fails in practical cases where the number of inducing points must grow with the size of the data, such as for time series [49].
GPs are commonly used for
When the data lie on a grid, separable kernels exhibit Kronecker structure which can be exploited for efficient inference [39]. This approach has been generalised to the partial grid setting [53]. and to structured kernel interpolation (SKI, [52]) which requires only that inducing points be on grid. Generally, these approaches are limited to the conjugate case, although Laplace- ppatio-tempora modelling that performs an orthogonal projection of the data to enforce independence between the latent processes.
extensions
exist [19]. Bruinsma et al. [10] present an approach
It has been shown that variational GPs be computed in linear time either by exploiting sparse precision structure [18] or via filtering and smoothing [11]. I]. Other inference schemes such as Laplace and expectation propagation have also been proposed [35, 51]. In the spatio-temporal case, sparsity has been used in the spatial dimension [24, 43]. These methods historically suffered from the fact that i) filtering was not amenable to fast automatic differentiation due to its recursive nature, and ii) state-of-the- the inference schemes had not been developed to make them directly comparable to other methods. The first is no longer an issue since many machine learning frameworks are now capable of efficiently differentiating recursive models [11]. We address the second point with this paper. A similar algorithm to ours that also sparse in the temporal dimension has been developed [50, 2]. and relevant properties of the spatio-temporal model presented here are also analysed in [45]. Fourier features [28] are an alternative approach to scalable GPs, but are not suited very long time series with high variability due to the need for impractically many inducing features.

---

2 Background
We consider data lying on a spatio-temporal grid comprising input-output pairs, (x(xt) RNxxNNDD, Y(st) E R Nx where N₂ is the number of temporal points, Ng the number of spatial points, and D 1 + Dg the input dimensionality (with DN being the number of spatial dimensions). We use and s to represent time and space respectively. The assumption of the grid structure is relaxed via the introduction of sparse methods in Sec. 3, and by the natural handling of missing data.
For consistency with the GP literature we let X vec(X("*)) E ER Y -ec(Y((st) E R where N = NrN, is the total number of data points. We use the operator vec(.) to simply convert data from spatio-tempone then space. For notational convenience we define Xn, = n.k+ Yn,k = y which indexes data at time index n and spatial index k. We use in to denote the n 'tth time point, S E NDN. to denote all spatial grid points and S₁ the k'th one. Let f .R. + R to be a random function with a zero-mean GP prior, then for given likelihood p(Y (X)) the generative model is,
grid into vector form, whilst keeping observations ordered by time and (st)
n.K
f(z) ~ GP(O, (P(,,,(Z,,))) x(x,2')), Y TI P(Y n.k fn,k); where fr, = f(X..)), and we let f. be the function values of all spatial points at time t.. When the kernel k is evaluated at given inputs we write the corresponding gram matrix as Kxx- = x(X, X'). To make explicit that takes spatio-temporal inputs we also abuse the notation slightly write f(z) - [tt,,) () and k(z, x²) k(t, s, f', s'). A summary of all notation used is provided in App. A. For Gaussian likelihoods the posterior, p(f Y). is available in closed form, otherwise approximations must be used. In either case, inference typically comes at cubic cost of O(N)
(1)
2.1 State Space Spatio-Temports Gaussian Processes
One method for handling the cubic scaling of GPs is to reformulate the prior in Eq. (1) as a state space model, reducing the computational scaling to linear in the number of time points [43]. The enabling assumption is that the kernel is both Markovian and separable between time and space: s(t,s,t', t's'''ss) = Hig (t, t') Kg ss, We use the term Markovian kernel to refer to kernel which can be re-written in state-space form (see [44] for overview). First, we write down the GP prior as stochastic partial differential equation (SPDE, see [15)) (t, = Ag ft, + L. w(t, s). where w(1, s) is a (spatio-temporal) white noise process and A, a suitable (pseudo- snreudd-xdifferentia )differential operator Isee 42]. By appropriately defining the model matrices and the white noise spectral density function, SPDEs of this form can represent a large class of separable and non-separable GP models.
When the kernel is separable, this SPDE can be simplified to a marginalising to a finite set of spatial locations, S E IR. N, giving, (f) = F f(t) dt L (t),). where f(t) is the Gaussian distributed state at the spatial points S at time t, with dimensionality d Nodi, where de is the dimensionality of the state- state-space model induced by Kg <e( ). dd(tt has spectral density Qc. and the matrix H extracts the function value from the state: f, = Ht.). F and L are the feedback and noise effect matrices [42]. This simplification to SDE is possible due to the independence between spatial points at time t and all other time steps. given the current state [45]. This follows from the fact that for any separable kernel, [t, s) and f(l', s') are independent given f(t', s) [36]. For a step size An = lnt1 - in. the discrete-time model matrices are,
SDE [24] by
An = (FFA,,
Q, - ((A,. - 7) L Q. LT &(An ( dr,
where (--) is the matrix exponential. The resulting discrete model is,
f(tn+1) A, f(t..) + qr: where In ~ N(O, Q.)) If p(Y,, (t,)) is Gaussian then Kalman smoothing algorithms can be employed to perform inference in Eq. (3) in O(Nid)) O(N.Ngd)).
Y, [(tr,) ~ p(Y,, Hf((),))
(3)
Markovian GPs with Spatial Sparsity Sparse GPs re-define the GP prior over a smaller set of M inducing points: let u = fZ) E IR g be the inducing variables at inducing locations Z E RMX then the augmented prior p(f. u) = p(f u)p(u), where p(u) = N(u 0. Kzz), and with Gaussian conditional p(f u). If the inducing points are placed on a spatio-temporal grid, with Z, E RR M,xD, being the spatial inducing locations, the conditional p(f u) can be simplified to (see App. D):
uf N(f [I & (K
(9) SS X
K (e) ]
(t) u, Lt (9 Qs).
(4)
3

---

where Q3 K K dent training conditional (FITC) method [37] approximates the full conditional covariance with its diagonal, leading to the following convenient property: gFire ()) = II QFITC |u) ss X K z. is block diagonal. This factorisation across time allows the model to be cast into the state-space form of Eq. (3), but where the state f(t) defined over the reduced set of spatial inducing points [24]. Inference can be performed in (N,M3&
Zg Z.S (see App. A for notational details). The fully indepen-
K(b)
IhT Qfrrc(ff, 1,,), where the last equality holds because I % (Kss
2.2 Sparse Variational GPs
To perform approximate inference the presence of sparsity or non- non-Gaussian likelihoods, variational methods cast inference as optimisation through minimisation of the Kullback-Leiblee divergence (KLD) from the true posterior the approximate posterior [8]. Although direct computation of the
KLD is intractable, it can be rewritten as the maximisation of the evidence lower bound (ELBO).
Unlike FITC, the sparse variational GP VSGG,, [46]) does not approximate the conditional p(f u). but instead approximates the posterior as q(f, u) = p(f u) g(u), where (uu) = N(u m, P) is Gaussian whose parameters are to be optimised. The SVGP ELBO is:
LSVGP (gg(u) [Ey(e] u) [logp(X Y []]
KL (g p[)),
(5)
which can be computed in O(N + M SVGP has many benefits over methods such as FITC, including: Gau-Saussian likelihoods can be handled through quadrature or Monte-Carlo approximations [27, 33], it applicable to big data through stochastic VI and mini- batihbg [26],
and the inducing locations are 'variationally protected' and hence prevent overfitting [7].
Natural Gradients Natural gradient descent calculates gradients in distribution space rather than parameter space, and has been shown to improve inference time and quality for variational GPs [26, 40]. A natural descent direction is obtained by scaling the standard gradient by the inverse of the Fisher information matrix, (-g()) [ log 2()] [5]. For a Gaussian approximate posterior, the natural gradient of target with respect the natural parameters À can be calculated without directly forming the Hessian, since it can be shown to be equivalent to the gradient with respect to the mean
parameters f [m, mm + P] [25]. The natural parameter update, with learning rate 33, becomes,
OC
x x 9
0μ
To update the approximation posterior, À can be simply transformed to the moment parameterisation [m, PJ. A table of mappings between the various parametrisations is given in App. G.
CVI and the Approximate Likelihood Khan and Lin [31] show that when the prior and approx- imate posterior are conjugate (as in the GP case), further elegant properties of exponential family distributions mean that Eq. (6) is equivalent to a two step Bayesian update:
i 3 - (1 - B) X x (Egr)]log p( Y [f)
A 4 n
(7)
where 77 are the natural parameters of the prior and À are the natural parameters of the likelihood contribution. Letting g(.) 82, po f)] the gradients can be computed in terms of the mean and covariance via the chain rule: g(u) = g (m) 2 gg)) m, g(P)). Eq. (7) shows that, since the prior parameters n are known, natural gradient variational inference completely characterised by updates to an approximate likelihood, which we denote N(Y f, V). parameterised by covariance V ==6 (-2x())) and mean Y = ( (see App. A). The approximate posterior then has the form,
a.
N(Y f. P() (f) J (. df
(8)
Computing (f) amounts to performing conjugate GP regression with the model prior and the approximation likelihood. This approach is called menjugateecomputation variational inference (CVI, [31). To re-emphasise that the CVI updates compute the exact same quantity as Eq. (6), we provide an alternative derivation in App. H by directly applying the chain rule to Eq. (6).

---

3 Spatio-Temporal Variational Gaussian Processes
In this section introduce a spatio-temporal VGP that has linear complexity with respect time whilst obtaining the identical variational posterior as the standard VGP. We will then go on to derive this method's sparse variant, which gives the same posterior SVGP when the inducing points set similarly (i.e., on a spatio-temporal grid), but is capable of scaling to much larger values M.
3.1 The Spatio-Tempores VGP ELBO
We first derive our proposed spatio-temporal VGP ELBO. We do this by exploiting the form of the approximate posterior after natural gradient step in order to write the ELBO as a sum of three terms, each of which can be efficiently computed through filtering and smoothing. As shown in Sec. 2.2, after a natural gradient step, the approximate posterior q(f) N(Y f. v)p(f) decomposes Bayesian update applied to the model prior with an approximate likelihood. Following Chang et al. [11] substitute Eq. (8) into the VGP ELBO:
log P (()) (f))
log P p(F) / N Y f, V) p(F) df
LVGP = (g(())
= Eg(()
(9)
N(Y V) pff
No
Byg(r...)) x)) IYg n. [[n,k]] ByyC)) log N(Y N(Y [f, V)] + log p(r) [N [f, V)]-
The first term is the expected log likelihood, the second is the expected log approximate likeli- hood, and the final term is the log marginal likelihood of the approximation posterior, log (Y) log Ey(r) [N(Y f, V) Naïvely evaluating LVGP requires (N computation for both the update to (f)) and the marginal likelihood. We now show how to compute this with linear scaling in No.
We observe that after natural gradient update, V, the approximate likelihood covariance, has the same form the gradient g(P) because, seen Eq. (7), À only additively updated by g(u). Since the expected likelihood, (y(() [logp(Y f)], factorises across observations, g(P) is diagonal
and hence is V. The approximate likelihood therefore factorises in the same way as the true one:
N.
log N Y (
Σ
log N( n,k [fnkk: Vn,k)-
(10)
We now turn attention to computing the posterior and the marginal likelihood. Recall that if the kernel is separable between time and space, k(t, s, s'', = (t, t') (s, then the GP prior can be re- ce-written as Eq. (3). This separability property further results in the state-space model matrices having a convenient Kronecker structure,
(11) where qn N(0, K/ & Q!")) and H IN, & H Here A E R.dedd Q0 R4 and (()) 1x are the transition matrix, process noise covariance, and measurement model of the SDE (see Sec. 2.1) induced by the kernel k,(-,)) respectively.
F(ta+1) IN. 2 A00) f(t...) + 9n
Y, f((n)) p( (f((),)),
Because the GP prior is Markov and the approximate likelihood factorises across time, the approx- imate GP posterior is also Markov [45]. Hence marginals q(f,) can be computed through linear filtering and smoothing applied to Eq. (11). Furthermore, the marginal likelihood of linear Gaussian state-space model, p pYY p(Y,) IIN- P(Y,, Y I1), can be computed sequentially by run- ning the forward filter, since p( Y, Y 11m-1)= = S n. Hf,)t)) P(((tt) Y Yim-1) (f where p(f(tn)) Y Y1m-11 is the predictive filtering distribution. By combining of the above properties we can now write the ELBO as
No N. Egie, Mogg(Ynn. [x.))]
N₂ N
LST-VGP =
Σ
Σ
EygC,
[log N(Y fn,ko Vn.k)]
n=1A=1
log E p(f(t.) R1mn-1) [Y.Y. Hf,t,), V.)]..
(12)
5

---

Algorithm/ Spatio-temporal sparse VGP
Algorithm Sparse spatio-temporal smoothing
Input: Data: Y), Initial params.: Y, V).
Input: Likelihood: YY, V). Space prior: ( 2,2,
Learning rates: rates:(3,p)
Time prior: (9, Q(t), ,))
while ELBO (L) not converged do
P Construct model matrices:
> CVI natural gradient step:
An = 1M, 3 (9)
q(u), f= Alg. 2( V)
Q = K(a) 2, 2 Q
E (B 18 [xxr u) [lglp(( ( [f]]]
H = IM, @ H(
Ā -1 (1- 3) 3) ++++ $ ae 8u
b Filtering and smoothing:
V (-2X(2))-
Y V
if using parallel filter / smoother then
> Hyperparameter gradient step:
q(u). L = Alg. Y V. H)
(u). / Alg. 2(Y,
else
(Btt) w) (Y [f]
qu)), L Alg. 3 V.
L -= Bu( llog N(Y u, V) + ELBO
end if
> Return posterior marginals and log likelihood: return q &
= A + P00
end while
This ELBO can be computed with linear scaling N₁: O(N,N We now show that the natural
gradient step for updating the parameters of N(Y f, V) can be computed with the same complexity.
3.2 Efficient Natural Gradient Updates
As discussed in Sec. 2.2, a natural gradiens update to the posterior, g(f) xx p(f) f, V), has superior convergence properties to gradient descent, and is completely characterised by update to the approximate likelihood, N(Y f, V), whose mean and covariance are the free parameters of the model, and implicitly define the same variational parameters as VGP. Since the likelihood factorises across the data points, these updates only require computation of the marginal distribution gffn,)) x) to obtain E. (t..,. [logp((X ,., k2,
and gradients.
As we have shown, computation of the marginal posterior amounts to smoothing over the state, f ~ N(f m, P). with the model in Eq. (11). The time marginals are given by applying the measurement model to the state: q(f.)) ((,,,)) = f [fm) df can then be obtained by integrating out the other spatial points. Given the marginal, Eq. (7) can be used to give the new likelihood parameters Y and V. The full learning algorithm iterates this process alternately with a hyperparameter update via gradient descent applied to the ELBO, Eq. (12), and has computational complexity NN,NNdd) We call this method the spalio-temporal variational GP (ST-VGP).
N(f, m, Hm,, P,, HP H), after which
3.3 Spatial Sparsity: from 2 d to O(N.M44)
We now introduce spatial inducing points, Ze in order to reduce the effective dimensionality of the state-space model. Whilst we maintain the same notation for consistency, it should be noted that the sparse model no longer requires the data to be spatio-temporal grid, only that the inducing points are. In this case, letting q(u)= N(u m (u), p(w) be the sparse variational posterior, the marginal q(f.,) = N(f mn, P..) only depends on m.
(u) p(()) due to the conditional independence property for separable kernels discussed in Sec. 1. We compute the posterior q(u) via filtering and smoothing over the state [(t) in a similar way to ST-VGP by setting,
A I M, @ A
Qn = k Z.
Q = M, 8 (())
(13)
Alg. gives the smoothing algorithm. However, the natural gradient update, Eq. (7), now becomes,
(1-1)) θ 3E (gg(u
[Ep(t] u) [log p(Y [f]]]
(14)
Op(n)
which results in ((2) and hence also V.,. being a dense matrix (i.e., V is block-diagonal) due the conditional mapping, pff, un). Therefore the approximate likelihood for the sparse model factorises
across time, but not space (see App. J for details): log N(Y u, V) E log N(Y,, un V.).

---

10
ST-SVGP. Hold
Log run time (seconds)
40000
SVGP- 2 Held
ST-SVGP- 2 Trained
SVGP-
10
ELBO
20000
10
ST-SVGP
MF-ST-SVGP- SVGP
200
400
600
800
1000
0
100
200
300
400
500
Number of time steps
leraatin number
(a) Wall-clock time comparison
(b) Negative ELBO across training itcrations
Figure 2: (a) Log wall-clock time, including any startup costs, across synthetic spatio-temporal datasets with an increasing number of time steps (average across 5 runs). (b) Negative ELBO during training for the small-scale NYC-CRIME dataset.
The Spatio-Temporte Sparse VGP ELBO Adding inducing points in space is equivalent to placing the inducing points on spatio-temporal grid (i.e., inducing points exist at time steps), and hence the variational objective directly follows from LSVGP using similar argument Sec. 3.1:
p(Y () yfF at plun / N(Y V) du
SVT. [E log
N(Y u, p(FFaTTpur
N. By(u [Epc. un) [logp( (orgp(Ynn n., [max]]] Σ Bo(u.) [log N(Yn. uns Vn)]
N log E (((O..)) Va 1) N( Hf,(t), V.)),
(15)
where the final term is given by the forward filter.
Efficient Natural Gradient Updates The marginal required to compute the ELBO and natural gradient, y(f,,x)) = [f p(f u) q()) du dfzn, = [f p(f u.,) q(u,) du,, df is the predictive distribution at input X,, n, from the posterior q(u). Because the inducing points have only been placed in space, this can be simplified through the Kronecker structure given by the state-space model. As shown in App. I, the marginal mean and covariance are,
M = =K(+) K
S. Z.Z. m
(16)
PA. K S. 2,2
K (x) K, 7, K Z,8. + K XX,XX.,, x.,,
(K! S. S. K 2, 2.S.
where m Hm,. p(") HP, are given by filtering and smoothing. By combining the above properties we see that all the terms required for the natural gradient updates and hyperparameter learning be computed efficiently in O (N,M) We call this approach the spatio-temporal sparse variational GP (ST-SVGP). The full algorithm is given in Alg.
Further Improving the Temporal and Spatial Scaling
We now propose two approaches further improve the computational properties of ST- VGP and ST-SVGP. First, we show how parallel filtering and smoothing can be used for non-conjugate GP inference, which results in theoretical span complexity of Ollog N.d3). We then present spatial mean-field approximation, which can be used independently, or in combination with sparsity.
4.1 Parallel Bayesian Filtering and Smoothing
The associative scan algorithm [9] uses divide-and-conquee approach combined with parallelisation convert sequential associative operations into log N' sequential steps (for an operator *, associa- tivity implics (a * b) (b c)). This algorithm has been made applicable conjugate Markov
7

---

PMI0 Hackaey Old Street (HK6)
30
2019-11-01
201004-01
SAAR ST.SYGP
00-
2019-00-05
3019-00-15
2019-47-20
3019-02-25
Figure 3: Observations of PM10 at site HK6, showing rich short-scale structure (top). Mean and 95% confidence of ST-SVGP trained with 30 spatial inducing points (totalling 64,770 inducing points) and SVGP with 2000 inducing points (minibatch size 100). Both models have similar training times. ST-SVGP captures the complex structure of the time series whereas SVGP smooths the data.
GPS by deriving new form of the Kalman filtering and smoothing operations that are associative [41, 13]. We give the form of these associative operators in App. and show, for the first time, how these methods can be adapted the non-conjugate setting. This follows directly from the use of the CVI approach to natural gradient VI. which requires only conjugate computations, i.e,, the linear filter and smoother. Consider the ST-SVGP ELBO:
LST-SVGP Bq(()) [Epce]
log P
Egtt) No V) + (p N(Y u, V)
factorises across time and space. compute in purallel
factorises across time. compute in paralle)
compute with parallel filter
The first two terms can be computed in parallel, since they decompose across time given the marginals g(u,). The final term can be computed via the parallel filter, and the required marginals (u,,) via the parallel smoother, which makes ST-SVGP highly parallelisable algorithm. Alg. 4 gives the filtering and smoothing algorithm and Alg. 2 shows how this method can be used in place of the sequential filter when performing inference. One drawback of the parallel filter is that when the state dimension is large, many of the available computational resources may be consumed by the arithmetic operations involved in a single filtering step, and the logarithmic scaling may be lost. Fortunately, the spatial mean-field approximation presented in Sec. 4.2 helps to alleviate this issue. In App. E provide more details on the method as well as a detailed examination its practical properties.
4.2 Spatial Mean-Field Approximation
We reconsider the state space model for the spatio-temporal GP derived in Sec. 2.1, which has process noise 9, N(O, K
(.)). This Kronecker structure implies N, independent processes
linearly mixed using spatial covariance K The linearity of this operation makes it possible to reformulate the model to include the mixing as part of the measurement, rather than the prior:
f(tn+11) A, f(t,) qa.
Yn,
pY,,
[C X ((),)((t))
(17)
where q. NO,,IN, @ Q and C(a) is the Cholesky factor of K ys (see App. F for the derivation).
This has the benefit that now both A, and Q IN, @ Q
(t) are block diagonal, such that under the prior the latent processes are fully independent. This enables mean-field assumption between N₂ latent posterior processes: (f((t) ~ II (fa (4)), where (t) is the d, di-dimensional state corresponding to spatial point Sy. This approximation enforces block-diagonal structure in the state covariance, such that matrix operations acting on the full state can be avoided. Dependence between the latent processes is still captured via the measurement model (likelihood), and experiments show that this approach still approximates the true posterior well (see Sec. 5 and App. K.6) whilst providing significant computational gains when N, is large.
5 Experiments
We examine the scalability and performance of ST-VGP and variants. Throughout, we use Matémn-/2 kernel and optimise the hyperparameters by maximising the ELBO using Adam [32].

---

Table I: NYC-CRIME (small) results.
Table 2: Test performance on matching average run time in seconds for the NYC-CRIME (large) count dataset.
ST-SVGP = SVGP when Z is fixed.
TRAIN 2 MODEL
RMSE
NLPD
MODEL
TIME (CPU)
TIME (GPU)
RMSE
NLPD
ST-SVGP
3.02
= 0.13
1.72
= 0.04
ST-SVGP
20.86
0.46
0.61
0.00
2.77
+ 0.06
1.66
+ 0.02
SVGP
1.02
1 0.13
1.72
4 0.04
MP-ST-SVGP
20.69
= 0.86
0.32
= 0.00
2.75
= 0.04
1,62
4 0.02
ST-SVGP
2.79
= 0.15
1.64
0.04
SVGP 1500
12.67
0.11
0.13
0.00
3.20
1 0.14
1.82
1 0.05
SVGP 2.94
= 0.12
1.65
= 0.05
SVGP-3000
80.80
3.42
0.45
2 0.01
3.02
= 0.18
1.76
+ 0.05
We use learning rates of p= 0.01, 8 = in the conjugate case, and p = 0.01, 3 = 0.1 in the non-conjugate case. We use 5-lold cross-validation (i.e., 80 20 train-test split), train for 500 iterations (except for AIR-QUALITY where we train for 300) and report RMSE, negative log predictive density (NLPD, see App. K.1) and average per-iteration training times on CPU and GPU. When using GPU, the parallel filter and smoother are used.
Synthetic Experiment We construct 7 toy datasets with rich temporal structure but smooth spatial structure (see App. K.2) and varying size: Nº 10, 82, 155, 300, 500, 1000, and construct 500x 100 grid that serves as a test set for all cases. As the dataset size increases we expect the predictive performance of all methods improve the expense of run time. We compare against SKI and SVGP (see Sec. 1.1). Fig. 2a shows that whilst S SVGP becomes infeasible for more than 300 time steps, the ST-SVGP variants scale linearly with time and are faster than SKI (except for the very small datasets, in which the model compile time in JAX dominates). In App. K we show the test performance of all models.
NYC-CRIME - Count Dataset We model crime numbers across New York City, USA (NYC), using daily complaint data from [1]. Crime data has seasonal trends and is spatially dependent. Accurate modelling can lead to more efficient allocation of police resources [20, 3]. We first consider a small subset of the data to highlight when our methods exactly recover SVGP. We bin the data from 1 to 10th of January 2014 (N, 10) into spatial grid of size 30 30 and drop cells that do not intersect with land (N, = 447, N = 4470). We run ST-SVGP and SVGP with inducing points initialised to the same locations. We plot the training ELBO in Fig. 2b and performance in Table 1. For fixed inducing points, both models have the same training curve and provide the same predictions (up to numerical differences). Optimising the inducing points improves both methods. A comparable inference method for non-conjugate likelihoods has not yet been developed for SKI.
We next consider observations from January to July 2014, with daily binning (N, = 182) and the same spatial grid (.gg 447, N 81,3544. We run ST-SVGP and its mean-field variant (MF-ST- SVGP) with 30 spatial inducing points (equivalent to SVGP with M = 30x = 5460). Table 2 shows that our methods outperform SVGP (with M 1500, M 3000 and batch sizes 1500, 3000 respectively) because they can afford more inducing points for the same computational budget.
Regression: AIR- AUR-QUALITY Finally, model PM10 (ug/m²) air quality across London, UK. The measurements exhibit periodic fluctuations and highly irregular behaviour due to events like weather and traffic jams. Using hourly data from the London air quality network [29] between January 2019 and April 2019 (N, = 2159), we drop sensors that are not within the London boundaries or have more than 40% of missing data (N₂ 72, N 155,448). We ST-SVGP and MF-ST-SVGP with 30 inducing points in space (equivalent to SVGP with M = 30x2159 = 64,770 inducing points). To ensure the run times are comparable on both CPU and GPU, we run SVGP with 2000, 2500, 5000, and 8000 inducing points with mini-batch sizes of 600, 800, 2000, and 3000 respectively. We run SKI with N₂ temporal inducing points and 6 inducing points in each spatial dimension.
Table 3: AIR- AIR-QUALITY regression. ST-SVGP fits the fast-varying structure well, whereas SVGP smooths the data. Average run time and standard deviation in seconds shown for a single training step. ST-SVGP and MF-ST-SVGP use 30 spatial inducing points, equivalent to SVGP with 30x2159 64,770 770 inducing points. Number of inducing points chosen make run time comparable.
MODEL (BATCH SIZE)
TIME (CPL)
TIME (GPL)
RMSE
NLPD
ST-SVGP
16.79 : 0.63
4.47
0.01
9.96 :: 0.56
8.29 :: 0.80
full spatio-temgont model
MF-ST-SVGP
13.74
1 0.49
0.85
0.01
1042 + 0.63
8.52 1 0.91
<- with mean-fich field assumptioo
SVG-20000 (600)
20.21
- 0.28
0.17 0.00
15.46 0.44
12.93 4 0.95
SVGP-2500 (800)
40.90 1 1.11
0.25 0.00
15.53 1. 1.09
13.48 1 1.85
S SVGP-5000 (2000)
119-20.00
14.20 0.44
12.73 2 0.73
baselines
SVGP-8000 (3000)
4.09 0.00
13.83 0.47
12.40
4 0.75
SKI
23.36 1.01
3.01
12.01 0.55
10.32
0.79

---

Our methods significantly outperform the SVGP baselines because they can afford considerably more inducing points. As shown in Fig. 3 the SVGP drastically smooths the data, whereas ST-SVGP fits the short-term structure well. The mean-field approach is significantly more efficient, especially when using the parallel algorithm, but we do observe a slight reduction in prediction quality.
Conclusions and Discussion
We have shown that variational inference and spatio-temporal filtering can be combined in a principled manner, introducing approach to GP inference for spatio-temporal data that maintains the benefits of variational GPS, whilst exhibiting favourable scaling properties. Our experiments confirm that ST-SVGP outperforms baseline methods because the effective number of inducing points grows with the temporal horizon, without introducing significant computational burden. Crucially, this leads to improved predictive performance, because fast varying temporal information is captured by the model. We demonstrated how to apply parallel filtering and smoothing in the non-conjugate GP case, but empirical analysis identified maximum state dimension of around d zs 50, after which the sub-linear temporal scaling is lost. However, our proposed spatial mean-field approach alleviates this issue somewhat, making the combined algorithm extremely efficient even when both the number of time steps and spatial points are large. The resemblance of our framework to the VGP approach suggests many potential extensions, such as multi-task models [4] and deep GPs [16]. We provide JAX code for all methods at https: / gi thub. com/ Aww./alll//////////a./..aaaa/.ll spatio- temporal-GPs.
6.1 Limitations and Societal Impact
We believe our work takes an important step towards allowing sophisticated GP models to be run on both resource constrained CPU machines and powerful GPUs, greatly expanding the usability of such models whilst also reducing unnecessary consumption of resources. However, when using predictions from such methods, the limitations of the model assumptions and potential inaccuracies when using approximate inference should be kept mind, especially in cases such as crime rate monitoring where actions based on biased or incorrect predictions can have harmful societal consequences.
Acknowledgments and Disclosure of Funding
This work was supported by the Academy of Finland Flagship programme: Finnish Center for Artificial Intelligence (FCAI). Parts of this project were done while OH was visiting Finland as part of the FCAI-Turing initiative, supported by HITT/FCAL OH acknowledges funding from The Alan Turing Institute PhD fellowship programme. AS and WJW acknowledge funding from the Academy of Finland (grant numbers 324345 and 339730). NAL contributed under the NVIDIA AI Technology Center (NVAITC) Finland program. TD acknowledges support from EPSRC (EP/T004 134/1), UKRI Turing AI Fellowship (EP/V02678X/I), and the Lloyd's Register Foundation programme Data Centric Engineering. We acknowledge the computational resources provided by the Aalto Science-IT project and CSC-IT IT Center for Science, Finland.
References
[1] 2014-2015 crimes reported in all 5 boroughs of New York City. adanschroederr/rrinss new- yowwyork--oo city.
http:://ww kaggle com/
[2] V. Adam, S. Elettheriadis, A. Artemev, N. Durrande, and J. Hensman. Doubly sparse variational Gaussian processes. In Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AISTATS), volume 108 of Proceedings of Machine Learning Research, pages 2874-2884. PMLR, 2020.
[3] V. Aglietti, T. Damoulas, and E. V. Bonilla. Efficient inference in multi-task Cox process models. In Proceedings of the 22nd International Conference Artificiad Intelligence and Statistics (AISTATS), volume 89 of Proceedings of Machine Learning Rexearch, pages 537-546. PMLR, 2019.
[4] M. A. Álvarez and N. D. Lawrence. Computationally efficient convolved multiple output Gaussian processes. Journal of Machine Learning Research, 12:1459-1500 July 2011.
[5] S.-i. Amari. Natural gradient works efficiently in learning. Neural Computation. 10(2):251-25 1998.
10

---

[6] S. Banerjee. Hierarchical Modeling and Analysis for Spatial Data. Chapman & Hall/CRC, Boca Raton, FL, 2004.
[7] M. Bauer, M. der Wilk, and C. E. Rasmussen. Understanding probabilistic sparse Gaussian process approximations. Advances in Neural Information Processing Systems 29 (NIPS), pages 1533-1541. Curran Associates, Inc., 2016.
[8] D. M. Blci, A. Kucukelbir, and J. D. McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 11((518):59-98) 2017.
[9] G. E. Blelloch. Scans primitive parallel operations. IEEE Tronsactions Computers, 38(1 1):1526 15388, 1989.
[10] Bruinsma, Perim, W. Tebbutt, Hosking. A. Solin, and Turner. Scalable exact inference in multi- output Gaussian processes. In Proceedings of the 37th International Conference on Machine Learning (ICML), volume 119 of Proceedings of Machine Learning Research, pages 1190-1201, PMLR, 2020.
[11] P. E. Chang, W. J. Wilkinson, M. E. Khan, and A. Solin. Fast variational learning in state-space Gaussian process models. In 30th IEEE International Workshop Machine Learning for Signal Processing (MLSP), pages 1 6, Espoo, Finland, 2020. IEEE,
[12] R. Condit. Tropical Forest Census Plots. Springer-Verlag Verlag and R. Landes Company, Berlin, Germany, and Georgetown, Texas, 1998.
[13] A. Corenflos, Z. Zhao, and S. Särkkä. Gaussian process regression in logarithmic time. arXiv preprint arXiv:2102.0990 2021.
[14] N. Cressie. Statistics for Spatio-Tempont Data. Wiley, Hoboken, NJJ, 2011.
[15] G. Da Prato and J. Zabezyk. Stochastic Equations in Infinite Dimensions, volume 45 of Encyclopedia of Mathematics and its Applications. Cambridge University Press, Cambridge, 1992.
[16] A. Damianou and Lawrence. Deep Gaussian processes. Proceedings of the Sixteenth International Conference on Artificial Intelligence and Statistics (AISTATS), volume 31 of Proceedings of Machine Learning Research, pages 207- 2052215, PMLR, 2013.
[17] A. Datta, S. Banerjee, A. O. Finley, and A. E. Gelfand. Hierarchical nearest-neighhor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, 1 1116514)800-812 2016.
[18] N. Durrande, V. Adam, L. Bordeaux, S. Eleftheriadis, and J. Hensman. Banded matrix operators for Gaussian Markov models in the automatic differentiation era. In Proceedings of the 22nd International Conference Artificial Intelligence and Statistics (AISTATS), volume 89 of Proceedings of Machine Learning Research, pages 2780-27 2789. PMLR, 2019.
[19] S. Flaxman, A. Wilson, D. Neill, H. Nickisch, and A. Smola. Fast kronecker inference in Gaussian processes with non-Gaussian likelihoods. In Proceedings of the 32nd International Conference on Machine Learning (ICML), volume 37 of Proceedings of Machine Learning Research, pages 607-616. PMLR, 2015.
[20] S. Flaxman, M. Chirico, P. Pereira, and C. Loeffler. Scalable high-resolution forecasting of sparse spatiotemporal events with kernel methods: winning solution to the NIJ "Real-time crime forecasting challenge". arXiv preprnt arXiv:1801.0255 2019.
[21] Gardner, G. Pleiss, K. Q. Weinberger, Bindel, and A. G. Wilson. GPyTorch: Blackbox matrix matrix Gaussian process inference with GPU acceleration. In Advances Neural Information Processing Systems 31 (NeurIPS), pages 7576-7586. Curran Associates, Inc., 2018
[22] A. E. Gelfand and F. M. Schliep. Spatial statistics and Gaussian processes: A beautiful marriage. Spatial Statistics, 18:86 104, 2016. Spatial Statistics Avignon: Emerging Patterns.
[23] O. Hamelijnck, W. J. Wilkinson, N. A. Loppi, A. Solin, and T. Damoulas. Data for NeurIPS 2021 submission: Tomporal Variational Gaussian Processes, 2021. URL https://doi..oooo wttp///...rr/1...2/....0.. zenodo.4531.3 4531304.
[24] J. Hartikainen, J. Riihimäki, and S. Särkkä. Sparse tempral Gaussian processes with general likelihoods. In Artificial Neural Networks and Machine Learning ICANN 2011, pages 193-200 Berlin, Heidelberg, 2011. Springer.
11

---

[25] J. Hensman. M. Rattray, and N. D. Lawrence. Fast variational inference in the conjugate exponential family. In Advances in Neural Information Processing Systems 25 (NIPS), pages 2888-2896. Curran Associates Inc., 2012.
[26] Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence (UAI), pages 282-290.. AUAI Press, 2013.
[27] J. 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 Research, pages 351 -60. PMLR, 2015.
[28] J. Hensman, N. Durrande, and A. Solin. Variational Fourier features for Gaussian processes. Journal of Machine Learning Research, 18(1), 2017.
[29] Imperial College London. Londonair London air quality network (LAQN). www.o/awww.....iii....... org. uk, 2020.
[30] E. J. Keogh and M. Pazzani. An indexing scheme for fast similarity search in large time series databases. In Proceedings of the 11th International Conference on Scientific and Statisticol Database Management, SSDBM '99, page 56, USA, 1999. IEEE Computer Society.
[31] M. E. Khan and W. Lin. Conjugate-commuuatiin variational inference: Converting variational inference in non-conjugate models to inferences in conjugate models. In Proceedings of the 20th International Conference Artificial Intelligence and Statistics (AISTATS), volume 54 of Proceedings of Machine Learning Research, pages 878 887. PMLR, 2017.
[32] D. Kingma and Ba. Adam: method for stochastic optimization. arXv preprint 111.991412.6980, 2014.
[33] K. Krauth, F. V. Bonilla, K. Cutajar, and M. Filippone. AutoGP: Exploring the capabilities and limitations of Gaussian process models. In Proceedings of the Thirty-Third Conference on Uncertainty in Artificial Intelligence (UAI). AUAI Press, 2017.
[34] A. G. d. G. Matthews, M. van der Wilk, T. Nickson, K. Fujii, A. Boukouvalas, P. Lcón- Villagil,, Z. Ghahra- mani, and J. Hensman. GPflow: A Gaussian process library using TensorFlow. Journal of Machine Learning Research, 18(40):1 6, 2017.
[35] H. Nickisch, A. Solin, and Grigorievskiy. State space Gaussian processes with non-Gaussian likeli- hood. In Proceedings of the 35th International Conference on Machine Learning (ICML), volume 80 of Proceedings of Machine Learning Research, pages 3789-3798. PMLR, 2018.
[36] A. O' O'Hagan. A Markov property for covariance structures. Statistics Research Report, 98:13, 1998.
[37] J. Quiñonero Candela and C. E. Rasmussen. A unifying view of sparse approximate Gaussian process regression, Journal of Machine Learning Research, 6:1999-1959, Dec. 2005.
[38] C. Rasmussen and C. Williams. Gaussian Processes for Machine Learning, MIT Press, Cambridge, MA, USA, 2006.
[39] Y. Saatçi. Scaluble Inference for Structured Gaussian Process Models. PhD thesis, University of Cambridge. 2011.
[40] H. Salimbeni, S. Eleftheriadis, and J. Hensman. Natural gradients in practice: Non-conjugate variational inference in Gaussian process models. In Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS), volume 84 Proceedings of Machine Learning Research, pages PMLR, 2018.
[41] S. Särkkä and A. F. on Automatic Control, 66(1):299-38) 2021.
Temporal parallelization of bayesian smoothers. JEEE Transactions
[42] S. Särkkä and A. Solin. Applied Stochastic Differential Equations. Cambridge University Press, 2019.
[43] S. Särkkä, A. Solin, and J. Hartikainen. Spatiotemporal learning via infinit--iieensional Bayesian filtering and smoothing. IEEE Signal Processing Magazine, 30(4):51-61, 2013.
[44] A. Solin, Stochastic Differential Equation Methods for -Teapoall Gaussian Process Regression. Doctoral dissertation, Aalto University. Helsinki, Finland, 2016.
12

---

[45] W. Tebbutt, A. Solin, and R. E. Turner. Combining pseudo-point and state space approximations for sum-separable Gaussian processes. In Proceedings the 37th Conference Uncertainty in Artificial Intelligence (UAI), Proceedings Machine Learning Research. PMLR, 2021.
[46] M. Titsias. Variational learning of inducing variables in 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-5744 PMLR, 2009.
[47] J. Wainwright and Jordan, Graphical Models, Exponential Families, and Variational Inference. Now Publishers Inc., Hanover, MA, USA, 2008.
[48] Wang, G. Pleiss, J. Gardner, S. Tyree, K. Q. Weinberger, and A. Wilson. Exact Gaussian processes on a million data points. In Advances in Neural Information Processing Systems 32 (NeurIPS), pages 14648- 14449,1499, Curran Associates, Inc., 2019.
[49] W. Wilkinson, M. Andersen, D. Reiss, D. Stowell, and A. Solin. End-to-end probabilistic inference for nonstationary audio analysis. In Proceedings of the 36th International Conference on Machine Learning (ICML., volume 97 of Proceedings of Machine Learning Research, pages 6776-6785 PMLR, 2019.
[50] 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), Proceedings of Machine Learning Research, pages 1747-1755, PMLR, 2021.
[51] W. J. Wilkinson, P. E. Chang. M. R. Andersen, and A. Solin. State space expectation propagation: Efficient inference schemes temporal Gaussian processes. In Proceedings of the 37th International Conference on Machine Learning (ICML), volume 119 of Proceedings of Machine Learning Research. PMLR, 2020.
[52] A. Wilson and Nickisch. Kernel interpolation scalable structured Gaussian processes (KISS GP), In Proceedings of the 32nd International Conference on Machine Learning (ICMI.), volume 37 Proceedings of Machine Learning Research, pages 1775-1784, PMLR, 2015.
[53] A. G. Wilson, E. Gilboa, A. Nehorai, and P. Cunningham. Fast kernel learning for multidimensional pattern extrapolation. In Advances in Neural Information Processing Systems 27 (NIPS), pages 3626-3634. Curran Associates, Inc., 2014.
13


Reviews:
[{"review_id": "mw_Nm9Y7ih5", "reviewer": "Reviewer_HABe", "summary": "The authors propose a novel combination of existing strategies for scalable inference for spatio-temporal GPs. The techniques utilized were Kernel separability for convenient Kronecker structure, inference by filter-smoother in the  time domain, inducing point methods, and what appear to be 2 original contributions, a mean-field approximation to their spatial latent processes, and some computational strategies leveraging parallelization of the optimization. The authors make some camparisons to other methods in terms of computation time and accuracy.", "questions": "", "limitations": "", "rating": 7, "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": "Overall, the paper is well written and organized, although there are a number of minor details which I will list below that should to be addressed where the authors give insufficient details, introduce notation that is not explained (other than in the supplement), or possibly have a typo. Otherwise, the authors do a good job of a difficult task, synthesizing a wide range of GP literature to set the proper context for the their approach and to understand the basic advantages that the techniques confer.\n\nMy main complaint in terms of content is with respect to Section 5 (Experiments). The authors conduct a number of experiments, 1 with synthetic data and 2 on real data. An obvious contender with their methods is Wilson & Nickisch's SKIP method, which also makes use of a regular grid for the inducing points. The authors show that SKIP compares favorably in terms of wall-clock time but do not compare it's performance otherwise except for an example in the supplement, which is insufficient for this analysis. It would be better if the authors could include SKIP performance in the main-paper results. \n\nMinor points\n- Have the authors made a comparison with Hensmen et al.' variational Fourier feature method? This also seems to be an obvious contender with the author's approach.\n- Line 66- What is D_s?\n- Line 88 - \"Where d_t is the dimension of the state-space model induced by \\kappa_t...\" Could the authors elaborate on this? Why should the dimension at each spatial point be larger than 1?\n- Equation (4) The authors introduce notion in the superscripts of matrices that is not explained. Although they describe the notation in the supplement an explanation of the superscripts needs to appear in the main text as well.\n- Lines 134-135 - \"This shows that the natural gradient is completely characterized bu updates to an approximate likelihood\" I'm sorry but you kind of lost be here. It's not obvious at all from this explanation why that the latter should follow from the former. Please explain.\n- Line 136 - What do the superscripts on the \\lambda 's mean?\n- Equation (11) - What do the superscripts (t) on the A, K, and Q matrices mean? Please explain in the main text.\n- line 249 - The authors are making a claim about a mean field approximation, which I would have expected to translate into factorization (i.e. product) of _densities_ over spatial points but the have written a product of functions over spatial points. Did I just miss something huge or is this a typo?\n ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have not addressed any limitation and potential negative societal impact.  They might consider how scalability reduces resource consumption.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "XPNKZMGDtu", "reviewer": "Reviewer_MX1g", "summary": "The paper presents a novel approach to variational inference in Gaussian Processes for spatio-temporal data (including non-gaussian likelihood), using state-space representation for linear scaling in the temporal dimension. The authors show how using natural gradients and conjugate-computation VI approach leads to decomposition of the ELBO which, together with linear-time filtering and smoothing of the state-space representation, results in significant computational gains.\n\n\nThe results are supported experimentally on one synthetic and 2 real datasets (conjugate and non-conjugate data). The proposed models outperform the baselines (sparse variational GP).\n\nOverall this is a very solid paper.", "questions": "", "limitations": "", "rating": 7, "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": "\nThe problem of computational complexity in GP inference is still an important one and, as this work shows, there is still more room for improvement. This paper, while not groundbreaking, definitely provides substantial value to the community.\n\n### Strengths: \n\nThe results of the paper are useful. The typical inducing point approach does not work well in the temporal dimension and the number of inducing points needed increases proportionally, still resulting in effectively cubic complexity. With this in mind, the linear complexity in the temporal dimension proposed in this work provides significant computational benefits.\nThe authors propose further improvements to the method via spatial mean-field approximation and a parallel version of the filtering algorithm.\n\nThe theoretical results are clear. I have not checked all the math thoroughly, but the results look correct. Detailed proof and derivations are given in the appendix.\n\nThe experimental results are thorough. Both RMSE and predictive likelihood are shown with error bars, the computation time is given and the experimental settings are described in detail.\n\nThe writing is very clear and easy to follow.\n\n\n### Weaknesses:\n\nRelated work is focused on spatio-temporal modelling and missing some other relevant references. Notably, Vincent et al. [1] proposed doubly stochastic variational GP, merging the best of two main approaches to GP inference: variational methods and state-space representation, by introducing inducing states in the state-space representation. Since the submitted work addresses the same problem of making inference in GP more efficient, the connection to [1] should be explained.\nAnother possibly relevant line of work is the Variational Fourier Features approach [2]. Unlike typical inducing points, VFF would be suitable for capturing temporal dependencies. Some discussion on this is desirable.  \n\nExperimental comparison is only provided to the sparse GP. If the (missing) relevant work is directly applicable to spatio-temporal data, a comparison is needed.\n \n\nThe work is quite dense in theory, which resulted in placing important details in the supplement. Algorithms 3 and 4 are referred to in the text and in alg 2, but can only be found in the appendix. At least explicitly refer to the appendix when talking about alg 3, 4.\n\n\n\n[1] Adam, Vincent, et al. \"Doubly sparse variational Gaussian processes.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2020.\nAPA\t\n\n[2] Hensman, James, Nicolas Durrande, and Arno Solin. \"Variational Fourier Features for Gaussian Processes.\" J. Mach. Learn. Res. 18.1 (2017): 5537-5588.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The limitations of the work are discussed sufficiently. Societal impact is not addressed due to the general theoretical nature of the work. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "V6OAjZaMlry", "reviewer": "Reviewer_64f6", "summary": "The paper presents an approach for GP-Regression in the case of spatio-temporal data. The paper discusses a new scalable approach in the temporal domain by utilizing Bayesian filtering and smoothing in conjugation with conjugate-computation variational inference. Further scalability in the spatial domain is established by the use of an inducing point method with and without a mean field assumption.", "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": "All in all, I enjoyed reading the paper, and the approach seems like a nice addition to the computational library for GP-Regression.\n\nOriginality: The use of filtering and smoothing for GP-Regression is a well-known technique. Spatio-temporal GPs have been discussed in this context extensively as well. However, as the authors pointed out the scalable approach in the spatial domain is missing.\n\nQuality: I think in general the paper is technically sound. However, when first reading the paper I was not so sure about it. The reason is that there are a lot of mathematical blunders and ambiguities, which should definitely be fixed. The claims are evaluated empirically, which is standard in this field of research. The method itself makes definitely sense in this context.\n\nClarity: I think this is one of the biggest weaknesses of the paper. Organization could definitely be improved and I oftentimes had a bit of a hard time following the discussed steps. But in general, I think the included background is informative and well selected. Though, I could see people having trouble understanding the state-space GP-regression when coming from the more machine-learning like perspective of function/ weight space GP-regression.\n\nSignificance: I think there is definitely a significance in this work, as GP-Regression is usually a bit problematic because of the scaling, though it is still used extensively in certain areas, such as Bayesian Optimization or for modelling dynamical systems in robotics. \n\n\n•\tBackground: $f$ is a random function which describes a map e.g. $f:  \\mathbb R \\times \\mathbb R^{D_s} \\rightarrow \\mathbb R$ and not a function of $\\mathbf X \\in \\mathbb R^{N_t \\times N_s \\times D_s}$ as described in eq 1. At least, when one infers the inputs from the definitions of the kernel functions.\n\n•\tIn general, the definitions are confusing and should be checked,e.g. check if $f_{n,k}=f(X_{n,k})$ is correct and properly define $X_{n,k}$.\n\n•\tThe operator $\\mathcal L_s$ is not mentioned in the text\n\n•\t2.1: The description of the process $\\bar{f}$ is confusing as the relationship to the original process $f$ is established just at the end. \n\n•\tIt would be helpful to add a bit more background on how the state space model is constructed from the kernel $\\kappa_t(\\cdot,\\cdot)$, e.g. why it induces the dimensionality $d_t$ and also describe the limitations  that a finite dimensional SDE can only be established, when a suitable spectral decomposition of the kernel exist.\n\n•\tIt should be mentioned that $p(y \\mid H \\bar{f}(t_n))$ has to be chosen Gaussian, as otherwise Kalman Filtering and Smoothing and CVI is not possible. Later on in the ELBOs this is assumed anyway.\n\n•\t$p(u)=N(u \\mid 0, K_{zz})$ is a finite dimensional Gaussian not a Gaussian process and p(u) is not a random variable ($=$ not $\\sim$).\n\n•\tThe notations for the covariances, e.g. $K_{zz}$ are discussed in the appendix. I am fine with it; however, it should be referenced as I was confused in the beginning.\n\n•\t2.2: The $\\log$ is missing for the Fisher information matrix.\n\n•\tThe acronym CVI is used in the paragraph headline before the definition.\n\n•\tSome figures and tables are not referenced in the text, such as figure 1. \n\n•\t3.1: In line 173 the integration should be carried over $\\bar{f}$ and not $s$, I guess?\n\n•\tI had a bit of a hard time establishing the connection $\\mathbb E_{p(f)}[N(\\tilde Y \\mid f, \\tilde{V})]=p(\\tilde Y)$ which is the whole point why one part of the ELBO can be calculated using the Kalman filter. Adding this to a sentence to the text would have helped me a lot. \n\n•\tOne question I had was that for computing the ELBO the matrix exponential is needed. When backpropagating the gradient for the hyper parameters, is this efficient? As I am used to using the adjoint method for computing gradients of the (linear) differential equation.\n\n•\tReference to Appendix for the RMSE and NLPD metrics is missing.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A", "review_text": "", "extra_scores": {"ethical_concerns": "N/A", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "8"}, "extra_text": {}}, {"review_id": "IHeoFVs8Zln", "reviewer": "Reviewer_LdLE", "summary": "This paper presents an efficient algorithm for variational inference of spatio-temporal processes using Gaussian processes as the variational family. The paper leverages so-called conjugate-computation variational inference to compute/represent distributions in the approximating variational family. By leveraging Markov, separability of the covariance kernel and prior sparsity assumptions on the spatio-temporal process, the paper demonstrates that the cubic computational complexity can be improved to linear in terms of the time variable. ", "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": "Overall I found the paper to be an acceptable contribution to the literature. The ideas are pretty standard, but used in a sensible way. In terms of originality, I would argue that some of the ideas, particularly around using smoothing to compute the posterior and natural gradients, are not new, though they have not been applied to the current setting. The paper is well written and reasonably clear. The computational speed up noted above should be of significance, particularly for large spatio-temporal data sets. \n\nOne thing I would ask the authors to improve upon is the exposition on the sparse priors, and in particular improve the description of the supporting/inducing points. I would also urge them to be more specific in how the SPDE is related to the GP prior (how is f related to \\bar{f})? Also, provide a citation for why the SPDE reduces to an SDE in their setting. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors do not address the limitations of their method or any negative societal impact. This could be important since their example dataset is of crime statistics in NYC. How does the bias in this method affect any conclusions drawn from the method?", "review_text": "", "extra_scores": {"ethical_concerns": "None.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors of this paper improve on the computational complexity of inferring spatiotemporal GPs. Specifically, via a combination of sparse priors, inducing points, and Markov separability they reduce a cubic cost to a linear cost. The reviewers all agree that the work is clear and stands to benefit the community. Moreover the few remarks raised by the reviewers seem to be readily addressed in revision. I therefore am happy to recommend this paper be accepted at NeurIPS.
Decision: Accept (Poster)