Robust and Conjugate Spatio-Temporal Gaussian Processes
Authors: William Laplante
Matias Altamirano
Andrew B. Duncan
Jeremias Knoblauch
Francois-Xavier Briol
Venue: ICML 2025 poster (2025)
Keywords: Gaussian Processes
Robustness
Spatio-Temporal Analysis
Generalised Bayes
Primary area: probabilistic_methods->gaussian_processes

Abstract: State-space formulations allow for Gaussian process (GP) regression with linear-in-time computational cost in spatio-temporal settings, but performance typically suffers in the presence of outliers. In this paper, we adapt and specialise the *robust and conjugate GP (RCGP)* framework of Altamirano et al. (2024) to the spatio-temporal setting. In doing so, we obtain an outlier-robust spatio-temporal GP with a computational cost comparable to classical spatio-temporal GPs. We also overcome the three main drawbacks of RCGPs: their unreliable performance when the prior mean is chosen poorly, their lack of reliable uncertainty quantification, and the need to carefully select a hyperparameter by hand. We study our method extensively in finance and weather forecasting applications, demonstrating that it provides a reliable approach to spatio-temporal modelling in the presence of outliers.

Robust and Conjugate Spatio-Temporae Gaussian Processes
William Laplante 123
Matias Altamirano Andrew Duncan" 34
Jeremias Knoblauch
François-Xavin Briol
Abstract
t=0.53
t=0.67
0.80
Data
State-space formulations allow for Gaussian pro- cess (GP) regression with linear-in-tinn computa- tional cost in spatio-temporal settings, but perfor- mance typically suffers in the presence of outliers. In this paper, we adapt and specialise the robust and conjugate GP (RCGP) framework of Altami- rano et al. (2024) to the spatio-tempone setting. In doing so, we obtain an outlier-robust spatio- temporal GP with computational cost compa- rable to classical spatio-temporal GPs. We also overcome the three main drawbacks of RCGPs: their unreliable performance when the prior mean is chosen poorly, their lack of reliable uncertainty quantification, and the need to carefully select hyperparameter by hand. We study our method extensively in finance and weather forecasting applications, demonstrating that it provides a re- liable approach to spatio-temporal modelling in the presence of outliers.
1.0
Truth ST-RGGP STGP
0.0
00.5
1.5
-20
2.5
3.0
-1
Figure 1. Spatio-temporal GPs in the presence of outliers. The top row shows the observed data function of spatial covariates IR? Outliers are highlighted in red and are uniformly distributed [-8, -6] U [6, 8). The second row gives the fit of regular STGP, whilst the third row gives our proposed ST-RCGP fit. The last row shows the true latent function ( $2 1) (in()nt) cos(2rt))s) Further details are provided in Appendix C.4.
& Witmer, 2018). Their popularity arises from their ability to encode spatial and temporal properties such as smooth- ness, periodicity, and sparsity (Duvenaud, 2014), allowing them model phenomena such as local weather patterns or seasonality. Crucially, GPs also have an exact, closed- form posterior when using a Gaussian likelihood. However, naive implementations have a cubic computational cost in the number of data points N, limiting their scalability, To address this issue, a plethora of approximations have been proposed (Drineas et al., 2005; Titsias, 2009; Hensman et al., 2013; Wilson & Nickisch, 2015). While effective, these do not typically recover the exact GP, require careful tuning, and can degrade performance for complex datasets (Bauer et al., 2016; 6; Pleiss et al., 2018)
1. Introduction
Gaussian processes (GPs; Williams & Rasmussen (2006)) flexible probabilistic models used in vast class of prob- lems from regression (Williams & Rasmussen, 2006) to emu- lation (Santner et al., 2018) and optimisation (Garnett 2023). GPs originated in spatial statistics, where their use for re- gression was known as kriging (Krige, 1951; Stein, 1999), but they have more recently been used widely in spatio- temporal settings, including in epidemiology (Senanayake al., 2016), neuroimaging (Hyun al., 2016), object track- ing (Aftab et al., 2019), and psychological studies (Kupilik
'Deparmment of Physics and Astronomy, University Col- lege London, London, United Kingdom *Department of Sta- tistical Science, University College London, London, United Kingdom 'The Alan Turing Institute. London, United Kingdom "Department of Mathematics, Imperial College London, Lou- don. United Kingdom. Correspondence to: William Laplante <william. Cailliam.lapaant..A 24w.el.c..k....
In spatio-temppral settings, an alternative strategy is to re- formulate the GP as a state-space model (SSM) (Reece Roberts, 2010; Hartikainen Särkkä, 2010; Sarkka & Har- tikainen, 2012; Solin, 2016; Nickisch et al., 2018; Hameli- jnck et al., 2021). This gives rise to class of models- spatio-temporaa Gaussian processes (STGPs)-wih linear cost in the number temporal observations. But, as with standard GPs, STGPs lack robustness to model misspeci- fication, such with outliers arising from extreme events
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Robust and Conjugate Spatio-Temports Gaussian Processes
(Heaton et al., 2011), measurement errors that are spatially correlated (Tadayon & Rasekh, 2019), and other hetero- geneities (Fonseca et al., 2023). Figure I illustrates this in the case of outliers. Clearly, the STGP (2nd row) fails to align with the ground truth (4th row).
2. Background
GP Regression Let (xx, yh J be observations, where X X S R² are covariates and ya E y S RR are responses. For observation noise Ek. GP regression considers
To resolve this issue, existing work on STGPs has focused on using likelihoods corresponding to distributions which are more expressive than Gaussians, such as mixtures or heavy-tailed distributions. Doing so breaks conjugacy, and the posterior must typically be approximated (Hartikainen et al., 2011; Solin & Särkkä, 2014; Hamelijnck et al., 2021). We refer the reader to Nickisch et al. (2018) for a com- prehensive overview and to Wilkinson et al. (2023) for a Python package providing those algorithms. Although these methods have been efficiently implemented, they typically use additional optimisation steps at each time-point and are hence significantly more costly than conjugate STGPs. We also highlight small body of work that uses outlier- rejection Kalman filters with STGPs to improve robustness in inference tasks (Bock et al., 2022; Waxman & Djurić, 2024). While these outlier-reeection methods offer robust- ness at lower computational cost compared to mon-conjugate STGPs, they are generally considered less expressive and not as strongly supported by theoretical foundations.
- (x+)) + Exs for k = I...., N
(I)
where the latent function f X' + R is modelled by a GP prior f ~ GP(m, k:) with mean ance (or kemne)
X - R and covari- X x X + R. such that for any inputs X = (x1..., xx)] the vector f = (f(x))..., f [xx is distributed as N-dimensional Gaussian (mm, K) with m = (m(x1),.... m and (K), = x(x,,,,) x,). The kernel K typically depends on parameters @ E that we omit from the notation for brevity. When the observa- tions y [,1,..,,,)) have independent Gaussian noise (EI, (EI,...,EN) EN) ~ (N0 LLL),, the posterior predictive for ft = fxx)) at x+ E X is faly. X ~ NN EGp) with
par = m. + k, (K+ + o INNNNN (y (y-m) m)
(2)
ΣGp = ker - k. (K + o ININ) -lk,
where IN is the N x N identity matrix, k. (I((2,,1) .1), (())),,,,(((,,................. K(E,,, eN)) his. = x(I,,, and m. m(x.). The orange colouring can be ignored for now, but will be used to highlight differences with RCGPs. To obtain this mean and covariance, we must invert an N x N matrix- an operation with computational complexity O(N3).
Recently, Altamirano et al. (2024) introduced a method, called robust and conjugate Gaussian processes (RCGPs), that uses generalised Bayesian inference (Bissiri et al., 2016; Knoblauch et al., 2022) to confer robustness to standard GPs. Their approach is highly attractive since it provably provides robustness outliers whilst maintaining conjugacy, but it shares the cubic cost of GPs. Existing work on RCGPs also has three main limitations: RCGPs underperform when the prior mean is chosen poorly (see Ament et al. (2024) and the Appendix of Altamirano et al. (2024)), their uncer- tainty quantification properties have not been well-studied, and they have an additional hyperparameter compared to standard GPs. Further, the existing heuristic for tuning this parameter relates to the proportion of outliers in the data, which is unknown in practice and generally has to be hand- picked on case-by-case basis.
State-Space Formulation An alternative approach, which has linear- in---iintttmee cost, is to use a state-space representation. Consider Equation (1) with spatio-tempore inputs. At time ta E T S R. we now have imputs Xx., = (s,, fx) E X S T on spatial grid S = (Si,...,,,,, sn,)) S C
R",xd, with n. points and spatial dimensionality d, (i.e.
d = dg + 1). The observations corresponding to X2- (x ) are denoted yı (WVk...... m,.,,,) ) E R leading to a total number of data points N = n44 where n, is the number of time steps. If standard GPs were used here, the cost would be On?in?). which may become impractical.
In this paper, we show how to refine and specialise the RCGP framework for spatio-temporal data, Our algorithm, denoted spatio-temporal RCGP (ST-RCGP), inherits the computational and memory efficiency of STGPs, as well as the robustness properties of RCGPs (see 3rd row of Fig- ure 1). The sequential aspect of the state-space formulation also allows us to overcome the three main limitations of RCGPs (sensitivity to prior mean, lack of reliable uncer- tainty quantification, and the additional hyperparameter). Overall, we observe that the ST-RCGP provides inferences comparable state-of-the-art non-Gaussian STGPs with computational cost similar to classical STGPs.
A solution to this issue is to reformulate GPs as SSMs. We first collect partial derivatives of f with respect to time in a state vector z = z(s, t) E g(++1) where (z(s, f)); = o² / () for 1, 1,..., v+ 1. We assume stationary and separable kernel that K(s, t; t') = ks s')xx(t-t'),, where K2 and Kr are spatial and temporal kernels respec- tively. For a large class of kernels (Solin, 2016), we can represent the GP prior as the solution to a stochastic differ- ential equation (SDE):
0z(s, Dt
Fzz(s + L.w(s,t)
(3)
2

---

Robust and Conjugate Spatio-Temports Gaussian Processes
where Fe E R(+ L, E R(++1)x1 and w(s.) () is a spatio-temporal white noise process corresponding to the derivative of Brownian motion with spectral den- sity Qe E I [x((++1)) (Solin, 2016). The matrices Ft, L, Qu, and the constant v depend on K. Equation (3) defines a continuous latent process, but given a finite collection of points, this becomes a SSM with states Z):- = (z(s), tx), z(S,,, tx)) R=,(++1) and time steps A4, = f - (Hartikainen & Särkkä, 2010; Sarkka & Hartikainen, 2012; Särkkä & Solin, 2019):
Varying y
Varying /
Varying C
1.0-
0.6-
0.0
9
0
21 (3L)
c- 0.1
0.5
c 20
3 0.1
0.5
3 1.0
Figure 2. Behaviour of y ++ wIMO y) as we vary c, B, and 7 We emphasize 1yke - (xxx)) ) and Wino (Xm,(XX, Bk) with datum yku
20 (0, Eo); and for k >>0 0,
Z 1 + 9k-1. 9k-1 ~ (O, 11-
(4)
Generalised Bayes for GP Regression Conjugate GPs and STGPs rely on the fragile assumption that c1...., CN are Gaussian. While non-Gaussian likelihoods can address this, they also break conjugacy- and thus require expensive and potentially inaccurate approximations. Fortunately, gen- eralised Bayes (GB) (Bissiri et al., 2016; Knoblauch et al., 2022) can produce robust and conjugate GP posteriors. GB mitigates model misspecification by constructing posterior distributions through a loss L : J.N x y.V x X + RR via
Yk = Hz; c₁
where HE E (n, defined such that Hz;- = fa == (f(11, ..., (e),,,,)...,,........ f(sn,, +4)) The matrices A and E depend on Fr, L, and Qe,t in the SDE formulation. We de- fine them in Appendix C.2, and provide information on how to compute them in practice for a list of common kernels k.
Filtering and Smoothing Solving the SSM from Equa- tion (4) via sequential inference amounts to first retrieving the filtering distribution ((zuyy1:)) and then the smooth- ing distribution p(zz ly).,) (Särkkä, 2013). Taking Equa- tion (4), and assuming the observation noise is Gaussian that cg- ~ (0, I,,n), ), the predict and update equations are conjugate and given by the Kalman filter (Kalman, 1960) and Rauch-Turre Stricbel smoother (Rauch et al., 1965); see Section 8.2 of Särkkä (2013). The resulting distribution has densities p(zz V1:k-1) = NV(z majk-1, Majk-1,PajE-1)) and p(zaly 1x:k) N majk; PAJH)), with predict step:
yf x p(f exp NNC
(7)
where x denotes proportionality and & is typically an esti- mator of a statistical divergence between the likelihood model and the data (Jewson et al., 2018). The loss C used for RCGPs (Altamirano et al., 2024) is a modifi- cation of the weighted score-mutching divergence, first proposed by Barp et al. (2019), to regression. It gener- alises score-matching (Hyvärinen, 2005), and is a special case of the framework in Lyu (2009); Yu et al. (2019); Xu et al. (2022). Importantly, it uses weight function w : X x y + IR \ {0) that down-wett weights unreasonably extreme data points. This choice yields the RCGP posterior predictive f. (MRCOPP: ERCGP) with
m == A Je-1m6-1k-1
(5)
P PAJK-1 == P 1jt- A Ex-1, and (Bayes) update step:
MRRGPP = m. + k. (K + 0² =JJw) -1 (y (y-man mm))
P = (PH)- + H -Ho-21 - HI,
(8)
DROP = ki. - k, (K + )w)w) -2k.,
K₁ : PAI.HTT IIn,
(6)
where w = (w(x1, M1): , (XN,, )) and
milk == m + KR - fe);
mw m + o²V log(w
where i. = Hm Ka E (=) majk E (n.((+1) and PæIK E (M+1)) The posterior p(zxlyi:w is obtained by marginalising the smoothing dis- tribution, and matches exactly the original GP posterior (Solin, 2016). For a new input t.). ). the predictive is obtained by including (s., t.) in the filtering and smoothing algorithms and predicting without updating. Importantly, this algorithm only inverts matrices whose sizes scale with where typically v < 10 (Hartikainen & Särkkä, 2010) or n,, but not with n4- In contrast to the standard GP's (nii?) time and On) memory cost, the STGP only requires ( time and (nn??) memory (Solin, 2016).
0
Jw diag
2
The expressions look similar to those in Equation (2), with differences highlighted in orange. This also reveals that standard GPs are recovered for the constant weight w(x. v) a o/V2. The similar forms also suggest that, much like classical GPs, naive implementation of RCGPs has a complexity of O(V3) due to the associated matrix inversion. However, under mild conditions on w, RCGPs benefit from provable robustness to outliers (sce Proposition 3.2 of Altamirano et al., 2024).
3

---

Robust and Conjugate Peatio-eemmpral Gaussian Processes
1.0
2.5
Coverage
0.8
0.6
0.0
0.4
2.5
0.2
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.2 0.4 0.6 0.8 1.0 Quantiles
mi(2) 0 (Bad prior mean)
RCGP (m,, c= Qx(1 - c))
RCGP (ma, c= 0.8)
ma(a) sin(z) (Good prior mean)
RCGP (mg Qx(1 - -))
Data
Figure 3. Failure Modes with Existings RCGPs. We generated N = 80 data points from the true function f(x) = 3 s (2rzz using (0, 0.5) additive noise, then corrupted 10% of observations around
0.76. The lefimost figure shows that RCGP's posterior mean is strongly affected by the choice of prior mean, highlighting Issue #1. In the two middle plots, we demonstrate that the recommended c= @1 - <) = 3.7 (from Altamirano et al. (2024)) is for from optimal emphasizing Issue #3-even when fitting simple data and knowing the true proportion of outliers These two plots also show RCGP's inaccurate uncertainty quantification when hyperparameters are chosen poorly, which is further supported by the rightmost coverage plot and illustrates Issue #2
In prior work, Altamirano et al. (2024) suggested to take w of the form of an inverse multi-quadri (IMQ) kernel depending on centering function y X + R. shrinking
mance of RCGPs when m is not chosen carefully. Indeed, setting 7xx) == m(x) can be undesirable when m is not a good approximation of f. In that case, observations close to mxx) but far from f(x) will have large weights despite likely being outliers, whereas observations close f(x but far from m(x) will be down-weighted despite likely having little noise. We show this in Figure 3. Altamirano et al. (2024) suggest fixing this problem by choosing a prior mean using a simpler regression model, but this solution requires an additional fit of the data, which can be cumbersome.
function e: X + R2 and constant / > 0:
(y (x)) -
wIMO(x, v) = (1+
(9)
c(x)²
The IMQ has been used frequently to robustify score-based divergences (Barp et al., 2019; Key et al., 2021; Matsub- ara al., 2022; Altamirano et al., 2023; Matsubara al., 2024; Liu & Briol, 2024). It is a bump function centred at 7(x) that grows ly [x)] decreases, and shrinks ly - (x)] increases (see Figure 2). The centering function (x) determines where the bump is maximised, and any far from >(x) is down-weighted (see Figure 2). The shrink- ing function c(x) determines how quickly we weeght observations that deviate from 7(x), and |3 determines the maximum of the weight function (see Figure 2). Note that 3 is multiplicative constant which is equivalent to GB learning rate (Wu & Martin, 2023).
Issue #2 (Poor uncertainty quantification): The values of 33 and have a significant influence on the predictive variance, but it is not clear that the suggested choices are sensible when it comes to uncertainty quantification, and Al- tamirano et al. (2024) did not study this question. This could mean that RCGPs are consistently under- or overconfident- problem we also illustrate in Figure 3. Sinaga al. (2024) proposed computation aware extension of RCGPs to im- prove uncertainty quantification, but their approach still depends on having good values for ,3 and c.
While IMQ-based weights have shown much promise, they rely on three hyperparameters that can be difficult to select: y, 3, and c. Altamirano et al. (2024) proposed centering at a o/v2 to guarantee that observations are never assigned larger weight than they would in standard GP. Finally, they sug- gested a heuristic choice for based on the assumed propor- tion of outliers: if E E [0, 1] is the expected proportion of outliers, they suggest setting the shrinking function to the constant c(x) = (1 - <), for QN - the (1 e)th quantile of (ly+ - )
Issue #3 (Selection of shrinking function c): The pro- posed heuristic for selecting e requires the user to speculate the proportion e of outliers. Not only is this difficult do reliably if we do not know but, shown in Figure 3, the approach can also fail to select good values of even when the correct is used. This can lead to unreliable posterior estimates and either under- or overconfidence in the RCGP posterior. addition, setting x + c(x) to be constant can be sub-optimal when outliers are clustered in time or space, in which case adapting the shrinking function according to its input could lead to improved uncertainty quantification.
the prior mean via (x) m(x). and setting
These choices can provide good performance in many set- tings, but introduce several failure modes:
3. Methodology
Issue #1 (Sensitivity to the prior mean m): Altamirano et al. (2024); Ament et al. (2024) pointed to the poor perfor-
RCGPs We now show that, similarly to GPs, inference for RCGPs can be performed using an SSM

