Venue: NeurIPS 2021 Conference

Reviews:

final_score: rating: 6: Marginally above the acceptance threshold
review_content: 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. 
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. 

One 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. 
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_title: null

final_score: rating: 7: Good paper, accept
review_content: 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.


The 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).

Overall this is a very solid paper.
main_review: 
The 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.

### Strengths: 

The 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.
The authors propose further improvements to the method via spatial mean-field approximation and a parallel version of the filtering algorithm.

The 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.

The 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.

The writing is very clear and easy to follow.


### Weaknesses:

Related 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.
Another 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.  

Experimental comparison is only provided to the sparse GP. If the (missing) relevant work is directly applicable to spatio-temporal data, a comparison is needed.
 

The 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.



[1] Adam, Vincent, et al. "Doubly sparse variational Gaussian processes." International Conference on Artificial Intelligence and Statistics. PMLR, 2020.
APA	

[2] Hensman, James, Nicolas Durrande, and Arno Solin. "Variational Fourier Features for Gaussian Processes." J. Mach. Learn. Res. 18.1 (2017): 5537-5588.
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_title: null

final_score: rating: 7: Good paper, accept
review_content: 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.
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.

My 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. 

Minor points
- 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.
- Line 66- What is D_s?
- 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?
- 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.
- 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.
- Line 136 - What do the superscripts on the \lambda 's mean?
- Equation (11) - What do the superscripts (t) on the A, K, and Q matrices mean? Please explain in the main text.
- 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?
 
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_title: null

final_score: rating: 7: Good paper, accept
review_content: 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.
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.

Originality: 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.

Quality: 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.

Clarity: 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.

Significance: 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. 


•	Background: $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.

•	In 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}$.

•	The operator $\mathcal L_s$ is not mentioned in the text

•	2.1: The description of the process $\bar{f}$ is confusing as the relationship to the original process $f$ is established just at the end. 

•	It 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.

•	It 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.

•	$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$).

•	The 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.

•	2.2: The $\log$ is missing for the Fisher information matrix.

•	The acronym CVI is used in the paragraph headline before the definition.

•	Some figures and tables are not referenced in the text, such as figure 1. 

•	3.1: In line 173 the integration should be carried over $\bar{f}$ and not $s$, I guess?

•	I 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. 

•	One 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.

•	Reference to Appendix for the RMSE and NLPD metrics is missing.

limitations_and_societal_impact: N/A
review_title: null

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)