---

Robust and Conjugate Spatio-Tempore Gaussian Processes
representation and filtering/smoothing updates. To achieve this, we the weighted score-matching loss
Proposition 3.1 offers two key advantages over standard RCGPs: firstly, it operates with significantly smaller matri- ces of size n. (v + 1) x n.(v + 1) rather than nt X R as for vanilla RCGPs; speeding up computations. Secondly, whilst vanilla RCGPs require the weight function to be fixed during inference, we can adapt it as filtering is performed, thereby improving inference quality. This novelty powers much of our developments, and allows us to address Is- sues #1, #2, and #3 mentioned in Section 2. Our notation highlights this adaptivity by indexing weights with time.
1
n
YR. = The
1 ( 3k. (a (HK.s) l12
i=1
2V (x. yk. ) fc. (wk=)) (10)
where sy sy( := (f- y)/o /r is the score of the Gaussian error, and the dependency on 2y is through fk.j == (H Then, the generalised filtering posterior equations are
The standard formulation of RCGPs arises as particular case of Proposition 3.1 where Wa are specified as outlined in Section and are only allowed to depend on the cur- rent observation (xt: filtering solutions of the ST-RCGP match the RCGP batch solutions exactly. This emphasizes the importance of adap- tivity: ST-RCGPs can fix prevailing issues of vanilla RCGP precisely because they allow adaptive weights, which in turn causes filtering and smoothing distributions to differ. Propo- sition 3.2 formalizes this result, and identifies conditions for which ST-RCGP and vanilla RCGP match.
pe pe(2xly (1)) x
exp (n ((x,, Yk. zx)), pe pe(Zalyi---1) = / p(zz
dzz-1
In this setting, smoothing and
for p(z 2 1) defined as for Equation (4).
Proposition 3.1. Suppose we model E1. EN....,CN ild NN(00,0)) o?) and choose w : X x y + \ such that y +> w(y,x) is differentiable. Then, the generalised posterior Pc Pclzgly (Zg)yia)) associated to the SSM formulation of RCGPs (Zg: PGlk) klk) and inference can be performed through Kalman fitering/smoothing with predict step
is N N(z;; milk
Proposition 3.2. Consider constructing a generalised pos- terior using L as defined in Equation (10), and assume that weights are chosen non-adaptively. Then, whenever the GP prior is a Gauss-Marke ing and smoothing distributions constructed sequentially as in Proposition 3.1 will be equal to those constructed in a single batch from using the vanilla RCGP form defined in Proposition 3.1 of Altamirano et al. (2024).
GB mijk-1 == Ak-1-- k-1kk
process, it will hold that filter-
P PaR-1:== = A
GB k- I|k
AK + Σx-1,
and update step
DGB P k1/k 1 + H σ 2J H
See Appendix B.2 for more details and the proof.
KG : PGR H -22 2--1
W
GB
== m + KCB (y
Weight Function Proposition 3.1 imposes some con- straints on the choice of weight function: it should be strictly positive, differentiable over its domain (this ensures that all quantities in Proposition 3.1 are well-defined). We also need the weight to be bounded from above to ensure robustness, and to avoid attributing weight to arbitrarily large E IR we require that / wxx, y) = The IMQ satisfies all these conditions, and has been well-studied and recom- mended in prior literature (Matsubara al., 2022; Riabiz et al., 2022; Chen et al., of weight function. This choice is further justified by the fact that the IMQ hyperparameters 7, f3 and e are related to concepts of robustness that we will exploit when specifying the IMQ. In particular, by selecting 7, 3 and e, we want to
m blk
where fux = fg + 0² V²V, y log(w?) and Jws da
diag( o2 ww 2.). Moreover, the associated predictive is
pe(yslyra-1)= - / = =N(yy;; Sw):
it an ideal choice
where p(yxlzz) = and Sn := HP HP + I
ha = Hma
See Appendix B.1 for the proof. We refer to the SSM formulation of RCGPs as ST-RCGPs, and note that the above closely relates to existing GB filtering updates (Duran- Martin et al., 2024; Reimann, 2024); but differentiates itself by using the weighted score-matching objective of Barp et al. (2019)). In contrast to what we propose in this paper, however, these previous filtering methods do not deal with hyperparameter optimisation and smoothing, or notions of down- woinhtiigh optimally.
1. wewn- observations far from where we expect the data to be, as specified by the center of the data (via 7);
2. match the rate at which we down-weight observations with how confident we are in our estimate of the center- something that was not considered for RCGPs (via c);
3. be able to recover the STGP when there are no outliers, that is in well-specified settings (via !))
5

---

Robust and Conjugate Spatio-Tempone Gaussian Processes
First, since y H+ WiMQ [xx v) peaks the centering func- tion ((xx), we let (xxx) == By By--pc((/us--)) [Yg], which is the mean of our GB filtering posterior predictive, thus re- specting principle 1 as observations far from our prediction are down-weighted. This approach removes the vulnerabil- ity of the ST-RCGP to a poor prior mean and does not incur additional cost as filtering is part of the inference procedure of the ST-RCGP, thus resolving Issue #1. We also note that this estimator of the center is robust, since we use the GB posterior predictive, which we will see is itself robust.
compute the PIF in closed form follows:
(IF( D) = KL (Pc((ID))lpe((IPmm.)
We call a posterior robust if the PIF is bounded in its first input, this implies that even as lym., - -mmil - co, the contamination's is effect on the posterior is bounded. Our next results formalises the robustness of ST-RCGP.
Proposition 3.3. For w wIMQ with B = xx) fg and c² (xx) - -iia(Ssk) the PIF of ST-RCGP satisfies
Second, we observe in Figure 2 that reducing с narrows WIMO. which in turn increases the rate V WMM x 1/c² at which observations are down-weighted To respect principle 2, we let
sup PIIF D)I sc,
Mo. ER
making it robust to outliers.
e(xx)) == Eyyrettll nl(Yy - (xx
See Appendix B.3 for the proof.
which is the variance of the GB filtering posterior predictive and reflects our confidence in This selection of e is systematic, adaptive, and relates to lyk= - -(xx)), as it is extracted from the same distribution as y is. These two features are lacking in RCGPs and resolve Issue #3.
Robust Hyperparameter Optimisation The parame- ters f we need to optimize consist of the noise level o² and kernel parameters such as lengthscale, amplitude and smoothness. Based on ply* Zzg, 9) = .N(y;; fr.. I!I..) from Section 2 and the GB posterior pc (2 O) N (m ( a standard approach would be to minimise the objective
Finally, we set / a/ v/v2, because as seen in Proposi- tion 3.1, 33 is superfluous hyperparameter when a is con- sidered, and we wish to recover the STGP in the event where ((xx)= y/ for k = 1.... n4 or equivalently w = 31, ac- complishing the purpose of principle 3.
(0) = Σ log 1ogpc(yalyna-1,.)
In the spatio-temporal setting, ? and c are vec- tors Th = (7(x),)),.... x)) and Cj "(x,,,))...,,,,.....
which has a closed form provided in Appendix C.3. This ap- proach is a version of leave-one-out cross-validation, which has been proposed in Altamirano al. (2024) to fit RCGP's hyperparameters (see Appendix C.3). Although objectives like ip are common, they are problematic for robustness. In particular, uses form of log- likeeliodood loss on the poste- rior predictive-a loss that is well-known to be susceptible to outliers. Unsurprisingly, this makes yp a poor criterion for hyperparameter optimisation-it will overfit the hyperpa- rameters to outliers, and undo most of the robustness gains we achieved through Proposition 3.3 (see Appendix C.3). Fortunately, there is simple and elegant way around this: weighted likelihood approaches (see e.g. Field & Smith, 1994; Windham, 1995; Dupuis & Morgenthaler, 2002; De- waskar et al., 2023). Specifically, we define
Here, using the posterior predic- tive PC (y for adaptivity, we take Y == f. and ci = diag(S (s.) In Figure 11 (presented in the appendix), we illustrate how our selection of hyperparameters improves the posterior. Although we favour the above specification of Yk. there are other ways to estimate the center of the data that can be appropriate in special cases (see Appendix C.1
Robustness To study the robustness of ST-RCGP to out- liers, we use the classical framework of Huber (1981). Con- sider the dataset D- {(xx, yx)] K= and suppose that for some m E (1...., ne) and j E {1, (1...., [.....,,,,, we replace single observation Vm..j in the vector ya by an arbitrar- ily large contamination ymmj. resulting in the contaminated dataset De m..j = (D) \ {Xm,, ym))) ))) U {(xm: Y }} where ym = (U.,,...., Dm ..., ). The influence of such contamination is measured using the divergence between a posterior based on D and the posterior based on Dim..
The
PPB (0)
To log wow.ocyylyy........
where we := w (Wk) E R is a function that maps the weight vector w& at time t5 to type of summary indicating whether or not the observations at time 4 jointly constitute an outlier along the temporal dimension of the observation process. In the absence of a spatial dimension, this value can be chosen as the weight itself so that Ty := Wy, as it reflects our belief that the corresponding point is an outlier.
As function of Iym,, y. this divergence is called the posterior influence function (PIF) and has been studied in Ghosh & Basu (2016); Matsubara et al. (2022); Altami- rano et al. (2023; 2024); Duran-Martin et al. (2024). For Gaussian posteriors, Altamirano et al. (2024); Duran-Marthn et al. (2024) used the Kullback-Leible (KI.) divergence to
6

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
1.0
robustness outliers and statistical efficiency- defined in this paper as model's ability recover the standard GP in well-specified settings we report the expected weighted ratio (X) == (x-)p( [w(x, "//wwrGP(x, Y)]. We expand on these metrics in Appendix C.I.
ROGP
0.8
vn(x)
0.6
0.4
0.2
1.0
ST-ROGP
0.8
Even when investigating temporal tasks without spatial di- mensions, we keep the name ST-RCGP to clarify that in- ference proceeds (i) via the state-space representation of Proposition 3.1, and (ii) via hyperparameters chosen using the methods of Section 3 3-rather than those used for vanilla RCGPs in Altamirano et al. (2024).
Contrggo
0.6
0.4
0.2
0.2 0.1
0.66
0.5
1.0
0.2 0.1 0.6 0.8 1.0
Quantiles
Figure 4. Simulated Temporal Data With Focussed Outliers. Left:
Fixing vanilla RCGP In previous sections and Figure 3, we found that vanilla RCGP is vulnerable to prior mean misspecification (Issue #1), can produce poor uncertainty estimates (lssue #2) and cannot correctly select (Issue #3). In our first experiment, we show how ST-RCGP im- proves on those issues. To this end, we simulate data from a GP and compare the fit produced by both algorithms in Figure 4. The plotted predictives show that RCGP is com- promised by outliers due to Issue #1 and Issue #3, while the coverage plots point unreliable uncertainty estimates in- duced by Issue #2. In contrast, using adaptive centering and shrinking function allows ST-RCGP to produce reliable uncertainty estimates and predictions.
NN 200 points are generated from GP with Matérn kernel (v 3/2). and contaminated with 5% outliers. RCGP uses y(x) = m(x) 1. Σ Y/ and Qx (0.95). While predictions are made on data with outliers, we fit the hyperparameters of RCGP on decontaminated data to make as strong a competitor as possible. In contrast, ST-RCGP uses the robust hyperparameter optimisation and selection from Section 3 to and predict from data with outliers. Despite this, ST-RCGT performs much better than RCGP. Right: The coverage plots measure how often the uncontaminated data falls within the predictive distribution's 's 9-th quantile.
N
In the spatio-temporal setting, however, we might instead define it as the mean, lower quantile, or the minimum of weights across spatial locations. We refer the reader to Appendix C.3 for more details and experiments showing how PGB improves on yp in temporal and spatio-temporal settings. On top of the better parameter estimates, we also find that the weighted approach reduces the variance in gradient computation (Wang et al., 2013), and seems to yield faster convergence. The computational cost of (au is ((nmg) the same as for inference with STGPs.
ST-RCGP Well-Specified Settings While robust meth- ods offer protection from contaminated data, some can only do so at the expense of statistical efficiency. Here, we show that ST-RCGP remains statistically efficient in well- specified settings and robust to outliers when required. The setup we use is described in Appendix C. 13, and results are reported in Table 1. While RMSE and NL.D are compara- ble across methods in well-specified settings, STGP suffers from a clear drop in NLPD and RMSE when outliers are introduced. In contrast, compared to other ST-RCGP maintains the lowest RMSE and NLPD in both cases, Also, its EWR is high in well-specifted settings and drops in the presence of outliers to obtain robustness. The ST-RCGP thus exhibits the properties we seek out of robust methods in well-specified settings.
4. Experiments
methods, the
In the remainder, we study the advantages of ST-RCGP on numerical examples. First, we investigate how ST-RCGP improves upon RCGP. Second, we explore ST-RCGP in well-specified settings without outliers. Third, showcase its virtues financial time series with severe outliers and its superior performance relative to competitors. Our last numerical experiment studies the robustness ST-RCGP ex- hibits towards spatio-temporal temperature anomalies. The code to reproduce all experiments is available at ht tps: / //gi thub. com/w wiw.i..... 1 www.appaaaataaaaaee
Table I. EWR, RMSE, and NLPD for different methods in well- specified settings and with data containing outliers.
Outliers
Method
EWR
NL.D
RMSE
Throughout Section 4, we evaluate experiments based root mean squared error RMSE(X.) 9) /By--u((Ixx) [YY - e and evaluations of the neg- ative log predictive distribution NLPD(X, g) == Ey-po((1x [- log Ps Vl].,6)) on the test data, where Po is the Gaussian density with mean 9 and variance 8² specified by the model. To capture the tradeoff between
STGP
1.0 + 0.0
30 + 2
0.38 + 0.02
Yes
RCGP
0.886
0.001
6.6 0.3
0.21 + 0.01
ST-RCGP
0.863
0.002
5.5 = 0.2
0.19 + 0.01
STGP
1.0
6.5 i 0.2
0.19 1 0.01
No
RCGP
0.894
0.001
5.6 0.1
0.19+0.01
ST-RCGP
0.903 + 0.001
5.7
0.19
7

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
April 22nd, 2013
April 23rd, 2013
Execution Time (s)
14700
40
DJIA Index
14650
14700
30
14600
14650
20
14550
14600
10
14500
12:57
13:12
13:26
0
14450
10
11
12
13
14
15
16
10
11
12
13
14
15
5000 10000
Time
(Hour)
Time (Hour)
# Data Points
Data
X
Crash
GP
RCGP
RP
ST-RCGP
Prior Mean
Figure 5. Comparing Estimates of DJIA Index During Twitter Flash Crash Incident. RCGP uses prior mean of the data the days. We notice a dip in the posterior predictive of RCGP when the prior mean aligns with the crash event (red), whereas ST-RCGP and relevance pursuit algorithm (RP) from Ament al. (2024) remain unaffected. In the right plot, we show computational time vs the number of data points. Note that RP has cost K times the standard GP, where typically K >> 1.
Robustness During Financial Crashes On April 23rd 2013, the Associated Press's Twitter account was hacked and posted false tweets about explosions at the White House. This led to a brief but significant market drop, including in the Dow Jones Industrial Average (DJIA), which quickly re- bounded after the tweet was proven false. This data set was previously studied by Altamirano al. (2024) and Ament et al. (2024). We plot the GP, RCGP, and ST-RCGP fits in Figure 5, as well as the robust GP via relevance pursuit (RP) fit from Ament et al. (2024). While RP behaves as desired, the plot reveals that the GP not robust to the crash. Interestingly, the RCGP performs even worse since Y, which is chosen as the constant prior mean obtained by averaging two days of data, happens to be close to the outliers during the crash, implying that RCGP is centered around the outliers. This is another instance of Issue #1, and it is addressed by ST-RCGP. By using an adaptive center- ing function introduced in Section 3, ST-RCGP is centered around more reasonable values. The right-hand panel of Figure also highlights that ST-RCGP leads substantive computational gains relative RCGP and GP due its state-space representation from Proposition 3.1. We do not plot RP because it is implemented using different package, and computational time is thus not comparable. But the cost of RPs a multiple of that for GPs (Ament et al., 2024), and hence also significantly more than that of RCGPs and ST-RGCPs. is also unclear how easily RPs could use the spatio-temporal structure get linear-in-time cost. For this reason, do not explore RP beyond this experiment.
Table 2. Performance comparison N 46800 data points with outliers. Below, Total (s) denotes clock time for full inference, 1-Step (ms) the estimated slope of a linear model fitted to execution time data for different data set sizes. Further, RMSE and NLPD use 1000 test points around the induced crash that not outliers.
Total (6)
1-Step (ms)
RMSE
NLPD
STGP
8.0 = 0.7
0.17 = 0.01
0.54 + 0.01
13.9 0.3
ST-RCGP
9.4:0.
0.20 = 0.02
0,145 + 0.002
-0.02 0.01
MEP
29.1 0.5
0.60 = 0.01
0.15 0.01
-0.77 : 0.03
0.03
MVI
2.6
0.61 = 0.01
0.15
-0.53 :
0.04 0.04
MLa
29
0.62 = 0.03
0,17 + 0.01
-0,39
Appendix C.14). Naive GP implementations struggle with data of this size, and so we focus on inherently sequential methods. Table 2 summarises the results, and compares ST-RCGP against STGP, as well as several off-the-shelf in- ference methods for sequential GPs with Student's 's errors from the BayesNewton package (Wilkinson et al., 2023) that include Markov expectation propagation (MEP), Markov variational inference (MVT), and Markov Laplace (MLa). While STGP and ST-RCGP have similar computational cost, the robustness of ST-RCGP leads to superior performance. Conversely, ST-RCGP's performance is comparable to mod- els using Student's 's errors for this problem, but its compu- tational cost is substantively lower. This is true despite the fact that ST-RCGP is exact, while the other robust methods in Table 2 only produce approximate inferences.
Since the flash crash dataset only contains N 810 data points, we further explore the computational properties of ST-RCGP by considering trading day of index futures data with N = 46, 800 data points. We then synthetically induce a crash similar to that in the previous example (see
Forecasting Temperature Across the UK Having estab- lished the performance of ST-RCGP for simpler problems, our last experiment studies its behaviour on apatio-temporal temperature data with synthetically induced outliers. In particular, we induce what are often referred to as focussed
8

---

Robust and Conjugate Spatio-Temporey Gaussian Processes
September
October
November
December
CC
30.00
25.79
Table 3. Performance on temperature data. Models perform inference on contaminated data, but RMSE and NLPD use outlier-free data. November introduces focussed outliers, and De- cember onc-step forecast.
STOP
50"N
17.37
13.16
8.95
STGP
ST-RCGP
4.74
RMSE
NLPD
RMSE
NLPD
ST-RCGP
55"N
0.53
May
0.79
1.63
0.59
0.76
Jun
0.77
1.53
0.54
0.67
50"N
33.6
Jul
0.76
1.49
0.55
0.75
Aug
0.77
1.61
0,50
0.66
10"W5"W O' STE Wo O' STE W0WWSWW O" S WO W OF 5
Sep
0.82
2.13
0.49
0.71
Oct
1.23
7.50
0.50
0.63
Figure 6. Temperature fit across the UK between September and December 2023. Focussed outliers appear Novemher and shown as red crosses. We forecast the month of December. September the last month included in parameter optimisation. include September show performance on previously seen data.
Nov
3.30
41.62
1.19
2.44
Dec 3.84
3.30
1.38
2.13
outliers, and which simulate the impact of rare natural phe- nomena that affect neighbouring weather stations. The data we use collected by the Climate Research Unit (see Harris al., 2020) and measures temperature from 16/01/2022 to 16//1/2022 at 571 locations, containing N = 11, 191 data points in total. Hyperparameter optimisation is per- formed from 16/01/2022 to 30/09/2023, with later dates serving as test data. October and November 2023, we perform in-sample predictions, and December 2023 used for one-month temperature forecast. Figure and Table display the results, and Appendix C.15 provides coverage analysis. The two models perform similarly when there are no outliers. However, the ST-RCGP outperforms the STGP slightly because the outliers impact the STGP posterior dur- ing the months before they are introduced due to smoothing. With outliers, the STGP loses prediction accuracy (NI.PD and RMSE), whereas ST-RCGP maintains consistent RMSE and NLPD over time, offering more reliable predictions at comparable computational cost.
Our method builds on classical STGPs, which remain vul- nerable to scaling in the spatial dimension and could be addressed similarly to Hamelijnck et al. (2021) through vari- ational approximations. Furthermore, computationally suboptimal if parallel computing is avail- able, in which case parallel-scan algorithms (see also Särkkä & Farefa-Fernández (2020)) could be interesting to adapt to ST-RCGPs if possible. However, we consider these inves- tigations to be outside the scope of this paper. An entirely different but equally relevant endeavour for future research is explore the use of non-Gaussian likelihoods within the RCGP framework. This would extend its use beyond the standard regression setting, and allow its application to classification problems, count data, and various bounded data domains.
method may be
Acknowledgements
WL was supported by UCL's Center for Doctoral Training in Data-Intensive Science and by the Alan Turing Insti- tute. MA was supported by a Bloomberg Data Science PhD fellowship. FXB and JK were supported by the EPSRC grant EP/Y01 and JK was additionally supported by EP/W005859/1.
5. Conclusion
We proposed ST-RCGPs based on an overhaul of the RCGP framework of Altamirano et (2024) that addressed some major drawbacks of vanilla RCGPs, further improved their computational efficiency and paved the way for their use in spatio-temporal problems. ST-RCGPs have the computa- tional complexity of STGPs. but additionally provide robust- ness to outliers. Further, we empirically demonstrated that ST-RCGPs match the performance of the relevance pursuit algorithm from Ament al. (2024) and of various robust GP methods from the Bayes/eewon package (Wilkinson et al., 2023) at fraction of their computational cost and without approximation error.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences our work, none which we feel must be specifically highlighted here.
9

---

Robust and Conjugate Spatio-Tempore Gaussian Processes
References
Dupuis, D. J. and Morgenthaler, S. Robust weighted likeli- hood estimators with application to bivariate extreme value problems. Canadian Journal of Statistics, 30(1): 17-36, 2002.
Aftab, W., Hostettler, R., De Freitas, A., Arvanch, M., and Mihaylova, L. Spatio-temporal Gaussian process models for extended and group object tracking with irregular shapes. IEEE Transactions on Vehicular Technology, 68 (3):2137-2151, 2019.
Duran-Marting G., Altamirano, M., Shestopaloff, A., Sanchez-Betancour L.,, Knoblauch, J., Jones, M., Briol, F.-X., and Murphy, K. P. Outlier-robust Kalman filtering through generalised Bayes. In International Conference on Machine Learning, pp. 12131-12171. PMLR, 2024.
Altamirano, M., Briol, F.-X., and Knoblauch, J. Robust and scalable Bayesian online changepoint detection. In International Conference on Machine Learning, pp. 642- 663. PMLR, 2023.
Duvenaud, D. Automatic model construction with Gaussian processes. PhD thesis, Apollo University of Cambridge Repository, 2014.
Altamirano, M., Briol, F.-X., and Knoblauch, J. Robust and conjugate Gaussian process regression. Interna- tional Conference on Machine Learning, pp. 1155-1185. PMLR, 2024.
Field, C. and Smith, B. Robust estimation: A weighted maximum likelihood approach. International Statistical Review/Revue Internationale de Statistique, pp. 405-424, 1994.
Ament, S., Santorella, E., Eriksson, D., Letham, B., Ba- landat, M., and Bakshy, Robust Gaussian processes via relevance pursuit. Advances in Neural Information Processing Systems (to appear), 2024.
Fonseca, T. C.,, Lobo, V. G., and Schmidt, A. M. Dynamical non-Gaussian modelling of spatial processes. Journal of the Royal Statistical Society Series C: Applied Statistics, 72(1):76-103, 2023.
Barp, A., Briol, F.-.., Duncan, A., Girolami, M., and Mackey, L. Minimum Stein discrepancy estimators. Ad- vances in Neural Information Processing Systems, 32, 2019.
Garnett, R. Bayesian optimization. Cambridge University Press, 2023.
Ghosh, A. and Basu, Robust Bayes estimation using the density power divergence. Annals of the Institute of Statistical Mathematics, 68(2):413 20166
Bauer, M.. Van der Wilk, M., and Rasmussen, C. E. Un- derstanding probabilistic sparse Gaussian process approx- imations. Advances in Neural Information Processing Systems, 29,2016 2016.
Golub, G. H. and Van Loan, C. F. Malrix computations. JHU press, 2013.
Bissiri, P. G., Holmes, C. C.,, and Walker, S. G. A general framework for updating belief distributions. Journal of the Royal Statistical Society: Series B (Statistical Method- ology), 78(5):1103-111) 2016.
Hamelijnck, O., Wilkinson, W., Loppi. N., Solin, A., and Damoulas, T. Spatio-temporal variational gaussian pro- cesses. Advances in Neural Information Processing Sys- tems, 34:23621 -2663331 2021.
Bock, C., Aubet, F.-X,, Gasthaus, J., Kan, A., Chen, M., and Callot, L. Online time series anomaly detection with state space gaussian processes. arXiv preprint arXiv:2201.0677 2022.
Harris, I., Osborn, J., Jones, P., and Lister, D. Version 4 of the cru ts monthly high-resolution gridded multivariate climate dataset. Scientific data, 7(1):109, 2020.
Hartikainen, J. and Särkkä, S. Kalman filtering and smooth- ing solutions to temporal Gaussian process regression models. In IEEE International Workshop on Machine Learning for Signal Processing, pp. 379-384, 2010.
Chen, W. Y., Barp, A., Briol, F.-X., Gorham, J., Girolami, M., Mackey, L., and Oates, C. Stein point markov chain monte carlo. International Conference on Machine Learning, pp. 1011-1021. PMLR, 2019.
Hartikainen, J., Riihimäki, J., and Särkkä, S. Sparse spatio- temporal Gaussian processes with general likelihoods. In Artificial Neural Networks and Machine Learning, pp- 193-200. Springer, 2011.
Dewaskar, M., Tosh, C., Knoblauch, J., and Dunson, D. B. Robustifying likelihoods by optimistically re- re-ieintting data. arxiv preprint arXiv:2303. 10525, 2023.
Drineas, P., Mahoney, M. W., and Cristianini, N. On the Nyström method for approximating a gram matrix for improved kernel-base learning. Journal of Machine Learning Research, 6(12), 2005.
Heaton, M. J., Katzfuss, M., Ramachandar, S., Pedings, K., Gilleland, E., Maanshardt-Shamseldin E., and Smith, R. L. teatio-temporal models for scage-scale indicators of extreme weather. Environmetrics, 22(3):294-3333 2011.
10

---

Robust and Conjugate Opatio-Temports Gaussian Processes
Matsubara, T., Knoblauch, J., Briol, F-X.,, and Oates, C. J. Hensman, J., Fusi, N., and Lawrence, N. D. Gaussian processes for big data. In Uncertainty in Artificial Intelli- gence, pp. 282, 2013.
Generalized bayesian inference for discrete intractable likelihood. Journal of the American Statistical Associa- tion, 119(547):23452255
2024.
Huber, P. J. Robust statisties. Wiley Series in Probability and Mathematical Statistics, 1981.
Nickisch, H., Solin, A., and Grigorevskiy, A. State space Gaussian processes with non-Gaussian likelihood. In International Conference on Machine Learning, Proceed- ings Machine Learning Research, pp. 3789-3798, 2018.
Hyun, W., Li, Y., Huang, C., Styner, M., Lin, W. Zhu, H., Initiative, A. D. N., et al. Stgp: Spatio-temporal Gaussian process models for longitudinal neuroimaging data. Neuroimage, 134:550-563 562500662,10 2016.
Pleiss, G., Gardner, J., Weinberger, K., and Wilson, A. G. Constant time predictive distributions for Gaussian pro- cesses. In International Conference on Machine Learning, pp. 4114 4123. PMLR, 2018.
Hyvärinen, A. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(24):695-709, 2005.
Rauch, H. E., Tung, F., and Striebel, C. T. Maximum likeli- hood estimates of linear dynamic systems. AIAA Journal, 10001445-145000 1965.
Jewson, J., Smith, J. Q., and Holmes, C. Principles of Bayesian inference using general divergence criteria. En- tropy, 20(6):442, 2018.
Reece, S. and Roberts, S. An introduction to Gaussian processes for the Kalman filter expert. In International Conference on Information Fusion, pp. 1-9, 2010.
Kalman, R. E. A new approach to linear filtering and pre- diction problems. Transactions of the ASME-Joural of Basic Engineering. 82(Series D):35-45, 1960.
Reimann, H. Towards robust inference for Bayesian filtering of linear Gaussian dynamical systems subject to additive change. Master's thesis, Universität Potsdam, 2024.
Key, O., Fernandez, T., Gretton, A., and Briol, F.-X. Com- posite goodness-of-fit tests with kemelss. In NewrlPS 2021 Workshop Your Model Wrong: Robustness and Misspecification Probabilistic Modeling, 2021.
Riabiz, M., Chen, W. Y., Cockayne, J., Swietach, P., Niederer, A., Mackey, L., and Oates, C. J. Optimal thinning of meme output. Journal of the Royal Statis- tical Society Series B: Statistical Methodology, 84(4): 1059-1084 2022.
Knoblauch, J., Jewson, J., and Damoulas, T. An optimization-centric view on Bayes' rule: Reviewing and generalizing variational inference. Journal of Machine Learning Research, 23(132):1-109 2022.
Ruanaidh, J. O. and Fitzgerald, W. Numerical Bayesian methods applied to signal processing. Springer Science & Business Media, 1996.
Krige, G. A statistical approach to some basic mine valuation problems on the Witwatersrand. Journal of the Southern African Institute of Mining and Metallurgy, 52 (6):119-139, 1951.
Santner, T. J., Williams, B.J., and Notz, W. I. The Design and Analysis of Computer Experiments. Springer, 2nd edition, 2018.
Kupilik, M. and Witmer, F. Spatio-temporal violent event prediction using Gaussian process regression. Journal of Computational Social Science, ((2):437-45) 2,,,,451,1111 2018.
Särkkä, S. and García-Fernández, Á. F. Temporal paral- lelization of bayesian smoothers. IEEE Transactions on Automatic Control, 66(1):299-306, 2020.
Liu, X. and Briol, F.-X. On the robustness of kernel goodness-of-fit tests. arXiv:2408.05854 2024.
Sarkka, S. and Hartikainen, J. Infiniteddimensinaal Kalman filtering approach to spatio- spatio-temporal Gaussian process regression. In International Conference on Artificial Intel- ligence and Statistics, volume 22, pp. 993-1001, PMLR, 2012.
Lyu, S. Interpretation and generalization of score matching. In Proceedings of the Twenty-Fifth Conference Uncer- tainty in Artificial Intelligence, UAI '09, pp. 359-366, Arlington, Virginia, USA, 2009. AUAI Press. ISBN 9780974903958.
Särkkä, S. and Solin, A. Applied stochastic differential equations, volume 10. Cambridge University Press, 2019.
Matsubara, T., Knoblauch, J., Briol, F.-X., and Oates, C. Robust generalised Bayesian inference for intractable likelihoods, Journal of the Royal Statistical Society Series B: Statistical Methodology, 8402)997-1002, 2022.
Senanayake, R., O' O'Callaghan, S., and Ramos, F. Predict- ing spatio-temporal propagation of seasonal influenza using variational Gaussian process regression. In AAAI Conference on Artificial Intelligence, volume 30, 2016.

---

Robust and Conjugate Tempoo-e Gaussian Processes
Sinaga, M. A., Martinelli, J., and Kaski, S. Computation- robust Gaussian processes. In NeurIPS 2024 Workshop Bayesian mecision-making and Uncertainty, 2024.
Windham, M. Robustifying model fitting. Journal of the Royal Statistical Society. Series B (Methodological), pp- 599-609, 1995.
Wu, P.-S. and Martin, comparison of learning rate selection methods in generalized Bayesian inference. Bayesian Analysis, 18(1):105-132, 2023.
Solin, A. Stochastic differential equation methods for spatio- temporal Gaussian process regression. PhD thesis, Aalto University. 2016.
Xu, J., Scealy, J. L., Wood, A. T., and Zou, T. Gener- alized score matching for regression. arXiv preprint arXiv:2203. 2022.
Solin, A. and Särkkä, S. Gaussian quadratures for state space approximation of scale mixtures of squared ex- ponential covariance functions. In IEEE International Workshop on Machine Learning for Signal Processing, pp. 1-6. IEEE, 2014.
Yu, S., Drton, M., and Shojaic, A. Generalized score match- ing for non-negative data. Journal of Machine Learning Research, 20:1-70, 2019.
Stein, M. Interpolation of Spatial Data Some Theory for Kriging. Springer Science+Business Media, 1999.
Särkkä, Bayesian Filtering and Smoothing. Institute of Mathematical Statistics Textbooks. Cambridge University Press, 2013.
Tadayon, and Rasekh, A. Non-Gaussian covariate- dependent spatial measurement error model for analyzing big spatial data. Journal of Agricultural, Biological and Environmental Statistics, 24:49-72, 2019.
Titsias, M. Variational learning of inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics, pp. 567-574. PMLR, 2009.
Wang, C., Chen, X., Smola, A., and Xing, F. P. Variance reduction for stochastic gradient optimization. In Pro- ceedings of the 26th International Conference on Neural Information Processing Systems-Volume 1, pp. 181-189, 2013.
Waxman, D. and Djurić, P. M. A gaussian process-basee based streaming algorithm prediction of time series with regimes and outliers. In 2024 27th International Confer- Information Fusion (FUSION), pp. 1-8. IEEE, 2024.
Wilkinson, W. J., Särkkä, S., and Solin, A. Bayes-Newton methods for approximate Bayesian inference with psd guarantees. Journal of Machine Learning Reseurch, 24 (83):1-50,000
Williams, C. K. and Rasmussen, C. E. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
Wilson, A. and Nickisch, H. Kernel interpolation for scal- able structured Gaussian processes (kiss-gp). In Interna- tional conference on machine learning, pp. 1775-1784. PMLR, 2015.
12

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
Supplementary Material
The Appendix is structured as follows: We first introduce notation. Then, in Appendix A, we derive the loss function used in the main paper. Next, in Appendix B, we provide the proofs of all our theoretical results. Finally, in Appendix C, we provide additional details on our numerical experiments, as well as further experiments that complement those in the main text.
Notation
Suppose we have vector (V1, U2, UN) R. where vs RR R, and matrix M RNN, Then,
U1
diag(v) :=
and diag(M): : (M11, M22..... MNN) where My == (M)
UN-1
UN
v² = (uf..... vz)
log(v): = (log(11).. log(uN)
0 T uN By
V V
By
0
U1.....
where 0 By 04 is the partial derivative of o with respect to y.
w W1, U T for w ERN.
We now remind the reader of the dimensionality of matrices and vectors crucial to the main algorithm (Proposition 3.1) of this paper:
Ag-1, Eg-1 (=((ot)))nneettt)) are the transition matrices from Equation (4)
H (n,(x)) is the measurement matrix from Equation (4)
Yk: WK, PGB E th are the observations, weights and filtering predictives from Proposition 3.1.
PGB E (1)+ the covariance matrix from Proposition 3.1
&
m kij
E (n.(v+1) the mean from Proposition
KGB E (=(+)))),, is the Kalman gain matrix from Proposition 3.1
Jw, Km. is the weight-based matrix from Proposition 3.
Finally, when we use p(-) = .V ; E), we are referring to the density the normal distribution with mean f and covariance and when we use (H, E), we are referring the distribution.
A. Score-Matching & Loss Function
This section provides the derivation of the loss function used in the main paper based on the weighted Fisher divergence.
We define the weighted Fisher divergence at a fixed time tę- Let f.) = f(., tx) the spatial model at time t.- Let P& (y|s) ply|s, ss tx) be the density the true date-gengg process at time tg. and Pıx (yls) = f())f((s) be the density of our model at time 4g.
13

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
The weighted Fisher divergence for regression between the model and data-generating process time l2 depends on the corresponding scores ff, (y]s) Vy log (y)s) and 8h (ylss Vy log px(ss)) and is given by (Barp et al., 2019; Xu al., 2022; Altamirano et al., 2024):
D == Es-pE.
By (y
(sx (sx - S/4 (Y [
+C ESupa,
Ey [x>pp((1s) (.15)
2V () (Y (X (Y
where Pk,. p())tt) is the marginal, (x, (s, tx), y). and the qquality-which crucially does not depend on sp anymore- holds up to an additive constant C' not depending on fk-
Now, consider dataset such that X (x Xk,n,) )T for xk.j = (sj, (x) E X T. yv = (yks,,,.... yk,, ) and
let fg = f(xx) = (f( (xk.,.)) Morcover, let zy be related to as Hz The empirical loss we obtain for filtering is then
n, | 12 2V, ((wxg)²s
L(x:: Ski zR)
(11)
1=1
where Wk.j = wo(x Xk.j Dk.j), and 8s..j 85a (Un.))-
B. Proofs of Theoretical Results
B.1. Proof of Proposition 3.1
In the following, then show that the weighted score matching loss with Gaussian model yields quadratic loss. Finally, provide the derivation of the GB predictive.
derive the generalised Bayes filtering posterior distribution when the loss function is quadratic. We
Let us assume a quadratic loss in 24 or equivalently in f. == Hzg of the form:
1
L(x2, yk, 2x.)
(G. R 1f.
fr vR Ck),
2n
where Rs E RR is an invertible matrix, E "h and Ck E RR. The GB filtering update equations are then
pe(zaly1a) x p(zxly 1:k E Zk,)
1
P R (zz majk- exp (-1 (fr R fr - 2f vo.)
exp
(2
-p(k-1)
2
(12)
1
(P-1
M + H R, 1 H)zz 2n (P Mjk- - + H vz.) v=)))
exp
2
GB which implies that the mean & and covariance PGB of the GB posterior Pe
PGB A/K
(P-)
H R H
k|K
(13)
GB
PGB kijk
(P-
majk-1 + HVv))
m
k/k
klk-1
As with the typical Kalman filter, those equations can written in the form
PGB
PP- kjk-1 + H HRR," 'H)
KGGH
PGB H R
(14)
GB k/
majk-1 + K GB
RR Hmaj--),
where KCB is the Kalman gain matrix. The typical Kalman filter equations those used for STGPs are recovered when
Va R yx and Re == -
14

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
Weighted score matching and Gaussian likelihood. We now show that the loss function defined in Appendix A and Gaussian likelihood lead to quadratic loss in fg. We have Gaussian likelihood, which gives score function of the form:
2 (-1 To -,, (yk= fi:)) 8f.k = Vy (olog((yale)) = (fe- yk) o In, (15)
p(yn|zexXe) Xx) & exp
where fg = Hz 2x
1.0.11.18811 (phknn)) Then, the loss is given (2k,11-.n Fkones 1T yk (1 -.. v. J and sj
by
C(xx:Y yk: 2z)
wkk (Kk.s + vV (wRt))).ts)
(Jk. Uk.s) -huo) d²
10 (Jk. - yk.i)
Wk.j
+ +2V y
o²
j=1
2 fh.
1 whi + 0² vu(
2
2/
2 Vo(uh
T1.s
Wk
kkg
o²
Now let us group all the elements that do not depend on f and call C(yk.)
2 a² vo(wi,tts)) C(vn.)
L(x;, Yk, 2x)
W 1h.
g Sk.ykk. wkk + G2 Wkyskk Vywk. + ((nss)
11g
R2
04 T (yk.s 20 (wx,)) "Vyuks fk.s + C(yes)
W
04
j=1
Next, we rewrite the loss in terms of vectors and matrices as follows:
((x 1 2
2n, 04 fg diag(w af fr (diag(wī)yx 20 fo'ddag(w,, (W))) + C(va)))
where C(ys)) = - C(y.s), diag(w?) E Ru,xn, is the diagonal matrix of the vector w We O We for O the element-wise multiplication operator, and V, wk = (V WK,11.1.... V Wy.kk..) This leads to
R (y:: = diag(w?)
04
(16)
2 va(y::A)) = (diag(wi)) (wi))+ 20 dioddag (wx) V,)),,,))
GP predictive The GB predictive be written as:
Pe pe(yslyia-1) f pww.llenpppeeaaaliaaaaa
where p(yslz*) = -yy; Hzz II. ) is the likelihood and Pz = NZ majk-1, II)--) is the predict step defined in Proposition 1.1. Since both densities are Gaussian, this integral is known and the solution is also Gaussian of the form:
Pe(yxlyinn-1) / N(yxiHHu,
Hzn, I'In II,,,,,,( (Z m P PRIK- 1dzz -N(x;; Hmyj-1. HP PPPjj-1l .H. + I'In.)
15

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
B.2. Proof that ST-RCGP reproduces RCGP for Proposition 3.2
We first define the necessary quantities. Let the posterior density of RCGP be
o p(fix (etp://eeppttn...ccoo/aaaaa fix)),
where p(f (1ssk) is the prior density, and
CROOP(XLRK MIck fick))
wo (x(x) (y 2V 2V, v (w (Xi.) Vii)), (1),)))))
(17)
kn
where indexes time, j indexes spatial coordinates, and k E {1,..., me). Note that we now have dependence on fą instead of Zy for simplicity of notation; however the two definitions are equivalent since fg- = Hzz.
For fixed hyperparameters w (for example, c, 3 for Wisq), the loss from Proposition 3.1 relates to CRCUP as
CRCOP Vi.k. fi:k)
Σ v.,
(18)
which implies that the loss LRCGP is summable (can be broken down into summation over 1,..., k). Moreover, let the posterior filtering distribution of ST-RCGP on fy be
PSTRRGGP (fg X1:ks Y1:k) PSTRRCPP (f Xi.ks V1:k-1) cxp-a--)))ppp_____-11 fix)
PST-RCGP(fL XILs VIX-1) / p
where the loss C is specified as Proposition 3.1.
The above definitions such that the prior density is the same for PRCGP and PST-ROGP, and the loss L. which is summable, is specified identically. If we further impose on the Gaussian process to be Markovian (the so-called Markss process), which is requirement for the Gaussian process to be expressed as state-space model, these assumptions will allow us to postulate the following Proposition 3.2:
Proposition B.1. Suppose that
(i) the prior distribution p(fss is identical for PReGr and PST-RCGP;
(ii) the loss function C is defined in Proposition 3.1 and LaCGr (X Yicks fack) as in Equation (17), with weights that do not depend past observations;
(iii) the GP prior f ~ Nmm, R) is Gauss-Markov process and can be expressed as state-space model.
Then,
1. (fa Y1a) = (fT-RcGr(fs Vi.k) (fittering posteriors are equal)
2.
= PST-RCGG (fa (smoothing solutions are equal).
Proof. The proof divided into two steps: the first tackle claim 1., which is that the filtering posteriors are equal, and the second to claim that the smoothing distributions are equal (claim 2.). In both cases, proceed inductively, with base case and inductive step.
We first want to show that PRCGP(fi Xi.k. PST-RCGP (fal Vic)..
Step 1: PROGP (f& (ickk y1ck) = PST-ROGP (fr X1:ks X1.k) (claim 1.).
Base Case:
16

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
Consider 1. Then from the RCGP equations in Altamirano al. (2024) (Proposition 3.1) and in our Proposition 3.1, the RCGP posterior and the ST-RCGP posterior are:
PRCGP (fi v1) = (f1; |RCGP, MRCGP)
PST-RCGP (f (f1:///STRCCGPP EST-ROGP);
where
PRCGP mi Ki (Kı oJJ "2 - (yl mwn
(19)
ΣRCGP K1 (K₃ + o²Jw,, w.)) -1 o²Jun
and
PST-RCGP Hm + HP H (o -1
y1 funs
1|1
(20)
EST-RCGP H ((P98) -1
1
H σ
om
H
H
Then, expand εROGP follows:
εRGPP = Ki (K1 + 0² JJ)n)))
-1 =²J
-1 Ki + (o2J ww:) )
K1 - K₁ (o²Jw.)
K1 K₁ (K1 o2 J11))
K1,
where in the first line, apply the Woodbury identity (Chapter 2.1.3 Golub & Van Loan, 2013), and in the second, use the fact that for two invertible matrices A, B, we have (A-1 + B B-1)-1 A(A + B)-1 B. Moreover, since fū. = Hz
EST-RCOP = HPGB H so that 1p1
ESTRCOP H
100
+ H
JJ-H
H
H
= HHGB H 1|0
HPG8 1
H
HHPB 110 H + )oJw.)
HP
where again use the Woodbury identity. But, both matrices HPGB H and K₁ represent the covariance matrix of p(f x1), and thus must be equal by the assumption of identical prior (assumption (i)). Therefore, substituting K = HPGB H in cither of the expression for ΣROGP or ESTRRGP yields MROGP = ESTRCGP.
H
1
Now, we have
PRCGP = m, + K1 (Ki + )J),,) (y1 mw = mn + ERCGP (o²Jw.) (y1 m
where in the second equality, we substitute MROGP from Equation (19). Since HP, H is the covariance of PfT- yI), then, HT EST-ROGP = εRGGp (last equality holds from previous step). Moreover, Hmjjo fi and mj both the mean of pff)). Since the priors are identical, we then have Hmjjo = fi m, and by extension m.; = fu. since the two distributions are defined with the same loss function (and thus have identical weights). Then, as before, substituting these quantities yields HRCGP = USTRRGPP concluding the base case.
Inductive Step:
Suppose that PST.RGGP (f
for some 2. Then,
PROGP (fcl (folx14ss111)
(21)
α exp(-_(xxx Vick. f(se ((
17

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
(1:sk) But ((())) poffeffi-- x x1.+) p(f, (1 111s)) by the Markov assumption f (assumption (iii).. That is, fg is conditionally independent of fi:k-2 given f1-1-
Also, exp(-f(x1xx>) Vick f(a)) = exp(-f(xx. fx)) exp(- (000--((x-----00 11 118-11, since the loss is summable. Therefore.
(R V1.k) (xxx (i)) / p(fa lh xxx)) 111: (xx_)e)exxxxxxxss--iiae
fisk- 1
where we substitute p(fi:k x1:)) exp(-[(xx, yui - f(.x-)), Now, using the definition of the PRCGP for k 1:
p(f 1 x x1x))((f (f::))(ff--xxxii)) and
y VIsks fick))
PROGP X & exp(= /
11)
Splitting the integral into f1 and f obtain:
PROGP K1sk Vicke & exp(- yx. (e)) / p(fx/fe-1,,0 x111) / fin- y
exp((-((x..,, fc)) /
where in the last equality we integrate out f1 (which integrates since PROGP a density). Finally, using the inductive step assumption and the definition of the density PST.RGGP we have:
PRCUr Nicks V1.k) x exp(-I(xx, Yk, fx)) /
PSTERGOP (f: Vi:k),
Therefore, PRCGP(fa X1.ks V1:k) N1:k) Pf V1.k. However, since both sides densities. the proportionality implies equality; that is, PROGe (fy, X1.6, V1.k) = PST-RCGP (fg X1 V1.k), which concludes the first step of the proof. Now, we want to show that the smoothing solutions are equal (claim 2.). We use similar approach; however, the proof starts from the largest possible value of k: and then goes down in values.
Step 2: We want to show that PRCGP (fa Xum. ylm,) = PcC (fx X V1mm)).
Base Case:
Consider & ne- Then, the smoothing and filtering distributions are identically defined. Therefore, by Step 1, PRCGP (Ximes Iome)) PST:RCGP(Pn, Ximnes lome).
Inductive Step:
Suppose that PSTRCGP((G++ | Vimn)) = PROGP (fa++ X for k < ne - 1. We want to show that PSTROOP (fg X Vi:m.) PROGP(ff ( X1mee Vi.n.). Then, by Theorem 8.1 of Särkkä (2013), which requires that fş be independent of yk+1m, given f (satisfied by Markov assumption (iii),
/
(fxe:i/fFPPRcCFF (fe X Vi:ne) PROGP (Cr++ X1:441 Y1:k)
PROGP(Fe X1m, Vim.) PRCGP (x1I) Vis))
dfx+1
(22)
However, by Step 1, PRCGP (fa X1skk V1sk) PRCGP(fx+11 (X Vim,)) ) PST.RGP (fa XX ).m..) PST-RCGP (fo41 (1.641..1a Y14). But
PST-RCOP (f (fglx14 and the inductive step implies Therefore, there remains to show that PROGP (fx+ (x
/ p(fo+ (wc.p//CCGFf(eee)/ppee dff
PRCGP(fe+1 V1:k)
(23a)
/ p(fe++ (PSTR))SSSPGGPG ( Yos df,
(23b)
PSTTRGGG Co 1 Visk),
(23c)
18

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
where in Equation (23a) we integrate and expand PRCGP (f&, fi+1 Xisk+1s Y1:k), in Equation (23b), apply the result from Step 1, and in Equation (23c), the definition of PSTRRGP- We conclude that
Therefore, substituting the above quantities,
/
(fle+i(ff.)ssTRCCCO (fx++ Mximm dfs+1 PST-RCOP (fx+1 x1.41111.s))
PRCUP (fgg X1:m, 3aimmyyinm V1:m.) PST-RCOP X1.ks Y1.k)
PST-RRGPP X1mm, M ).
This completes our inductive step and the proof.
B.3. Proof Proposition 3.3
To prove the robustness of ST-RCGP, we rely on the fact that ST-RCGP and RCGP share the same distribution for spatio- temporal data; therefore, we could use the following result from Altamirano et (2024) adapted to the spatio-temporal setting:
Proposition B.2 (Altamirano et al. (2024)). Suppose f ~ GP(m, k). (cI,..,ENN) (N) NO, I and let C2. E R; k: 1,2 be constants independent of j. Then, for RCGP regression with supxext.e) w w(x, y) < co has the PIF
PI D) < Ci(w(Xxa) vm 3 C2.
(24)
Thus, vx X' supyey lg w(x 3) < oc RCGP is robust since supye, Me ey PIF RCGP (Vfm. D)I < co.
Hence, it suffices to verify that the proposed weight function satisfies the necessary conditions for robustness presented in Proposition B.2.
The weight function and the hyperparameter recommended are:
wiMo(X.9) (1 1+ (u (x-)(x))² - c(x)2
(25)
with 8 = 2- 7x) = Epc [v], and 2x()) = Epc [1y - (x))²
It is straightforward to verify that WIMQ(x, v) B for all E X' and y V. Since B x² +00, it follows that SUPxxXX w(x, y) < -00 Thus, with the recommended hyperparameters, satisfies the first condition.
Now, we need to check the second condition, which is that Vx E X, supyey lul w(x, y) < too. To show this, let us consider an arbitrary x E A' and two cases for y::
Case 1: |vl S [x)) [e(x)] Since wimo xx y) S 3 for all E and y E V, it follows:
lyl (xx y) < 10| < (>x(() + lc(x))
Which implies that, as long as [xxx) < oo and [c(x)] oo we have that:
|y| - w(x, v)² < xc
sup yEy s.t.lvs (ex)))
Case 2: 1ol > hx(x] + lc(x)]
c(x)²
c(x)²
|el - xxxxx))² S 1y1
S 1v
1w/3²
= -2
[x-' x(x)²
(x-)()))² o(x)²
(y - (x))²
ly (1 x(xx)
Now, we observe that this function is decreasing for |yl > [x(x) + lc(x)], and particularly:
c(x)2
lim +00
0,
2(x)
|
19

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
since [y] [x)) Therefore, attains its maximum at [y] [x(x)] |c(x)], which leads the following bound:
c(x)²
lul wx(x, S B²
S ( [c(x))).
lvl 2(x)
Which implics that, as long as [x(x) < oo and |c(x)] < oo we have that:
sup veyy
|y/ -l,w((x))))) y) <
..c.lp/>>/00
Finally, let check that [>(x] < oo and [xx) < xo. Since 7(x) = Pe [v] and (²() == Ep [(V {] are the mean
and the variance of pc (yí x) respectively, and pc (y) is Gaussian, we know that the 1>(x) < oo and le(x)] < co.
Putting all together, we have that:
Vo E X, sup lul w(xx u)² oo.
yey
C. Additional Numerical Experiments
C.1. Performance Metries
Let po(-lx) x) be the density of the true data-generating process on the spatio-temporal grid X = (xi,...,xs) xxs) E EX" The first performance metric we use is the root mean squared error (RMSE):
RMSE(X,0) BY x) [(Y - 0)}
(Jk.j
n4 k=1j=1
where N n,n, is the number of data points, y/ is the data, and Úi the model's prediction. The second performance metric use is the negative log predictive distribution (NLPD):
NLPD(X, (. == -oo X) [- log (Y)0, 0²)
Σ log 0s o
N
where ou is the model's variance on the prediction Di, and is the Gaussian density, Finally, for the methods with weights such as RCGP and ST-RCGP, we introduce the expected weight ratio (EWR):
(xx, Y) WSTGP (x, Y)
Bi
(Uk.j 7(x,))) c(xx,))
BWR(X) (px()
where WSTGP := 2/v2 is the constant weight for the standard spatio-temporal GP. Note that by construction, EWR < 1. In particular, if Wa = WSTGP for & = =1, 1,..., n4 then EWR = and we recover the STGP posterior exactly. Therefore, since the STGP not robust to outliers, EWRs that are near one are not necessarily optimal, since an EWR of 1 implies a solution-the vanilla STGP-whicl model's ability to recover the STGP well-specified settings-and robustness outliers. In well-specified settings, then want EWR be larger and closer to When there outliers, wish the opposite so that benefit from robustness.
not robust. This metric thus conveys a tradeoff between statistical efficiency-a
C.2. Implementation Details
All experiments are run on the CPU of 2020 13-inch MacBook Pro with MI chip and 8GB of memory.
Definition of Matrices AR-1 and Σe-1 In the following, provide few examples of the SDE matrices needed compute Ag-1 and Eg-1, and explain how they obtained in practice. We define A exp(FAtt and
Al F(At 1 F) dr,
Σ
Σ
(26)
20

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
where F- n, Fe (= L n, L (=*(v+)) and Qc Kg 2 Qe [R for (Ky) Ki (8t, 8j). The notation & refers the Kronecker product of matrices, and the term exp(FAtt) corresponds the matrix exponential. From Equation (26), see that the matrices Ak-1 and - k = 1,..,, are derived from F,, L, and Qc in the SDE formulation.
SDE Matrices for Common Kernels We provide in Table 4 common kernels and their SDE matrices. That is, for each kernel we select, we specify the matrix form of Fr, L,, along with Q We later explain how these can be used to compute A and 1 The parameters in Table 4 are the lengthscale E, the amplitude ø. and the period length wo: the inputs of the kernels are t, t' IR and == t t'; the functions used are the Gamma function I' and K₂ the modified Bessel function of the second kind.
Table 4. Table of kernels and their corresponding SDE matrices
Kernel
Formula
Parameters
F,
L
Qu
Wiener Pracess
Kwplt, (²) == o² min(t,
OR
Exponential
Kesp(T) oº2 exp(-F)
Ju.
-1/8
1
20
0
Matérn 1/2 kernel
KMex. (T) o2 Fooj
K2
Ou E; X v3/l
12902
-2-
2A
1
0
Matérn v 3/2 kernel
(same as above)
on, =- v5/8
0
16o
-13
-32
33
Periodic
Kipestodic (T) = of exp
2sin 8
Σ/
nI,
Waj
The matrix A can be straightforwardly computed most programming languages that offer linear algebra computations.
The covariance matrix E is rarely computed directly; instead, it is obtained by Matrix Fraction Decomposition (MFD) (see Chapter 6 of Särkkä & Solin (2019) for an overview of the method). Provided that the matrix F, is Hurwitz, that is, all its eigenvalues have strictly negative real parts, the procedure to compute E can be further simplified (Särkkä & Solin, 2019) (i.c., no need for MFD) to
Ego = Σ1,30
E E,,00 1
where the initial covariance Σgo can be found via the steady-state solution by solving for Et in the following continuous Lyapunov equation (Särkkä & Solin, 2019; Hartikainen & Särkkä, 2010)
F F, + F. + Q., =0
Note that we provide the full implementation for the ST-RCGP and the STGP in our code.
C.3. Robust Hyperparameter Optimisation
Issue With RCGP Although the RCGP method is robust in inference, it still has well-known problems with hyperparameter optimisation when there are outliers. The method used is outlined in Altamirano et al. (2024), which corresponds to the leave-one-out cross-validation (LOO--C). The optimisation objective is posed as follows:
à², == arg max 0².8
Σ log p bogpp"(,x,,,xx.. (uIxx,y-,
where y-; =(y1.... ... 94-1, 3:41. vn). Moreover, the pseudo marginal likelihood is given (" N( R a. + 02), where
fl = zi + m; - [K + -1 (y - mw)) [K o²Jw)
-17
0, R = [K + =JJw) ) 2 w(xi,
21

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
where the notation follows that of Equation (8).
In Figure 7, simulate data with outliers to highlight the issue. We generate 80 data points from GP with Squared Exponential kernel of amplitude position four outliers around the mark. To fit the data, we follow the suggested approach from Altamirano et al. (2024) and use a constant mean function equal to the arithmetic average of the data. In Figure 7, we show the results. Clearly, the RCGP fitting the data with outliers fails to produce adequate hyperparameters, whereas the fit on outlier-free data does well.
0.2, lengthscale e = 1, and variance x² = 0.01 in normally distributed data. We
Data
Outliers
RCGP fit (Outliers)
ROGP fit (No Outliers)
-1
2
0 x
2
Figure 7. RCGP Purumeter Optimisation With and Without Outliers. The fit with outliers yields kernel lengthscale of 0.19, a kernel amplitude of 0.99, and variance of 0.026. The fit without outliers yields kernel lengthscale of 1.46, kernel amplitude of 36.9, and variance of 0.0066.
Issue with STGPs Hyperparameter Optimisation the optimisation objective for STGPs is typically
(0) = log
log [2mS (0)s,
where S₁ (0) :0 II,, + HP k/k 1 ( and Ex y= (0).
However, with outliers, this objective does not perform well since reasonable estimates of the latent function will have large εş and thus small p(yk N1:k-1, effeectvivey fitting the outliers. This makes it so & may have global minimum #* vastly different from the one in well-specified settings, in which in Figure what happens when fit STGP with a yp objective. On the left side of the plot, we fit and make predictions with contaminated data. On the right side of the plot, we fit decontaminated data and make predictions on data with outliers. Both approaches yield undesired results, but most importantly, they are vastly different-a result that stems from introducing outliers in the data used for hyperparameter optimisation. Overall, Figure shows that şo is strongly affected by the presence of outliers and is not robust objective for finding hyperparameters.
could not recover the true latent function. We show
Improving ST-RCGP's Hyperparameter Optimisation: Temporal Setting For temporal data, we now show that when there are outliers the training data and we are using the ST-RCGP algorithm, choosing PGB from Section as an objective function for hyperparameter optimisation more reliable than using the regular yp It is given by:
PGB (0)
(log (2) [J) (0))
where definitions are as above with 1, and the weights wi are the ones from the ST-RCGP. In Figure 9, we fit the ST-RCGP on contaminated data with both objective functions to the best of our ability (that is, we adapt the learning rate and the number of optimisation steps as best as we can to get optimal results). We observe that using PGB improves the optimisation process drastically compared to using s, since way fewer steps were necessary, and doing so allows obtaining reasonable hyperparameter values for inference.
22

---

Robust and Conjugate
Gaussian Processes
0-
0-
-5
STGP
-5-
STGP
-10
Data
Data
Outliers
Outliers
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
(a) Fitting STGP on Outlier Data. The hyperparameters ob- tained are: Lengthscale 8 = 0.118, amplitude dx = 3.438, observation noise 0² 2.842.
(b) Fitting STGP on Outlier-Free Data. The hyperparameters obtained are: Lengthscale = 0.236, amplitude dx 3.187, observation noise o² 0.055.
Figure 8. Impact of outliers on hyperparameter optimisation with şo. The data is generated in the same way as for Figure but with 100 data points instead and five outliers from a (0, 10). Both models use Matérn 3/2 kernel.
0-
-5
ST-ROGP
Data
-5
ST-RCGP
-10-
Outliers
-10
Data
Outliers
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
(a) Fitting ST-RCGP on Oullirr Data With yp. The hyperpa- rameters obtained are: Lengthscale f = 0.683, amplitude ax = 2.078, observation noise o² = 1.909. The optimisation process takes roughly 100 steps rate of 0.009 and is unsta- ble (does converge and produces infinite values if learning rate increased).
(b) Fitting ST-RCGP on Outlier Data With PGB. The hyper- parameters obtained are: Lengthscale { = 0.112, amplitude 4.162, observation noise 0² = 0.094. The optimisation process takes roughly 25 steps rate of 0.2 and is consider- ably more stable.
Figure 9. Impact of objective function ST-RCGP's hyperparumeter optimisation. The data is generated the way for Figure 8. Both models use a Matérn 3/2 kernel and centering and shrinking function as specified in Section 3.
Improving ST-RCGP's Hyperparameter Optimisation: Temporal Setting We turn our attention to spatio- temporal data and make the same point as in the previous paragraph, which is that choosing FGB is more reliable than using the standard sp. We use the following objective:
nee
PGB (0)
wh (log (0) + (I]
where now, the representation ty (Wx) of the weights wa time step tu given by:
kk. (8)
Ük
Qi,n, (8)
where Qu.r, (6) the ō-quantile of the weights Wg. this experiment, we choose 8 0.05. The data we use identical to that for Appendix C.4, except that the outliers only introduced at steps h = 2 and 6, and the rest remains uncontaminated. As performance metric, we use the cumulative mean absolute difference (CMAD) between our estimate of the latent function and the true latent function (not the observations), which is given by CMAD E Σ n. j=1 [fE. fr.
We the data with an ST-RCGP that has a Matérn 3/2 kernel, and a centering and shrinking function as specified in Section 3. We compare using & and POB and optimise for 30 training steps with an Adam optimiser with a learning rate 0.1. We run the process five times with newly sampled data and outliers, keeping the latent function and the distributions from which the
23

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
data and outliers are sampled the same. We find that yp yields CMAD values of [1.9974, 1.6703, 1.9517, 1.6100, 1.7321]. whereas PCB produces CMAD values of [0.4633, 0.5385, 0.5058, 0.4479, 0.5321 The yPa objective yields substantially smaller CMAD values, and thus provides more reliable hyperparameter objective function. We illustrate this in Figure 10, where we show three time steps of the first run (out of the five) for both objective functions.
8Pn
yp
Prediction
Data
Latent Fanetion
Prediction
Duta
Latent Function
1=1227
1=0 27
t-033
1=3033 0-
0-
0-
-1
1-
t-0.40
1-0.40
0-
0-
0-
1
(a) Fitting ST-RCGP on Outlier Data with the Robust spoa.
(b) Fitting ST-RCGP on Outlier Data with the Regular sp
Figure 10. ST-RCGP Fits Using The Objective Functions yoR and
C.4. Synthetic Spatio-temporal Problem in Figure
grid of size n, 25 25 is produced and repeated through n4 = 10 time steps between to 0.2 and ty 0.8. We use latent function f(s $2, = sin(2rt)s} cos(2Rt)s and generate additive noise from a aN(0,, o²), where we set 0.2. At any time step, if data point is located in the region 81 0, with 10% probability, contaminate the data point with an outlier sampled from U[[-&, 6) J [6, 8)).
For both the STGP and ST-RCGP, we fit the data with the optimisation objective and use de-contaminated data (original data without outliers) for the objective. We use the Adam optimiser with 20 training steps and 0.4 learning rate. The two algorithms use a Matérn 3/2 kernel. The ST-RCGP uses the adaptive centering and shrinking function from Section 3.
C.5. RCGP Issues from Section 2
The data is generated by adding noise sampled from a NN o = 0.5) to a latent function f(x = 3 sin on a temporal grid e [0, 1.4] with ni = points. We substitute at locations outliers drawn from a W33 a = 0.2).
The RCGP optimisation process use for kernel hyperparameters and observation noise is the one recommended in Altamirano al. (2024). Note that we conduct optimisation the original data without outliers. We have three configurations: First, with constant prior mean m(x) equal to the data average and QN (0.9); second, with m(x) sin (2Tx) and с = Q (0.9); third, with m(x) sin(2rz) and c = 0.8. All configurations use a Squared Exponential kernel and are separately optimised (that is, they do not necessarily share hyperparameters) Note that inference is conducted on data with outliers.
C.6. ST-RCGP Posterior when Varying and |
We generate = n, = 100 data points from a GP prior with prior mean m = 0 and Matérn 3/2 kernel with f= 0.2, ox 2.0. We add noise from a (0, o²) where x² = 0.25. Outliers are generated by adding noise at temporal locations drawn from a (0 with 20.
24

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
The ST-RCGP Matérn 3/2 kernel with lengthscale & 0.198 and amplitude ox 3.01. We an adaptive centering function for all values of and Then, perform inference with varying values of , and to conduct sensitivity analysis of those hyperparameters. The results are shown in Figure 11.
We notice that increasing 3 will lead to overfitting the data. Conversely, decreasing & too much yields underconfident uncertainty estimates. The middle ground is 3 = e/v2, which yields mean and uncertainty estimates that are appropriate, thus supporting our selection of hyperparameters.
For c, Figure shows that a lower will be more robust to outliers, but may also overinflate the uncertainty estimate. On the other hand, that is too large will produce an algorithm closer to the STGP, which is not robust to outliers. Visibly, the adaptive choice for performs the best, confirming our choice of hyperparameter.
10a
Data
B 0.50
Data
5
0.01
B /=a/v2
0.05
adaptiec
100
-10
-10
0.4
0.6
0.8
0.4
0.6
0.8
Figure 11. Impact of e and B on the posterior: We keep all parameters other than 3, identical and choose me ft. Outliers are highlighted in red. When increases (or decreases), the rate at which we learn from data increases (or decreases), the confidence intervals narrow (or widen). When e decreases, increase robustness the cost of larger posterior uncertainty. These tendencies support and o0 + of where o the standard deviation of the predictive filtering posterior p(y) 111:
C.7. Well-Log Dataset
The well-log dataset, first introduced by Ruanaidh & Fitzgerald (1996), comprises 4,050 nuclear magnetic resonance measurements collected during the drilling of well. Change points in the sequence indicate transitions between sediment layers encountered by the drill. In addition to these distinct transitions, the data also includes outliers and noise caused by shorter-term geological events, such as floods, earthquakes, or volcanic activity. In Figure 12, we contrast the results obtained from the ST-RCGP and the STGP. They perform equally well in well-specified regions, but the STGP lacks robustness to outliers.
C.8. Sensitivity Analysis on IMQ Exponent
To understand how the shape of the weight function affects ST-RCGP's posteroo estimates, we conduct a sensitivity analysis. We examine how varying the IMQ exponent-ccurently a == 11/2- impacts results. We conduct this analysis because Proposition 3.3 shows robustness requires weights decay faster than 1 1/v yll i.e., a -1/4; but, overly fast decay can reduce statistical efficiency (overly robust) highlihhin a tradeoff worth exploring. The result of this analysis is shown in Figure 13.
C.9. Sensitivity Analysis on Centering and Shrinking Functions
To capture the impact of the centering and shrinking functions on the posterior results, we conduct an analysis where we compare ST-RCGP and RCGP when all parameters are kept as in RCGP, apart from the centering function. We further explore the impact of also altering the shrinking function. To do so, we generate data with outliers the same way as in Appendix 12. We then plot the posterior distributions both algorithms in Figure 14 with a confidence interval (CI) of 30.
25

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
Nuclear Magnetic Resonance Data
STGP
2
-4
ST-RCGP
-4
500
1000
1500
2000
2500
3000
Time Step
Figure 12. Well-log data. The top and bottom panels show the STGP and ST-RCGP fits the data, respectively.
C.10. Effect of Outlier Timing on ST-RCGP Posterior
We investigate the impact of changing the point in time at which outliers are introduced on the posterior of the ST-RCGP. The data is generated as in Appendix C.13. The results are shown in Figure 15.
C.11. Choices Centering Function in Special Cases
In Table 5, we highlight a few potential alternative choices other than ST-RCGP's filtering predictive f, for the centering function depending on whether the data is temporal spatio-temporal. Also, demonstrate which choices of centering functions recover the STGP and RCGP. The following points briefly explain the relevant features of each choice:
Data: This choice, which yields constant weights and thus recovers the STGP, implies that our best estimate of the center of the data is the current observation. Interestingly, if the observation is outlier, this implies we center the algorithm around the outlier. This provides intuitive understanding for why the STGP fails to be robust.
Prior Mean: Choosing for our centering function the prior mean m(x;) recovers the RCGP (assuming a constant for the ST-RCGP). However, as previously explored in Section 4, this choice can yield poor results when the prior mean aligns with the outliers.
This choice requires an additional hyperparameter (potentially many) that comes from K., The matrix K, dictates how much we want to smooth our data at time step tg. This can useful when there are unusual spatial structures in the data that K, capture. Alternatively, when there are few time steps, the filtering predictive might not be the best estimate of the center of the data, in which case good way to estimate the center that is more appropriate than simple average is to use K..
Spatial Smoothing:
Temporal Smoothing: This centering function employs the same concept as the "Spatial Smoothing". It requires a lookback period n; and weights v/ that determine how important the datum at step i is estimate the center of the data at step &.
Filtering Predictive: This is the choice we make in this paper for ST-RCGP and explained in Section 3. It involves
26

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
Impact on ST-RCGP Posterior When Varying IMQ Exponent
20
STGP (0)
1/32
10
1/4
10
20-
1/2
3/4
10-
10
0.0
0.2
0.4
0.6
0.8
1.8.0
0.2
0.4
0.6
0.8
Figure Impact on ST-RCGP when varying the IMQ exponent. The red-shaded plots indicate values of a that violate the robustness condition, ie. > -1/4. The figures are denoted by |ax), so that, for example, the top-right panel corresponds to -1/32. Outliers are highlighted red.
no additional hyperparameters, and is our model's best (and robust) prediction of yı given past observations. These reasons are why we choose this option over the other for the ST-RCGP.
Smoothed Predictive (S): Same concept "Spatial Smoothing," but the smoothing is applied on predictions instead of on the data.
Smoothed Predictive (T): Same concept as "Temporal Smoothing." but the smoothing applied on predictions instead of on the data.
C.12. Experiments Showing we Fix Vanilla RCGP from Section
Data We generate n, = 200 data points at evenly spaced inputs in x = [0. 1] from GP with Matérn 3/2 kernel with
lengthscale f 0.1 and amplitude o² = 2, and mean function m(x)) = 2e 2e-5x 5xe The sample function (f1,.... .J, fn,) drawn then centered: f + f. where f ne D ft. Noise is added and drawn from a (0, o2) with o² 0.25. We contaminate the data with 10 outliers lufl for i = 1, ..., 10 sampled from a A5, o² = 1). These outliers replace the data at specified locations. We keep both the original data and the contaminated data for further tasks.
Fit To fit the RCGP, we use the code from and follow Altamirano et al. (2024). We choose constant weight function equal to the mean of the data, and c (0.95) since there are 5% outliers. The hyperparameters are optimized on the original data (without outliers) since training RCGP on contaminated data would result in overfitting the outliers and unreliable predictions. However, RCGP predictions are made on contaminated data.
27

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
Identical Weights
Data
5
ST-RCGP
CI (RCGP)
CI (ST-RCGP)
Prior Mean
Weights Differ Centering Function
Weights Differ in Centering and Shrinking Functions
5-
0.2
0.6
1.0
Figure 14. ST-RCGP ond RCGP comparison. We contrast the two algorithms progressively change the specification of the weight function. The CI corresponds to the 30 confidence interval. The prior is the constant function used in the weight function of the RCGP and the ST-RCGP in the first plot. In the second plot, the centering function y of the ST-RCGP is the predictive mean. In the third plot, the shrinking function is the filtering predictive's covariance, as specified in Section 3.
Table Choice of Centering Function TYk For Weights wa Wi's for ni, ..., & are normalised weights selected from the ni past time step tx, (lookback period). K. E R Xna is row-wise normalized kernel matrix.
Algorithm Centering 7xx)
Description
STGP
Yk m(xx)
Data
RCGP
Prior Mean
K.yk
Spatial Smoothing
Σ, thay
Temporal Smoothing
ST-RCGP
f, = Hmaje-
Filtering Predictive
K.f.
Smoothed Predictive (S)
Smoothed Predictive (T)
To fit the ST-RCGP, we use the Adam optimiser and the robust scoring objective PGB from Section 3. The hyperparameter selection as in Section Our learning rate is 0.3, and the number of optimisation steps 70 (30 would be enough; we'r 've
28

---

Robust and Conjugate Tempio-Temppral Gaussian Processes
Observed Data
ST-RCGP
CI (ST-RCGP)
5
0.0
0.2
0.4
0.6
0.8
1.0
t
Figure 15. Impact of outlier timing on ST-RCGP Posterior.
used more to study convergence). In contrast to RCGP, the ST-RCGP is both training and predicting on contaminated data.
Coverage The coverage values computed given prediction and standard deviation For each quantile, find corresponding z-score, and determine the proportion of data points falling within + 20.
C.13. Experiments in Well-specified Settings from Section 4
select the outlier-free data the contaminated data. We use a Matérn 3/2 kernel for the GP prior. First, we perform hyperparameter optimisation on the dataset with outlier-free data. This involves 25 training steps using the Adam optimiser and learning rate of 0.3. The criterion we use our optimisation objective is the standard yo (since there are no outliers). Second, we obtain the performance metrics for each model. This done by generating new data as previously (from Appendix C.4). For each newly generated data, we compute statistical efficiency, RMSE and NLPD for each model. We take the average and standard deviation report our metrics.
We the same dataset in Appendix C.4, from which
C.14. Experiments with Financial Crashes from Section
Twitter Flash Crash We retrieve the "close" data from the DJIA index on April 23rd, 2013, and the previous day. This amounts to 810 data points. build an evenly spaced temporal grid from to with 810 points. The observations are then standardised.
The GP fit is implemented in Python's 's sklearn package and uses Matérn 3/2 kernel with amplitude d, 0.72, lengthscale = 0.0955, observation noise 0.02, and prior mean O. The RCGP fit uses Matérn 3/2 kernel with o, = 1, lengthscale & = 0.09, a and a constant prior mean equal to the average of the data. Also, RCGP has c 0.25, since offers a more robust posterior than
Q. The ST-RCGP fit uses a Matérn 5/2 kernel with amplitude d = 1., lengthscale = 0.1, -² = 0.02 and has an adaptive shrinking and centering function, as specified in Section 3. The RP fit is exactly the one from Ament et al. (2024) since is obtained using the same code.
29

---

Robust and Conjugate Spatio-Temporal Gaussian Processes
The execution time computed post-optimisation of each method, since we wish capture execution time at inference. Also, to avoid caching and establish a fair comparison, each model has second instance specifically for inference-making that basn't observed data yet but has the optimised hyperparameters.
Index Futures with Synthetic Crash The data is obtained from https: / http:/wwwww.g...... kaggle com/ c// caltech- cs15c-155---- 202002222ccc It captures an Index Futures price over time, measured at 500ms intervals. We select N 46800 data points, which amounts trading day. We built temporal grid between and 46800 and rescaled it by 0.5 (500 x 0.001). The observations are then standardised. The crash induced aims to mimic the Twitter crash incident, but with slightly more outliers. Therefore, we drop data points by (0.9995, 0.9994, 0.9992, 0.996, 0.994, 0.998, 0.998, 0.9998) of their original value (not standardised) create V-shaped outlier region and add random noise to the drop sampled from ((00 σ)) for d 0.0001. Note that the amplitude of this drop is roughly similar to that of the Twitter flash crash experiment.
The STGP has a Matérn 3/2 kernel with L = 6, oR = 6, and observation noise o = 0.14. The ST-RCGP has a Matérn 3/2 kernel with L = 6.5, Øx = 1, and observation noise a = 0.3. It also uses an adaptive centering and shrinking function as in Section 3. The BayesNewton methods all use a Student-T likelihood with df = (degrees of freedom), except for the Laplace method, which has df They are also optimized following the code from Wilkinson et al. (2023).
The execution times are computed 5 times and evaluated The one-step cost for each method is obtained by evaluating execution time for an increasing number of data points (N 5, 10, 100, 500, 1000, 2000, 2500, .... 35000, 40000, 46800), and then performing linear regression. The slope corresponds to the one-step cost. The RMSE and NLPD standard deviations are obtained by repeating inference on a newly generated crash event total of 20 times.
new instance of the model's class to avoid caching issues.
In addition Table 2, Figure 16 illustrates the fit of each method the data.
3450
33400
3820
3300
3750
Data
3760
Outlies
VI
3740
STGP
MEP
3720
MLGP
ST-RTGP
100
400
Time (300lme))
Figure Fitting the STGP, ST-RCGP and some methods from the BayesNewton package to the Index futures data with a synthetically induced crash.
C.15. Experiments with Temperature Forecasting from Section 4
The data is from the Climate Research Unit (CRU) and is available at https: httpa::://r/daaaaaaa uea. ac htt//ruudata.. hrg/. We select latitude and longitude ranges of [45, 60) and [-12, 8) respectively, which amounts to n, = 571 spatial locations per time step. The data monthly, starting in January 2022 and ending in December 2023, which is total of n, 24 time steps. This leads to N n4 n. 11, 9919 data points. We add focussed outliers on a patch with latitudes and longitudes in [51.25, 53.25), [-3, -1] respectively. The outliers are drawn from a A (120,, o²) with 10. Before fitting the data, we pre-process it by standardising.
To fit the data, we use the standard objective (since the data has been cleaned by the CRU beforehand, see Harris et al. (2020)) with Adam optimiser, 60 optimisation steps and 0.05 learning rate. Both the STGP and the ST-RCGP use a Matérn 3/2 kernel for the temporal and spatial kernels. For STGP, this yields the following hyperparameters: A Temporal amplitude
30

---

Figure 17. STGP and ST-RCGP Coverage during the month of Octoben. Robust and Conjugate Spatio-Temporal Gaussian Processes
of 0x 0.76, temporal lengthscale e 2.51, spatial amplitude xx, 0.76, spatial lengthscale L, 2.42, and variance a² = 0.15. For ST-RCGP, the hyperparameters
A Temporal amplitude of de, 1.12, temporal lengthscale & 2.13, spatial amplitude xx, 0.82, spatial lengthscale t 2.34, and variance o² 0.069. In Figure 17, we demonstrate the coverage of ST-RCGP and STGP on the temperature forecasting experiment for the month of October, which occurs before the introduction of outliers.
ST-RCGP
STGP
1.0
0.8
Coverage
0.6
0.4
0.2
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
Quantiles
Quantiles
31


Reviews:
[{"review_id": "3TeZFGWOrJ", "reviewer": "Reviewer_G2NE", "summary": "This paper expands the robust and conjugate GP framework to what it refers to as spatiotemporal GPs (sometimes referred to in other places as Markovian GPs, linear time GPs, or state space GPs). This is achieved by a generalized Bayes filtering solution, somewhat similar to other recent works on sequential generalized Bayes. Additionally, by considering the sequential nature of generalized Bayesian filtering, the authors are able to address some difficulties in \"vanilla\" robust and conjugate GPs that arise from finding appropriate weighting functions. The proposed method is tested on a number of synthetic and real-world datasets to show its effectiveness, both in spatiotemporal and 1-dimensional settings.\n\n## Update after rebuttal\n\nThe authors have clearly answered my questions, softening their claims and pointing to related work where appropriate. In addition, I find the additional notebooks provided to other reviewers in the rebuttal process nice. I thus raised my score from a 3 to a 4 post-rebuttal.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "We appreciate the thoughtful feedback and the expertise on the matter. We are glad you found the paper clearly written and easy to follow. We comment on the feedback below:\n\n**1) “provides inferences that are comparable to state-of-the-art non-Gaussian STGPs in the presence of outliers, but at a fraction of the cost\" (...) I find this claim a bit too general.**\n\nWe acknowledge that our claim is too broad and will adjust it as follows:\n“The ST-RCGP provides inferences comparable to state-of-the-art non-Gaussian STGPs with a computational cost similar to classical STGPs.”\nAdditionally, we believe further accelerating the ST-RCGP is an exciting direction for future work, so we will replace lines 422-426 (right side) with: \n“Our method builds on classical STGPs, which remain vulnerable to scaling in the spatial dimension and are not optimal if parallel computing is available. These issues were addressed by Hamelijnck et al. (2021) through variational approximations and parallel-scan algorithms (see also [1]), and could be adapted to ST-RCGPs, providing the same benefit.”\n\n**2) Regarding the EWR: it is stated that \"since the STGP is not robust to outliers, the closer EWR is to one, the less robust our posterior is to outliers, and vice-versa.\" I don't think this statement is strictly accurate**\n\nWe agree and will replace the highlighted sentence with: “Therefore, since the STGP is not robust to outliers, EWRs that are near one are not necessarily optimal, since an EWR of 1 implies a solution—the vanilla STGP—which is not robust.\"\n\n**3) I think it could avoid possible confusion to more clearly state that Proposition 3.1 assumes the modified Fisher divergence.**\n\nThank you for raising this point. We will ensure to make this clearer in the camera-ready version by replacing the current text in Proposition 3.1 with: “Then, we obtain a (generalised) posterior (…) with the score-matching loss $\\mathcal{L}$ given by X for Y” where X will be replaced by the loss as provided in lines 674-675, and Y the mathematical objects used in the loss, which are defined on line 677.\n\n**4) There is a small body of work using 3σ-rejection-esque filters for Markovian GPs that are not cited.**\t\t\t\n\n​​Thank you for bringing this to our attention. We will include the following statement after the last sentence of the paragraph on line 066:\n“We also highlight a small body of work that uses outlier-rejection Kalman filters with STGPs to improve robustness in inference tasks such as prediction and anomaly detection. While these methods offer robustness at lower computational cost compared to non-conjugate STGPs, they are generally considered less expressive and not as strongly supported by theoretical foundations.”\n\n**5) minor editorial comments**\t\t\t\t\t\t\t\t\n\nWe agree with these editorial comments and will make the necessary changes in the final version.\n\n**6) It is written \"assuming a stationary and separable kernel (...)\" but stationarity is not sufficient for representation as an LTI SDE (e.g., the squared exponential kernel).**\n\nThank you for highlighting this. We will update the final version accordingly.\n\n**7) Are there obvious technical barriers in also combining these works, i.e., having computation-aware ST-RCGPs?**\n\nWe believe there should be no fundamental technical barriers to combining these approaches because the filtering update in ST-RCGP closely mirrors that of the standard STGP. Given this structural similarity, the main arguments presented by Pförtner, Marvin, et al (2024) [2] should remain applicable.\n\n**8) How does the smoothing work?**\n\nThe generalised Bayes approach used on the filtering distribution produces a Gaussian posterior, and thus, the resulting smoothing distribution remains Gaussian. Therefore, as for classical STGPs, we can apply RTS smoothing using the GB filtering posterior. \nTo clarify ambiguities about smoothing, the above explanation will be added in the paragraph on lines 224-235 (Left side, Methodology section). \n\nFinally, we want to thank the reviewer again for the careful consideration given to our paper. We hope we have responded to any remaining questions or concerns in this rebuttal.\n\n**References**\n\n[1] Särkkä, Simo, and Ángel F. García-Fernández. \"Temporal parallelisation of Bayesian smoothers.\" IEEE Transactions on Automatic Control 66.1 (2020): 299-306.\n\n[2] Pförtner, Marvin, et al. \"Computation-aware Kalman filtering and smoothing.\" arXiv preprint arXiv:2405.08971 (2024).\n\n[3] Bock, Christian, et al. \"Online time series anomaly detection with state space Gaussian processes.\" arXiv preprint arXiv:2201.06763 (2022).\n\n[4] Waxman, Daniel, and Petar M. Djurić. \"A Gaussian Process-based Streaming Algorithm for Prediction of Time Series With Regimes and Outliers.\" 2024 27th International Conference on Information Fusion (FUSION). IEEE, 2024.", "claims_and_evidence": "The article makes a few main claims. The first is that the proposed generalized Bayes filtering solution is indeed a realization of the robust and conjugate GP framework. In particular, based on the chosen hyperparameters, it is claimed that the solution can be made robust (in the Huber sense), and that the vanilla GP can be recovered. A secondary set of claims state that several issues with vanilla RCGPs are ameliorated with some proposed changes to the weighting function. I find both of these sets of claims to be well-evidenced. \n\nThe only claim I find problematic is that the method \"provides inferences that are comparable to state-of-the-art non-Gaussian STGPs in the presence of outliers, but at a fraction of the cost\": Hamelijnck et al. show that their \"parallel\" formulation provide GPU implementations of MVGPs that are potentially orders of magnitude faster than the sequential formulation on CPU (c.f. Figure 4 in the Supplementary Material of [1]). To my knowledge, this formulation would not be possible with RC STGPs; the reviews indicate that all experiments are on CPU, in which case I find this claim a bit too general.\n\nI also have some minor technical comments (see the points in \"Other Comments or Suggestions\" below)", "methods_and_evaluation_criteria": "I believe the proposed methods and evaluation criteria make sense for the problem at hand. My only technical comment is regarding the EWR: it is stated that \"since the STGP is not robust to outliers, the closer EWR is to one, the less robust our posterior is to outliers, and vice-versa.\" I don't think this statement is strictly accurate. For example, if one downweights in the \"wrong\" locations, robustness is actually decreased with the EWR. Perhaps a more correct statement is just that EWRs that are near one are not necessarily optimal, since an EWR of 1 necessarily implies a solution which is not robust. I understand the sentiment of this statement, however, and find it an intuitive and appropriate measure for experiments.", "theoretical_claims": "I have checked Proposition 3.1 carefully. My only concern is that the proof depends on a particular choice of $\\mathcal{L}$, which is only briefly mentioned in the background section. I think it could avoid possible confusion to more clearly state that Proposition 3.1 assumes the modified Fisher divergence. I went through Proposition 3.2 more briefly, but did not spot any issues. It's not entirely clear to me how smoothing works (asked below).", "experimental_designs_or_analyses": "I believe the experimental design and analysis are sound, besides the potential omission of the \"parallel\" variational STGPs as mentioned above.", "relation_to_broader_scientific_literature": "Robust and scalable solutions to GPs are increasingly important. The submission does a good job of highlighting its relation to two main subsets of the GP literature: robust solutions -- particularly through generalized Bayes (such as RCGP), and scalable spatiotemporal solutions through state-space representations. Besides approximate or non-conjugate GPs, there is a small body of work using $3\\sigma$-rejection-esque filters for Markovian GPs that are not cited [2, 3]; these can be seen as creating a more naive (and degenerate) weighting function. Importantly, both of these works are \"adaptive\" like the ST-RCGP, in the sense that the current posterior predictive moments are used sequentially to detect anomalies/outliers.", "essential_references_not_discussed": "Besides approximate or non-conjugate GPs, there is a small body of work using $3\\sigma$-rejection-esque filters for Markovian GPs that are not cited [2, 3]; these can be seen as creating a more naive (and degenerate) weighting function. Importantly, both of these works are \"adaptive\" like the ST-RCGP, in the sense that the current posterior predictive moments are used sequentially to detect anomalies/outliers.", "other_strengths_and_weaknesses": "### Strengths\n- The paper is generally well-written and easy to follow.\n### Weaknesses\n- The paper might be seen as a straightforward combination of a few different ideas: STGPs and generalized Bayes filters. However, the benefits of adaptive weight centering and the thoroughness of experiments outweigh this weakness, in my opinion.", "other_comments_or_suggestions": "I have a few minor editorial comments:\n1. (Line 082, Right Column) It is stated that \"If standard GPs were used here, the cost would be $\\mathcal{O}(n_t^3 n_s^3)$, which would be impractical,\" which I find unnecessarily strong. I suggest that this \"may become impractical,\" instead.\n2. (Line 224, Right Column) I think this should read \"accomplishing the purpose of [principle] 3\".\n3. (Section 4) I understand that space constraints are an issue, but I think it would improve the exposition to motivate the EWR in a sentence or two, rather than just pointing the reader to a technical definition in the appendix.\n\nAnd a small technical comment:\n1. (Line 087, Right Column) It is written \"assuming a stationary and separable kernel [...]\" but stationarity is not sufficient for representation as an LTI SDE (e.g., the squared exponential kernel). Notably, however, any stationary kernel can be approximated arbitrarily well (c.f. [4]).", "questions_for_authors": "1. At several points, the computational-aware RCGP is mentioned. [5] also takes a similar step in proposing computation-aware STGPs. Are there obvious technical barriers in also combining these works, i.e., having computation-aware ST-RCGPs?\n2. How does the smoothing work? I see that there are smoothing solutions in the code, but it is not explained in the text. Does one obtain the \"correct\" smoothing solution by simply using the stored GB filtering statistics and applying RTS smoothing naively?\n\n### References\n[1] Hamelijnck, Oliver, et al. \"Spatio-temporal variational Gaussian processes.\" Advances in Neural Information Processing Systems 34 (2021): 23621-23633. \n\n[2] Bock, Christian, et al. \"Online time series anomaly detection with state space Gaussian processes.\" arXiv preprint arXiv:2201.06763 (2022). \n\n[3] Waxman, Daniel, and Petar M. Djurić. \"A Gaussian Process-based Streaming Algorithm for Prediction of Time Series With Regimes and Outliers.\" 2024 27th International Conference on Information Fusion (FUSION). IEEE, 2024.\n\n[4] Loper, Jackson, et al. \"A general linear-time inference method for Gaussian Processes on one dimension.\" Journal of Machine Learning Research 22.234 (2021): 1-36.\n\n[5] Pförtner, Marvin, et al. \"Computation-aware Kalman filtering and smoothing.\" _arXiv preprint arXiv:2405.08971_ (2024).", "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": "", "extra_scores": {}}, {"review_id": "Mq7e397IRv", "reviewer": "Reviewer_cvUj", "summary": "This paper introduces a methodology for spatio-temporal Gaussian Processes based on a state-space model and generalized Bayesian inference. Building on the robust and conjugate Gaussian processes (RCGPs) framework, it addresses and overcomes its key limitations, enhancing both robustness and computational efficiency.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "We thank the reviewer’s feedback and are glad our theoretical results and experimentation on our proposed method were appreciated. It has been mentioned that:\n\n**“The baseline models used for comparison vary across experiments due to data characteristics. It would be beneficial for the authors to include an additional state-of-the-art model beyond RCGP to better assess the overall effectiveness of the proposed method.”**\n\nWe agree that it would be beneficial to illustrate the fit of additional methods apart from RCGP. For this reason, we will include in the paper the following plot: https://anonymous.4open.science/r/ST-RCGP-21DD/experiments/financial-applications/high-frequeny-data/HFT-data-induced-crash-fit.pdf, which was produced by code available here:\nhttps://anonymous.4open.science/r/ST-RCGP-21DD/experiments/financial-applications/high-frequeny-data/HFT-index-futures-speed-comp.ipynb. This plot compares ST-RCGP against some of the state-of-the-art methods offered by the “BayesNewton” package and is an extension of Table 2. \n\n\nAgain, we are grateful for the consideration given to our paper and hope our answer satisfies the concern raised.", "claims_and_evidence": "The claims in this submission are well-supported by rigorous mathematical proofs and empirical experiments.", "methods_and_evaluation_criteria": "The proposed methods are evaluated against both standard Gaussian Processes and the state-of-the-art RCGP framework, using multiple datasets and diverse evaluation metrics.", "theoretical_claims": "Overall looks good to me.", "experimental_designs_or_analyses": "The experimental design looks good to me.", "relation_to_broader_scientific_literature": "This paper extends the setting of previous work to spatio-temporal modeling by incorporating a state-space formulation, as discussed in prior literature. This formulation addresses three key limitations of the earlier approach.", "essential_references_not_discussed": "None.", "other_strengths_and_weaknesses": "Strengths:\n- Provides strong theoretical results.\n- Evaluates the proposed ST-RCGP method across multiple experiments using various metrics.\n\nWeaknesses:\n- The baseline models used for comparison vary across experiments due to data characteristics. It would be beneficial for the authors to include an additional state-of-the-art model beyond RCGP to better assess the overall effectiveness of the proposed method.", "other_comments_or_suggestions": "None.", "questions_for_authors": "See above comments.", "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": "", "extra_scores": {}}, {"review_id": "3UBiz4vJpS", "reviewer": "Reviewer_vhF8", "summary": "This paper introduces Spatio-Temporal Robust and Conjugate Gaussian Processes (ST-RCGPs), which are an extension of robust and conjugate Gaussian processes (RCGPs) [1]. ST-RCGPs leverage the state-space formulation of spatio-temporal GPs (STGPs) to achieve computational efficiency while maintaining the robustness properties of RCGPs. The paper also addresses three key limitations of RCGPs: sensitivity to prior mean misspecification, poor uncertainty quantification, and reliance on manual hyperparameter tuning. Robustness to outliers is established by showing that the posterior influence function is bounded. Experiments on synthetic, financial, and weather data compare ST-RCGP to STGP, RCGP, and other baselines in settings with outliers.\n\n[1] Altamirano, M., Briol, F., Knoblauch, J. \"Robust and Conjugate Gaussian Process Regression\". ICML 2024.\n\n---\n## Update after rebuttal\nThe rebuttal has addressed my concerns and answered my questions. I have raised my score from 3 to 4.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "We thank the reviewer for their thoughtful review and positive remarks on our empirical evaluation and the ST-RCGP’s novel combination of robustness and efficiency. We comment on the feedback below:\n\n**Since the data is synthetic, I have some concerns about how meaningful the results really are.**\n\nRobustness studies often inject outliers due to the lack of datasets with labelled outliers and the ability to control testing across varied outlier types (as was done in our work). However, we will add an experiment in the paper with the well-log dataset from [8, 11], which consists of 4,050 nuclear magnetic resonance measurements recorded while drilling a well. Outliers were identified scientifically and correspond to short-term events in geological history. The results can be found here: https://anonymous.4open.science/r/ST-RCGP-21DD/experiments/well-log-data/NMR-data-fit.pdf with the code in the same folder.\n\n\n**The statement of Proposition 3.2 does not seem to be very strong.**\t\t\t\n\nWhilst we empathise with this comment, the study of the influence function (IF) dates from the seminal work of Hampel in 1968 [1] and has been the gold standard in robust statistics ever since [2,3,4]. The typical approach to show robustness is to bound the influence function. While one could work out explicitly the bounds, there is no clear gain in doing so as the values of the KL divergence are not very interpretable. ​​To prove the robustness of an algorithm, bounding the IF is then sufficient and remains one of the strongest (and in many cases the only) available criteria within the robust Bayesian literature [see e.g. 5,6,7].\n\n**Can you comment on any alternatives to the IMQ kernel for assigning the weights?**\n\nWe choose the IMQ for several reasons. To elaborate on those, in the final version, we change lines 237-243 with:\n“The statements of Propositions 3.1 and 3.2 impose some constraints on the choice of weight function. In particular, it should be strictly positive and differentiable over its domain to ensure quantities in Proposition 3.1 are well-defined, and it must have a bounded supremum to satisfy Proposition 3.2. Moreover, to avoid attributing weight to arbitrarily large $y\\in \\mathbb{R}$, we require $\\lim_{|y| \\rightarrow \\infty}w(\\mathbf{x}, y) = 0$ (decaying property). These requirements make the IMQ an appropriate choice of weight function, in addition to the fact that it has been well-studied and recommended in prior literature [6, 9, 10,11]. We further justify our choice on the basis that the IMQ hyper-parameters γ,$\\beta$ and c are related to concepts of robustness—relations we exploit to specify the IMQ. In particular, to select γ, β and c, we follow four guiding principles; we want to: …”\n\nWe acknowledge the need for a sensitivity analysis to understand how the shape of the weight function affects ST-RCGP's posterior estimates. In the camera-ready version, we will examine how varying the IMQ exponent—currently $\\alpha:=-1/2$—impacts results. We conduct this analysis because Proposition 3.2 shows robustness requires weights to decay faster than $1/\\sqrt{|y|},$, i.e., $\\alpha<-1/4$; but, overly fast decay can reduce statistical efficiency (overly robust)—highlighting a tradeoff worth exploring. The result of this analysis can be found at https://anonymous.4open.science/r/ST-RCGP-21DD/tests/IMQ-sensitivity/IMQ-exponent-testing.pdf, with the code in the same folder.\n\n\nThank you again for the thorough review. We hope this rebuttal addresses any remaining concerns and contributes to a stronger evaluation of our paper.\n\n[1] Hampel, Frank Rudolf. Contributions to the theory of robust estimation. 1968.\n\n[2] Huber, Peter J., and Elvezio M. Ronchetti. Robust statistics. 2011.\n\n[3] Maronna, Ricardo A., et al. Robust statistics: theory and methods (with R). 2019.\n\n[4] Hampel, Frank et al. Robust statistics: The approach based on influence functions. 1987.\n\n[5] Ghosh, A. and Basu, A. Robust Bayes estimation using the density power divergence. Annals of the Institute of Statistical Mathematics, 68(2), 2016.\n\n[6] Matsubara, T., Knoblauch, J., Briol, F.-X., and Oates, C. J. Robust generalised Bayesian inference for intractable likelihoods. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2022.\n\n[7] Duran-Martin, Gerardo, et al. \"Outlier-robust Kalman filtering through generalised Bayes.\" 2024.\n\n[8] Altamirano, Matias, François-Xavier Briol, and Jeremias Knoblauch. \"Robust and scalable Bayesian online changepoint detection.\" International Conference on Machine Learning. PMLR, 2023.\n\n[9] Chen, Wilson Ye, et al. \"Stein point markov chain monte carlo.\" International Conference on Machine Learning. PMLR, 2019.\n\n[10] Riabiz, Marina, et al. \"Optimal thinning of MCMC output.\" Journal of the Royal Statistical Society Series B: Statistical Methodology 84.4 (2022).\n\n[11] Ruanaidh, J. J. O. and Fitzgerald, W. J. Numerical Bayesian methods applied to signal processing. Springer Science & Business Media, 1996.", "claims_and_evidence": "The paper claims that\n- ST-RCGP is an \"outlier-robust spatio-temporal GP with a computational cost comparable to classical spatio-temporal GPs\"\n  - The robustness to outliers is demonstrated empirically (see e.g. Figures 4, 5, and 6, and Table 3) and theoretically (see Proposition 3.2). The improved computational efficiency is also demonstrated empirically (see e.g. Figure 5)\n- ST-RCGP can \"overcome the three main drawbacks of RCGPs: their unreliable performance when the prior mean is chosen poorly, their lack of reliable uncertainty quantification, and the need to carefully select a hyperparameter by hand.\n  - This claim is addressed in Section 3 by adapting the weight function throughout filtering steps, whereas RCGPs require it to be fixed.", "methods_and_evaluation_criteria": "In terms of methods, the main idea is to combine the existing RCGP framework with state-space GPs for improved computational efficiency in spatio-temporal domains, which seems sensible to me. The fact that this allows for an adaptive weight function with additional benefits is great and also makes sense.\n\nIn terms of evaluation, the experiments seem to be somewhat contrived but reasonable to test the model. For example, the outliers introduced in the temperature forecasting experiment seem to be kind of unrealistic. The quantitative evaluation uses RMSE, NLPD, and compute time as metrics, which can be considered as standard in the literature.", "theoretical_claims": "- Proposition 3.1 derives the state-space formulation of RCGPs. The resulting equations seem reasonable and the proof is provided in Appendix B.1, although I did not check the latter carefully.\n- Proposition 3.2 shows that the posterior influence function is bounded (in its first argument) to demonstrate theoretical robustness. The assumptions seem to be mild, but the result is also not very strong (just because a quantity is not infinite it could still be really, or even arbitrarily, large...). The proof is provided in Appendix B.2., but I did not verify its correctness.", "experimental_designs_or_analyses": "The paper presents four different experiments:\n1. \"Fixing vanilla RCGP\"\n  - This experiment compares ST-RCGP to RCGP on 1D synthetic data to demonstrate that ST-RCGP addresses the three issues outlined  in Section 2 (sensitivity to prior mean, poor uncertainty quantification, selection of shrinking function necessary). The visualization (see Figure 4) looks convincing, although it is only a single example which could be selected on purpose.\n2. \"ST-RCGP in Well-Specified Settings\"\n  - This experiment compares ST-RCGPs to regular STGPs and RCGPs on synthetic data, demonstrating that ST-RCGP performs better than both other methods in the presence of outliers and matches performance otherwise. Since the data is synthetic, I have some concerns about how meaningful the results really are.\n3. \"Robustness During Financial Crashes\"\n  - This experiment is taken from [1] and [2]. It considers historical financial data and demonstrates that ST-RCGP is robust to a sudden crash which produces outliers in the time series. The experiment also demonstrates the superior computational efficiency of ST-RCGP compared to RCGP, which is obtained due to the state-space formulation. In a later part of this experiment, a larger version of the data with a synthetically induced crash is considered.\n4. \"Forecasting Temperature Across the UK\"\n  - This experiment focuses on temperature predictions and demonstrates that ST-RCGP performs well in the presence of outliers, whereas the forecasting of the regular STGP is strongly influenced by the outliers, which results in incorrect predictions. Again, the outliers are synthetically introduced.\n\nWhile the presented results look quite strongly in favor of ST-RCGP, I am somewhat concerned about the amount of synthetic data used. Even the experiments involving real-world data have synthetically introduced outliers, which might not represent realistic scenarios.\n\n[1] Altamirano, M., Briol, F., Knoblauch, J. \"Robust and Conjugate Gaussian Process Regression\". ICML 2024.\n\n[2] Ament, S., Santorella, E., Eriksson, D., Letham, B., Balandat, M., Bakshy, E. \"Robust Gaussian Processes via Relevance Pursuit\". NeurIPS 2024.", "relation_to_broader_scientific_literature": "The paper is primarily combines [1] with a state-space formulation of GPs. The resulting update equations are somewhat related to the Kalman filter. An important component of ST-RCGP (and also the original RCGP) is the inverse multi-quadratic kernel (IMQ) which is used to set the weights for each data point, such that outliers can (ideally) be assigned negligible contribution. In terms of the general goal of designing a robust GP, [2] is a recently published method with the same goal.\n\n\n[1] Altamirano, M., Briol, F., Knoblauch, J. \"Robust and Conjugate Gaussian Process Regression\". ICML 2024.\n\n[2] Ament, S., Santorella, E., Eriksson, D., Letham, B., Balandat, M., Bakshy, E. \"Robust Gaussian Processes via Relevance Pursuit\". NeurIPS 2024.", "essential_references_not_discussed": "I am not aware of any essential references which were not discussed.", "other_strengths_and_weaknesses": "Strengths:\n- Novel combination of robustness (RCGP) and efficiency (state-space formulation)\n- Provides theoretical argument for robustness (Proposition 3.2)\n- Thorough empirical evaluation on four different experiments\n\nWeaknesses:\n- The statement of Proposition 3.2 does not seem to be very strong\n- Experiments make heavy use of synthetic data", "other_comments_or_suggestions": "N/A", "questions_for_authors": "1. Do you have any additional experiments / empirical evidence which does not rely on synthetic data? If I understand correctly, only Figure 5 considers outliers from a real-world scenario.\n2. Can you comment on any alternatives to the IMQ kernel for assigning the weights? This seems like a really important component for the method, yet the choice of IMQ seems somewhat arbitrary right to me.", "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": "", "extra_scores": {}}, {"review_id": "1wh9OiFRM4", "reviewer": "Reviewer_Gtm6", "summary": "The authors combine ideas from recent work on robust Gaussian processes with filtering ideas used in (spatio)temporal Gaussian process regression. They use temporal structure of the problem to set parameteres in the robust Gaussian process framework proposed by Altimirano et al 2024 (RCGP) in a sensible and automatic manner. They show that a computational speedup is possible using the filtering approach, allowing RCGPs to be scaled to spatio-temporal Gaussian processes with many time steps.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "We thank the reviewer for their careful consideration of our paper and helpful feedback. We address below the valuable suggestions made to further improve our work: \n\n**Scaling is shown to be linear in the number of time steps (but not in the number of spatial locations).**\n\nWe agree with this concern. In the camera-ready version, we will replace the claim with “An alternative approach which has linear-in-time cost is to use a state-space representation.”\n\n**It would be useful (...) if the authors ran an experiment in which the centering function was adaptive, but all other parameters were kept as in RCGP to show that this does in fact reduce sensitivity.**\n\nWe agree that isolating the effect of the adaptive centering function would strengthen our claim regarding reduced sensitivity. We will add an experiment comparing ST-RCGP and RCGP where all parameters are kept as in RCGP apart from the centering function. We further develop it by also altering the shrinking function. This analysis can be found here: https://anonymous.4open.science/r/ST-RCGP-21DD/tests/RCGP-vs-ST-RCGP-centering-function-sensitivity-analysis.ipynb. We welcome any further feedback.\n\n**Consider including coverage plots for at least one other dataset, for example the temperature dataset considered.**\n\nWe agree that having more datasets for which we demonstrate the (frequentist) coverage of ST-RCGP would help further strengthen our claim. We will include in the paper a coverage analysis of the ST-RCGP for the temperature dataset. This analysis can be found at https://anonymous.4open.science/r/ST-RCGP-21DD/experiments/weather-forecasting/model-fitting.ipynb (at the end of the notebook). \n\n**Why is decontaminated data used for selection of hyperparameters? Wouldn't it be more realistic to use the contaminated data?**\n\nFitting the RCGP on contaminated data leads to poor estimates due to overfitting outliers during hyperparameter optimisation. Using decontaminated data improves RCGP results. This choice thus strengthens our claim: even with this advantage, RCGP underperforms compared to ST-RCGP.\n\n**I'm not convinced this is a meaningful timing comparison.**\n\nWe don’t include optimisation time because it can be easily skewed by tempering with the optimisation procedure, especially for STGPs, which rely on gradient-based optimisation without clear stopping criteria, number of optimisation steps, learning rate, and even optimisers. Either way, the two optimisation objectives (for STGP and ST-RCGP) are extremely similar, and their computational cost are nearly identical.\n\nFocusing on inference time is meaningful in online settings where one-step inference matters. For example, real-time applications like stock price estimation may preclude methods that iterate at each time step (e.g. variational STGP) due to time constraints, favoring closed-form solutions such as the Kalman filter, the robust filter from Durán-Martín et al. (2024), or the ST-RCGP.\n\n**Appendix C.10 are the parameter presented the parameter selected after optimization, or the initialized values of parameters?**\n\nThe parameters presented are the ones post-optimization.\n\n**There should be a more clear description of the contribution relative to the prior work Duran-Martin et al, 2024. In particular, what is new?**\n\nThe key differences are:\n\n1) We deal with STGP problems such as hyperparameter optimisation and smoothing.\n\n2) We carefully specify the weight function (see answer to vhF8 on IMQ) and its parameters to improve robustness (see Weight Function paragraph on line 237), whereas Duran-Martin et al. (2024) provide few justifications for the centering and shrinking functions and learning rate, which are parameters intrinsically connected to robustness and crucial to improving performance.\n\n3) Duran-Martin et al (2024) use a weighted log-likelihood, while we use the weighted score-matching loss from Altamirano et al. (2024), leading to a different update step.\n\nWe will include in the paper, at line 226, “However, these methods do not deal with STGP problems such as hyperparameter optimisation and smoothing and do not investigate downweighting optimally.”\n\n**I don't understand the caption of figure 2.” “Comment on appendix A about abuse of notation.**\n\nWe thank the reviewer for bringing this to our attention. We will remove this sentence in the caption in the camera-ready version, and make clear that we are abusing notation. \n\n**How does the point in time at which outliers occur affect the sensitivity of your method?.**\n\nAn empirical analysis we conducted, available at https://anonymous.4open.science/r/ST-RCGP-21DD/tests/outliers-early-vs-late.ipynb, demonstrates that the method is not affected by the impact of outlier position, which we expect to be because of smoothing.\n\nFinally, we want to thank the reviewer again for the careful review and consideration given to our paper. We hope this rebuttal addresses any remaining questions and concerns.", "claims_and_evidence": "*Claim: The computational cost of the proposed method scale linearly in $N$ (made in paragraph state-space formulation)*\n\nThis claim is partially supported. Scaling is shown to be linear in the number of time steps (but not in the number of spatial locations). This claim is later stated correctly, and shown in proposition 3.1. The speed-up is also demonstrated in experiments. The others should avoid making the claim that the method scales linearly in the total number of observations, but the broader claim that the method improves scalability as compared to robust alternatives is well-supported.\n\n*Claim: Selection of the prior mean is important in RCGP. The proposed method improves selection of the prior mean using temporal structure in the problem.* \n\nBoth parts of the claim, that RCGP is sensitive to selection of the prior mean, and that the proposed method, ST-RCGP, reduces this sensitivity, seem well-supported by the figure and argument in the text. Figure 3 supports the first part of this claim. I found the second part of the claim slightly less well-supported. It would be useful, though I do not think essential, if the authors ran an experiment in which the centering function was adaptive, but all other parameters were kept as in RCGP to show that this does in fact reduce sensitivity. I didn't see such an experiment in the supplement. If it is already there I would appreciate a pointer to it, and encourage the authors to reference it when discussing issue #1 in section 2. \n\n*Claim: ST-RCGP improves upon (frequentist) coverage issues with RCGP*\nThis claim is supported by the coverage plot in figure 4. I think this claim has the least support. Consider including coverage plots for at least one other dataset, for example the temperature dataset considered. If this is not possible, could you explain why and consider including an additional simulation with spatial dimensions illustrating coverage properties of both methods?", "methods_and_evaluation_criteria": "The evaluation criteria and baseline make sense for the problem considered. Specific questions below:\n\n*Figure 1 (description in appendix C.4), \"For both the STGP and ST-RCGP, we fit the data with the optimisation objective φ and use de-contaminated data (original data without outliers) for the objective.\"* Why is decontaminated data used for selection of hyperparameters? Wouldn't it be more realistic to use the contaminated data?\n\n*Timing in figure 5, details in appendix C.10 \"The execution time is computed post-optimisation of each method, since we wish to capture execution time at inference.Also, to avoid caching and establish a fair comparison, each model has a second instance specifically for inference-making that hasn’t observed data yet but has the optimised hyperparameters..\"* \n\nI'm not convinced this is a meaningful timing comparison. It seems to me as though either optimization time should be included, as we are interested in the total cost of a procedure or caching should be used, because we only care about prediction at some new points after an algorithm has been trained. I would suggest including optimization time, or presenting both. Or could you clarify why the current approach might be meaningful for someone interested in using these methods in practice.\n\n*Appendix C.10* are the parameter presented the parameter selected after optimization, or the initialized values of parameters?", "theoretical_claims": "I looked over the proofs presented in appendix B. Although I did not go through them in sufficient detail to be confident of correctness, I did not see glaring issues, and the results seem very reasonalbe.", "experimental_designs_or_analyses": "I looked through the experimental design. Specific concerns were raised in a previous box on evaluation methods.", "relation_to_broader_scientific_literature": "Generally, discussion of prior work was good. There should be a more clear description of the contribution relative to the prior work Duran-Martin et al, 2024. In particular, what is new? It seems the methodology in your paper is based on combining the Kalman filtering approach used in spatial GPs with the robust GP method proposed in Altamirano et al 2024. Initially I thought this was the main contribution. How does this compare to what was already done in Duran-Martin et al? Should I really understand the main contribution of the paper to be the selection of better weighting and centering functions? I think either is a reasonable contribution for a conference paper, given the high quality of presentation. But I'd like the scope of what is new, and how the paper builds on current work to be a bit clearer.", "essential_references_not_discussed": "I am not aware of essential references that were missed.", "other_strengths_and_weaknesses": "Presentation in the paper is generally very good. Problems are well-motivated.", "other_comments_or_suggestions": "I don't understand the caption of figure 2. What is mean't be \"We emphasize...\" and where can I see this in the figure?", "questions_for_authors": "How does the point in time at which outliers occur affect the sensitivity of your method? In particular, it seems like at the first time step, your centering function doesn't improve on prior work since it relies on the filtering estimate. Does this mean the proposed method struggles more with outliers at early time steps than later time steps.", "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": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This work builds upon the methodology in spatio-temporal Gaussian process (GP) models where the problem is cast into a state-estimation problem. The authors extend the so-called robust and conjugate GP framework to obtain an outlier-robust spatio-temporal GP model with a comparable computation cost to vanilla spatio-temporal GPs. This submission was reviewed by 4 reviewers in the field, and the paper ended up with an average overall recommendation of 4.0 (min: 4, max: 4). The average score increased from 3.5 -> 4.0 during the discussion. All reviewers argue for accepting the paper.

This paper has clear merits. All the reviewers appreciated the presentation and proposed methodology. On the other hand, the paper might be seen as a straightforward combination of a few different ideas. However, the benefits of adaptive weight centering and the thoroughness of experiments outweigh this weakness (quoting reviewer G2NE).
Decision: Accept (poster)