Efficient Spatio-Temporal Gaussian Regression via Kalman Filtering

TL;DR

This paper introduces a scalable method for non-parametric spatio-temporal Gaussian process regression by combining Kalman filtering with assumptions of space/time separability, enabling efficient estimation from sparse, noisy data.

Contribution

It proposes a novel coupling of GP regression and Kalman filtering under specific assumptions, creating a finite-dimensional state-space model for dynamical GPs.

Findings

Kalman filter state is a sufficient statistic for minimum variance estimation.

The method scales efficiently to large spatio-temporal datasets.

Performance surpasses standard GP regression in synthetic and real data tests.

Abstract

In this work we study the non-parametric reconstruction of spatio-temporal dynamical Gaussian processes (GPs) via GP regression from sparse and noisy data. GPs have been mainly applied to spatial regression where they represent one of the most powerful estimation approaches also thanks to their universal representing properties. Their extension to dynamical processes has been instead elusive so far since classical implementations lead to unscalable algorithms. We then propose a novel procedure to address this problem by coupling GP regression and Kalman filtering. In particular, assuming space/time separability of the covariance (kernel) of the process and rational time spectrum, we build a finite-dimensional discrete-time state-space process representation amenable of Kalman filtering. With sampling over a finite set of fixed spatial locations, our major finding is that the Kalman filter state at instant $t_k$ represents a sufficient statistic to compute the minimum variance estimate of the process at any $t \geq t_k$ over the entire spatial domain. This result can be interpreted as a novel Kalman representer theorem for dynamical GPs. We then extend the study to situations where the set of spatial input locations can vary over time. The proposed algorithms are finally tested on both synthetic and real field data, also providing comparisons with standard GP and truncated GP regression techniques.

Tables

Table 1. Table 1: Comparison of the estimation Fit defined in ( 26 ) and the CPU times for the Gaussian time and space kernels ( 25 ).

	Fit [$\%$]	Memory [MB]	CPU time [sec.]
Kalman-based Alg.1 $r=6$	99.4	3	0.025
Truncated GP $q=5$ ($\sigma$)	95.1	2	0.008
Truncated GP $q=10$ ($2\sigma$)	98.5	8	0.043
Truncated GP $q=20$ ($3\sigma$)	99.3	32	0.151
Classic GP $q=\infty$	100	200	$\approx$19

Table 2. Table 2: Comparison of the estimation Fit defined in ( 26 ) and the CPU times for the Laplace time kernel ( 27 ).

	Fit [$\%$]	Memory [MB]	CPU time [sec.]
Kalman-based Alg.1 $r=1$	100	0.1	0.002
Truncated GP $q=15$	95.4	18	0.663
Truncated GP $q=30$	97.9	72	4.807
Truncated GP $q=40$	98.9	128	10.74
Classical GP $q=\infty$	100	200	$\approx$19

Table 3. Table 3: Memory and computational requirements per iteration for Algorithm  1 and both the classical and truncated GP methods applied to the Colorado Weather data-set.

	Fit [%]	Memory [MB]	CPU time [sec.]
Kalman-based Alg.1	100	4	0.02
Truncated GP (1 y. data)	99	150	15
Truncated GP (2 ys. data)	99.5	600	120
Truncated GP (3 ys. data)	99.9	1300	410
Classical GP (all data)	100	1.5 ${10}^6$	NA

Equations

$$y_i=f(a_i)+v_i,$$

$$\hat{f}(a):=\mathbb{E}\left[f(a)|\{a_i,y_i{\}}_{i=1}^N\right]=\sum _{i=1}^N{c}_{i}K(a_i,a),$$

$$\left[\begin{array}{c}{c}_1\\ ⋮\\ c_N\end{array}\right]=(\bar{K}+{\sigma }^{2}I{)}^{-1}\left[\begin{array}{c}{y}_1\\ ⋮\\ y_N\end{array}\right],\bar{K}\in {\mathbb{R}}^{N\times N},$$

$$\begin{array}{rl}V& (a)=\text{Var}\left[f(a)|\{a_i,y_i{\}}_{i=1}^N\right]=K(a,a)-\\ & \left[\begin{array}{ccc}K(a_1,a)& \cdots & K(a_N,a)\end{array}\right](\bar{K}+{\sigma }^2\mathbb{I}{)}^{-1}\left[\begin{array}{c}K(a_1,a)\\ ⋮\\ K(a_N,a)\end{array}\right].\end{array}$$

$$\begin{array}{rl}& s(k+1)=A(k)s(k)+w(k),\\ & y(k)=C(k)s(k)+v(k),\end{array}$$

$$\hat{s}(k+1|k)=A(k)\hat{s}(k|k)$$

$$\Sigma (k+1|k)=A(k)\Sigma (k|k)A(k{)}^T+Q$$

$$L(k+1)=\Sigma (k+1|k)C(k+1{)}^T$$

$${\left(C(k+1)\Sigma (k+1|k)C(k+1{)}^T+R\right)}^{-1}$$

$$\hat{s}(k+1|k+1)=\hat{s}(k+1|k)+$$

$$L(k+1)\left(y(k+1)-C(k+1)\hat{s}(k+1|k)\right)$$

$$\Sigma (k+1|k+1)=\left(I-L(k+1)C(k+1)\right)\Sigma (k+1|k)$$

$$\hat{s}(k)=\mathbb{E}\left[s(k)|y(0),\dots ,y(k)\right],$$

$$\mathbb{E}\left[s(k)|y(0),\dots ,y(k)\right]=\mathbb{E}\left[s(k)|s(k-1),y(k)\right],$$

$$S(\omega ):=\mathcal{F}[h(\tau )](\omega ).$$

$$W(\mathbf{i}\omega )=\frac{b_{r-1}(\mathbf{i}\omega {)}^{r-1}+b_{r-2}(\mathbf{i}\omega {)}^{r-2}+\cdots +b_0}{(\mathbf{i}\omega {)}^r+a_{r-1}(\mathbf{i}\omega {)}^{r-1}+\cdots +a_0},$$

$$\{\begin{array}{ll}& \dot{s}(t)=Fs(t)+Gw(t)\\ & z(t)=Hs(t)\end{array}$$

$$F$$

$$H$$

$$\mathcal{I}:=\left\{x_1,\dots ,x_{\mathrm{M}}|x_i\in \mathcal{X}\right\}.$$

$$y_i(t_k)=f(x_i,t_k)+v_i(t_k),v_i(t_k)\sim \mathcal{N}(0,{\sigma }^2),$$

$$K(x,x^{\prime },t,t^{\prime })=K_{\mathrm{s}}(x,x^{\prime })h(\tau ),\tau :=t-t^{\prime }.$$

$$\mathbf{f}(t):={\left[f(x_1,t),\dots ,f(x_M,t)\right]}^T.$$

$$\{\begin{array}{ll}& S_i:\{\begin{array}{ll}{\dot{s}}_i(t)& =F{s}_i(t)+G{w}_i(t)\\ z_i(t)& =H{s}_i(t)\end{array}i\in \{1,\dots ,M\},\text{vspace}0.2cm\\ & \mathbf{f}(t)={\bar{K}}_{\mathrm{s}}^{1/2}\mathbf{z}(t)\end{array}$$

$$\left\{\begin{array}[]{rcl}\dot{\mathbf{s}}(t)&=&\left(I\otimes F\right)\mathbf{s}(t)+\left(I\otimes G\right)\mathbf{w}(t)\\ \mathbf{f}(t)&=&\bar{K}_{\rm s}^{1/2}\left(I\otimes H\right)\mathbf{s}(t)\,.\end{array}\right.$$

$$\left\{\begin{array}[]{rcl}\mathbf{s}(k+1)&=&A(k)\mathbf{s}(k)+\mathbf{w}(k)\\ \mathbf{y}(k)&=&C(k)\mathbf{s}(k)+\mathbf{v}(k)\,,\end{array}\right.$$

$$\bar{Q}(k)={\int }_0^{T_k}\left(e^{F\tau }\right)G{G}^T{\left(e^{F\tau }\right)}^{T}d\tau ;$$

$$\hat{\mathbf{f}}(t):=\mathbb{E}\left[\mathbf{f}(t)|\{x_i,y_i(t_{\ell })\},x_i\in \mathcal{M}(\ell ),\ell =0,\dots ,k,t\ge t_k\right].$$

$$\hat{\mathbf{f}}(t)$$

$${\Sigma }^{\mathbf{f}}(t)$$

$$h(\tau )=\lambda e^{-|\tau |/{\sigma }_t}$$

$$S_r(\omega )=\sqrt{\frac{2\lambda }{{\sigma }_t}}\frac{1}{(1/{\sigma }_t+\mathbf{i}\omega )}\sqrt{\frac{2\lambda }{{\sigma }_t}}\frac{1}{(1/{\sigma }_t-\mathbf{i}\omega )}$$

$$K(x,x^{\prime },\tau )=K_{\mathrm{s}}(x,x^{\prime })h(\tau )=e^{-(x-x^{\prime }{)}^2/{\sigma }_x}\lambda e^{-|\tau |/{\sigma }_t}$$

$$W(\mathbf{i}\omega )=\sqrt{\frac{2\lambda }{{\sigma }_t}}(1/{\sigma }_t+\mathbf{i}\omega ),$$

Full text

1 Introduction

The “Big-Data” and “Machine-Learning” era we are living in during the last few decades has been supplied by the exponential growth of research fields like statistics and optimization. Within these areas, function estimation problems play an important role and many different regression approaches have been developed in recent years. In particular, in the Bayesian estimation context, Gaussian process (GP) methods [1], also known as kriging [2], have become the standard approach [3, 4] in application realms such as robotic networks, biomedicine, and system identification [5, 6, 7].

In the classical GP framework the process is assumed to be static so that only spatial locations are seen as input variables. However, due to the heavy computational requirements, characterized by a cubical growth in the number of input data, many efficient approaches have been developed. Some of these rely e.g., on the notion of pseudo input locations [8, 9, 10], the use of matrix factorizations [11] and approximations of the kernel function [12, 13] through the Nyström method or greedy techniques [14, 15, 16]. Different is the work [17] where the authors consider a state-space approach.

Recent research has instead focused on the use of such the classical methods in dynamical contexts. In fact, to capture many interesting time-varying phenomena, like wind and ocean currents, it is necessary to extend the methodology to the class of spatio-temporal processes. The simplest approach is to interpret time just as an additional input feature [4]. However, in dynamical scenarios this approach has important practical limitations mainly due to: i) heavy memory and computational requirements; and ii) the non iterative nature of the methodology. Indeed, the classical paradigm, being tailored for static processes, relies on batch implementations where data are processed at once, after they have been collected. In the dynamical context, two explored approaches to cope with computational complexity are sparse approximations [4, 18] or finite memory implementations [19, 20], often based on truncated observations. This paper instead takes inspiration from a different idea related to the use of Kalman filter [21]. In this context, the works [22, 23], whose inceptive idea can be traced back to [1], indeed focus on building state-space representations for Gaussian processes amenable to Kalman filtering. More specifically, [22] presents a preliminary result which applies to temporal GP regression models while [23] extends the approach to spatio-temporal GPs introducing also infinite dimensional state-space models.

Following the line of research inspired by [23], this work focuses on the practical implementability of dynamic GP regression. Our contributions rely on two assumptions. First, we assume that the kernel process is separable in space/time, with a rational temporal power spectral density. Then, we also assume that measurements are collected only on a finite set of fixed spatial locations (this hypothesis will be however removed in the second part of the paper). Our novel contributions can be summarized as follows:

• •

in [23] only approximated filtering schemes are proposed to deal with infinite dimensional operators, e.g. based on eigenfunctions expansions of the operators which govern the stochastic dynamics. Our estimation procedure is instead exact, i.e. it returns the exact minimum variance estimate on any spatio-temporal prediction location. This result comes from a novel result which we refer to as Kalman representer theorem. In the static scenario, classical representer theorems state that the optimizers of a wide class of variational problems admit a finite-dimensional representation [24, 25]. In particular, in the case of regularization networks the function estimate belongs to the subspace generated by the kernel sections centred on the observed input locations, with coefficients given by a linear output transformation [26]. The Kalman representer theorem here derived shows that in a dynamic scenario the estimate at instant $t_{k}$ is still the combination of the spatial kernel sections centred on the spatial locations but with time-varying coefficients which now depend linearly on the Kalman filter state;

• •

process stationarity is a key assumption in [23]. In this paper instead, inside the rich class of space-time separable kernel functions, no restriction is imposed on the nature of the spatial kernel. As for the temporal kernel, the only requirement is that it admits a state-space description which can be also time-varying. This restriction is really mild: in practice, due to the universal approximation properties of rational functions, one can just increase the state space dimension to approximate with arbitrary accuracy any temporal spectrum;

• •

we obtain a dynamic regression procedure which is computationally efficient. In particular, the complexity scales with the cube of the number of distinct measurement locations and only linearly on the number of prediction locations. Conversely, a naive implementation of a GP estimator exploiting a finite buffer, i.e. which uses only the most recent measurements, is not only unable to provide an exact solution but also has a complexity per iteration which grows cubically in terms of the buffer size (which could be much larger than the number of distinct measurement locations);

• •

we address also the situation where the sampling locations set is adaptive and changes over time. This set-up is interesting in many applications like aerial vehicle wind estimation and multi-robots exploration, where it is possible to keep in memory only a finite set of locations and measurements due to storage capacity limits. When memory limits are hit, old locations can be discarded (following a policy depending on the specific application) and the state-space can be accordingly modified. We design a new approach to perform such operations. In addition, we show that, after any change in the sampling set, if no other perturbations occur, our suboptimal estimate converges (with an exponential rate) to the optimal minimum variance estimate (obtainable only storing all the past measurements).

The remainder of the paper is organized as follows. In Section 2 we recall all the necessary preliminaries on GP estimation, Kalman filtering and spectral factorization theory. In Section 3 we formulate the problem at hand and present the necessary assumptions. In Section 4 we propose the solution to the estimation problem over the finite set of sampling locations. In Section 5 we extend the result to the prediction problem at any spatio-temporal location. In Section 6 we discuss on the computational complexity. In Section 7 we address the problem of time-varying sampling locations. In Section 8 we present some compelling simulations. Finally, in Section 9 we draw some concluding remarks. All the technical proofs are collected in the Appendix.

2 Preliminaries

In the following we recall all the necessary preliminaries on GP regression, Kalman filtering and spectral factorization of random processes, respectively.

2.1 GP regression

Let $f:\mathcal{A}\mapsto{\mathbb{R}}$ be a zero-mean Gaussian field with covariance, also called kernel, $K:\mathcal{A}\times\mathcal{A}\mapsto{\mathbb{R}}$, where $\mathcal{A}$ is a compact set. Assume to have a set of $N\in{\mathbb{N}}_{>0}$ noisy measurements of the form

[TABLE]

where $v_{i}$ is a zero-mean Gaussian noise with variance $\sigma^{2}$, i.e. $v_{i}\sim\mathcal{N}(0,\sigma^{2})$, independent from the unknown function. Given the data set of input locations $\{a_{i},y_{i}\}_{i=1}^{N}$, it is known [27, 3] that the estimate $\widehat{f}$ of $f$ is a linear combination of the kernel sections $K(a_{i},\cdot)$, i.e., the kernel sampled in the values corresponding to the available input locations. In particular, for any $a\in\mathcal{A}$, it holds that

[TABLE]

with expansion coefficients $c_{i}$ obtained as

[TABLE]

where $I$ denotes the identity matrix of suitable size and $\bar{K}$ is element-wise defined as $[\bar{K}]_{ij}:=K(a_{i},a_{j})$. Finally, the posterior variance of $\widehat{f}(a)$ evaluated at the generic location $a\in\mathcal{A}$ is given by

[TABLE]

Clearly, because of the matrix inversion in both (2) and (3), the method scales as $\mathcal{O}(N^{3})$. Moreover, in real-time applications, where a certain number of measurements are collected at each iteration, all the past measurements must be kept in memory. Thus, the method is more suitable for a batch and almost static implementation rather than for an iterative time-varying one.

Remark 1** (Spatio-temporal processes)**

In the following we consider spatio-temporal processes. Conversely to classical Gaussian processes, where the “location” $a$ usually denotes a spatial variable, in spatio-temporal processes $a$ represents both time and space. Hence, without loss of generality, we can write $f(a)=f(x,t)$. Accordingly, the domain $\mathcal{A}$ can be decomposed as $\mathcal{A}=\mathcal{X}\times\mathbb{R}_{+}$, with $\mathcal{X}$ and ${\mathbb{R}}_{+}$ denoting the spatial and temporal domain, respectively.

2.2 Kalman Filtering

Consider the following discrete-time dynamical system

[TABLE]

where, at each iteration $k$, $s(k)\in{\mathbb{R}}^{n}$ is the state vector, $y(k)\in{\mathbb{R}}^{m}$ is the output vector, $w(k)\in{\mathbb{R}}^{n}$ and $v(k)\in{\mathbb{R}}^{m}$ are i.i.d. zero-mean Gaussian random vectors with covariance matrices $Q\geq 0$ and $R>0$, respectively. $A(k)\in{\mathbb{R}}^{n\times n}$ and $C(k)\in{\mathbb{R}}^{m\times n}$ are the time-varying state and output matrices, respectively. As commonly done, we assume both process and measurement noises to be uncorrelated with respect to each other, i.e. ${\mathbb{E}}\left[w(k)^{T}v(h)\right]=0\ \forall_{k,h}$. Without loss of generality, we also assume the initial condition $s(0)$ being drawn from a Gaussian distribution with zero mean and covariance $\Sigma_{0}$, i.e., $s(0)\sim\mathcal{N}(0,\Sigma_{0})$.

The Kalman Filter [28] applied to the discrete-time linear state-space system (4) is described by the following recursive equations

[TABLE]

where $\widehat{s}(k|k)$ and $\Sigma(k|k)$ represent the filtered estimate of the state and the posterior error covariance, respectively; $\widehat{s}(k+1|k)$ and $\Sigma(k+1|k)$ represent the (one step) predicted state estimate and error covariance, respectively; $L(k+1)$ is the Kalman gain; finally, the filter is initialized assuming $\widehat{s}(0|-1)=\mathbb{E}[s(0)]=0$ and $\Sigma(0|-1)=Cov[s(0)]=\Sigma_{0}$.

We recall that, under the assumptions of normally distributed noises and perfect model knowledge, the Kalman filter is optimal, in mean square sense. Then, Eqs. (5) return the minimum mean square error estimate of the state, which corresponds to

[TABLE]

that is, the estimate of the state given all the measurements up to the $k$-th one. Moreover, in view of the Markovianity (memory-less of the system) property of the state, it holds that

[TABLE]

that is, the previous state and the last measurement represent the sufficient statistic to compute the optimal estimate of the state at the current time instant.

Finally, it is well known [29] that if the state and output matrices are constant, i.e. $A(k)=A$ and $C(k)=C$, under the additional hypothesis of stabilizability of the pair $(A,Q)$ and detectability of the pair $(A,C)$, the estimation error covariance of the Kalman filter converges to a unique value from any initial condition.

2.3 Spectral factorization of random processes

Here we want to show how a specific class of processes admits an equivalent exact state-space representation.

Consider a stationary random process $f(t)$ with covariance $h(\tau)$, $\tau=t-t^{\prime}$. Thanks to the Wiener-Khinchin theorem, it is known that the power spectral density (PSD) of the process is equal to the Fourier transform of its covariance $h$, i.e.,

[TABLE]

Moreover, in the particular case when $S(\omega)=S_{\rm r}(\omega)$ is rational of order $2r$, thanks to spectral factorization [30], its PSD can be rewritten as $S_{\rm r}(\omega)=W(\mathbf{i}\omega)W(-\mathbf{i}\omega)$ with

[TABLE]

where $\mathbf{i}$ denotes the imaginary unit. Finally, from realization theory, we have that rational functions of the form (6) are in correspondence to the equivalent continuous time state-space representation [31] (companion form) given by

[TABLE]

where $w(t)\sim\mathcal{N}(0,I)$, the model matrices are equal to

[TABLE]

and the initial state is $s(0)\sim\mathcal{N}(0,\Sigma_{0})$, with $\Sigma_{0}$ computed as solution of the Lyapunov equation $FX+XF^{T}+GG^{T}=0$.

3 Problem Formulation

Here, we formally state the main problem at hand and we introduce the necessary assumptions. As clarified later, the main assumption restricts our analysis on a particular yet sufficiently rich class of kernel functions separable in space and time. We already stress that, conversely to [32, 23], where the authors, dealing with infinite dimensional state-space systems, must resort to approximated approaches in order to practically implement their solutions, we develop an exact methodology without requiring stationarity of the kernel.

Consider a function $f:\mathcal{X}\times{\mathbb{R}}_{+}\rightarrow{\mathbb{R}}$ modeled as a zero-mean Gaussian Process with covariance $K$. Let $\mathcal{X}$ be any compact set. We define a finite dimensional subset ${\mathcal{I}}\subseteq\mathcal{X}$ consisting of a collection of given spatial input locations as follows

Definition 2** (Input Location Space)**

Consider the set $\mathcal{X}$. We denote with ${\mathcal{I}}\subseteq\mathcal{X}$ a finite collection of points containing $M$ locations from $\mathcal{X}$, i.e.

[TABLE]

As suggested by Definition 2, ${\mathcal{I}}$ represents our “observable” location space. Precisely, to consider the most general case, we assume to be able to collect noisy measurements of the form

[TABLE]

at non-uniformly distributed discrete-time instants $t_{k}$ only from a time-varying subset of spatial locations contained in ${\mathcal{I}}$, namely $x_{i}\in\mathcal{M}(k)\subseteq{\mathcal{I}}$ ($|\mathcal{M}(k)|=M_{k}$). To help the reader’s understanding, Figure 1 shows an illustrative representation of the considered non-uniform spatio-temporal sampling and measurements collection process.

The problem we want to solve is that of estimating $f$ over the entire “partially observable” domain $\mathcal{X}$, exploiting measurements coming from the “observable” subset ${\mathcal{I}}$. The problem could arise in diverse applications, e.g., in weather forecasting where, given a small number of weather stations which are able to collect measurements at certain discrete time instants, the goal is to estimate the weather conditions on a larger area.

To state our solution, we restrict the analysis on the following specific, yet sufficiently rich, class of kernel functions separable in space and time.

Assumption 3** (Generating Kernel properties)**

The kernel function $K$, covariance of the Gaussian process $f(x,t)$, is separable in time and space and stationary in time, namely,

[TABLE]

In addition, the power spectral density $S_{\rm r}(\omega)=W(\mathbf{i}\omega)W(-\mathbf{i}\omega)$ of $h(\tau)$ is a rational function of order $2r$, where $W(\mathbf{i}\omega)$ is like in (6).

We stress the fact that, differently from [32, 23], in Assumption 3 we do not require space stationarity of the process but just separability and time stationarity of the kernel function.

Our solution consists of two steps: first we show how to estimate the process $f$ over ${\mathcal{I}}$ (Section 4). Then, we extend our result to obtain a prediction of the process outside ${\mathcal{I}}$ (Section 5). Precisely, we show how our first solution can be exploited to reconstruct $f$ on any arbitrary spatio-temporal location $(x,t)\in\mathcal{X}\times{\mathbb{R}}_{+}$.

4 Kalman Regression on ${\mathcal{I}}$

Here we formally show how to built an exact state space representation for a certain class of GPs and we bridge GP regression and Kalman filtering over the observable finite-dimensional space ${\mathcal{I}}$, providing a clear and systematic methodology to implement the filter.

To implement the Kalman Eqs. (5), the first step is to build a state-space representation for the Gaussian process $f$. In particular, we are interested in reconstructing $f$ over the “observable” ${\mathcal{I}}$. To compactly represent the process over ${\mathcal{I}}$, it is convenient to define the vector

[TABLE]

The next proposition exploits Assumption 3 and the state-space realization for rational PSD given in (7) to show that the process $\mathbf{f}(t)$, admits an equivalent exact continuous-time state-space representation.

Proposition 4** (Equivalent CT-SS representation for $\mathbf{f}(t)$)**

Consider the process $\mathbf{f}(t):{\mathcal{I}}\times{\mathbb{R}}_{+}\mapsto{\mathbb{R}}^{M}$ generated by the spatio-temporal kernel $K$ satisfying Assumption 3. Let the triplet $(F,G,H)$ be a state-space representation for $S_{\rm r}(\omega)$ as described in Section 2.3. Then, $\mathbf{f}(t)$ admits the following strictly proper state-space representation

[TABLE]

where

• •

$\mathbf{z}(t)\coloneqq\left[z_{1}(t),\ldots,z_{M}(t)\right]^{T}\in{\mathbb{R}}^{M}$;

• •

$\bar{K}_{\rm s}\in{\mathbb{R}}^{M\times M}$* is obtained sampling the spatial kernel $K_{\rm s}$ over ${\mathcal{I}}$, i.e., $[\bar{K}_{\rm s}]_{ij}=K_{\rm s}(x_{i},x_{j})$, $x_{i},x_{j}\in{\mathcal{I}}$;*

• •

for $i\in\{1,\dots,M\}$, $w_{i}(t)$ and $s_{i}(0)$ are defined as for (7).

Observe that the subsystems $S_{i}$ in (9) are independent one from each other in the sense that one can easily verify that ${\mathbb{E}}\left[(s_{i}(t))^{T}(s_{j}(t))\right]=0\ \forall t,\forall i\neq j$. Basically, Proposition 4 states that, for each location $x_{i}\in{\mathcal{I}}$, the time evolution of $f(x_{i},t)$ admits a state space representation given by the system $S_{i}$ in Eq. (9). Then, these state space representations are “combined” through the sampled spatial kernel $\bar{K}_{\rm s}$ to build a representation for the overall process $\mathbf{f}(t)$. Figure 2 shows an illustrative representation of the process and measurements formation schemes for the case of five input locations.

Now, let be $\mathbf{s}(t)=\left[s_{1}^{T}(t),\ldots,s_{M}^{T}(t)\right]^{T}\in{\mathbb{R}}^{rM}$ and $\mathbf{w}(t)=\left[w_{1}(t),\ldots,w_{M}(t)\right]^{T}\in{\mathbb{R}}^{M}$, then (9) can be re-written in a more compact form as

[TABLE]

Observe that Eq. (10) gives a continuous-time state-space representation for the process. However, the Kalman Eqs. (5) as well as the measurements in (8) are expressed in discrete time. Thus, it is convenient to discretize the system (10) as

[TABLE]

where the measurements Eq. (8) has been already embedded in the model and, defining $T_{k}:=t_{k}-t_{k-1}$,111Here, with a little abuse of notation, with $k$ we refer to estimates at discrete time instant $t_{k}$. The same holds for the measurements which are collected according to Eq. (8). This is done to explicitly highlight the fact that at $t_{k}$ we perform the $k$-th iteration of the Kalman filter. we have

• •

$A(k)=\exp^{(I\otimes F)T_{k}}\in{\mathbb{R}}^{rM\times rM}$;

• •

$\mathbf{w}(k)\in{\mathbb{R}}^{rM}$ is a zero-mean random Gaussian noise with variance $Q(k)=I\otimes\bar{Q}(k)$, where

[TABLE]

• •

$\mathbf{y}(k)=\left[y_{i_{1}}(k),\ldots,y_{i_{M_{k}}}(k)\right]^{T}$ with $i_{1},\ldots,i_{M_{k}}$ the indexes identifying the locations which are available at $t_{k}$ according to $\mathcal{M}(k)\subseteq{\mathcal{I}}$;

• •

$\mathbf{v}(k)=\left[v_{i_{1}}(k),\ldots,v_{i_{M_{k}}}(k)\right]^{T}$,   $\mathbf{v}(k)\sim\mathcal{N}(0,R(k))$, with $R(k):=\sigma^{2}I\in{\mathbb{R}}^{M_{k}\times M_{k}}$;

• •

$C(k)=I_{k}\bar{K}_{\rm s}^{1/2}\left(I\otimes H\right)$ with $I_{k}\in\{0,1\}^{M_{k}\times M}$ the matrix selecting the input locations which are available at $t_{k}$ according to the set $\mathcal{M}(k)\subseteq{\mathcal{I}}$ (see Figure 2).

Before proceeding, it is important to remark the all the $A(k)$ are stable. Indeed, since the matrix $F$ derives from a state-space representation of a stationary power spectral density, $F$ is stable and, in turn, this implies the stability of $A(k)$. It easily follows that the pairs $(A(k),C(k))$ are detectable and the pairs $(A(k),Q(k))$ are stabilizable.

Now, based on the information gathered at discrete time instants $t_{k}$ and thank to the equivalent discrete-time representation (11), we show how to reconstruct the minimum variance estimate $\widehat{\mathbf{f}}(t)$ of $\mathbf{f}(t)$ at time instant $t\in{\mathbb{R}}_{+}$, $t\geq t_{k}$, $t_{k}$ being the last available sampling instant, defined as

[TABLE]

Proposition 5** (Kalman-based Regression on ${\mathcal{I}}$)**

Let Assumption 3 holds and assume to collect measurements of the form (8) at $t_{k}$ from a subset $\mathcal{M}(k)\subseteq{\mathcal{I}}$. Then, the estimate $\widehat{\mathbf{f}}(t)$ of $\mathbf{f}(t)$ with corresponding error covariance $\Sigma^{\mathbf{f}}(t)$ for any $t\in{\mathbb{R}}_{+}$, $t\geq t_{k}$, is computed according to Algorithm 1.

Remark 6

Note that, while the variable $\widehat{\mathbf{s}}(t)$ can be viewed as the state evolution of the temporal part of the process, it has no direct or easy physical interpretation with respect to the entire process $\mathbf{f}(t)$. Nevertheless, it contains all the information gathered from the collected measurements.

Figure 3 shows a block diagram of the estimation scheme corresponding to Algorithm 1. The block “State Regression” implements what the described in the “if-then-else” part of Algorithm 1 (lines 7$\div$27). Observe that measurements arrive only at $t_{k}$. Moreover, it is worth noting that in the case of uniform sampling from a constant set of locations, e.g., the entire set ${\mathcal{I}}$, all the matrices $A(k)$, $C(k)$, $Q(k)$ and $R(k)$ becomes constant. Thus, thanks to standard results on Kalman filtering [29], the filter gain $L(k)$ converges to a constant value which can be pre-computed offline. In this case the filtering corresponds to a static matrix multiplication thus alleviating the computational burden.

In order to help the reader’s intuition in the building process of the presented estimation procedure, we now present an example.

Example 7

Consider the exponential time kernel $h(\tau)$

[TABLE]

satisfying Assumption 3 since its PSD $S_{r}$ is equal to

[TABLE]

which is rational of order $2$.

Now, consider a zero-mean Gaussian process $f(x,t)$ with covariance

[TABLE]

that is, a Gaussian spatial kernel and an exponential time kernel. Moreover, for simplicity, we assume to collect measurements from the entire location set ${\mathcal{I}}$ at periodic time instants, i.e., $t_{k}=kT$. Thanks to Proposition 4, since $K$ satisfies Assumption 3, $\mathbf{f}(t)$ admits a state space representation. In particular, given $S_{r}$ as in (15) with

[TABLE]

it is easy to see the state-space model matrices are equal to

[TABLE]

while the matrix $\bar{K}_{\rm s}^{1/2}$ is computed as the Cholesky factorization of the sampled kernel $\bar{K}_{\rm s}$. Thanks to this, the discrete-time state-space representation follows from Eq. (11) with

[TABLE]

Observe that the state-space modeling is exact since $h(\cdot)$ satisfies Assumption 3. To conclude, thanks to the matrices $A$, $Q$ and $C$ above, Kalman filtering can be applied.

5 Kalman Regression on $\mathcal{X}$

In Section 4 we showed how to build an estimate $\widehat{\mathbf{f}}$ of the process $f$ over the observable finite-dimensional set ${\mathcal{I}}$. Here, we are interested in extending the result of Proposition 5 to build the minimum variance estimate

[TABLE]

of the process $f$ on any desired spatio-temporal location $(x,t)\in\mathcal{X}\times{\mathbb{R}}_{+}$, $t\geq t_{k}$, being $t_{k}$ the last time instant where measurements have been collected, see Figure 1.

To state our result we first introduce the following additional symbols

[TABLE]

where, with a slight abuse of notation, $\bar{K}_{\rm s}(\cdot,\cdot)$ denotes the space kernel $K_{\rm s}$ evaluated in all the locations contained in its arguments. More in details we have that

[TABLE]

while the $(i,j)-$th element of $\bar{K}_{\rm s}({\mathcal{I}},{\mathcal{I}})\in{\mathbb{R}}^{M\times M}$ is equal to

[TABLE]

Proposition 8** (Kalman Representer Theorem on $\mathcal{X}$)**

Consider the process $f:\ \mathcal{X}\times{\mathbb{R}}_{+}\mapsto{\mathbb{R}}$ generated by the separable kernel $K(x,x^{\prime},t,t^{\prime})=K_{\rm s}(x,x^{\prime})h(\tau)$ satisfying Assumption 3. Then, at the generic instant $t\geq t_{k}$, being $t_{k}$ the last available sampling instant, the estimate $\widehat{f}(x,t)$ of $f(x,t)$ is given by

[TABLE]

while its posterior variance is given by

[TABLE]

The result of the above proposition extends that of Section 4 to any additional input location $x\in\mathcal{X}$ by providing expressions for the estimate and the corresponding posterior variance. However, in view of Section 7, it is useful to introduce Eq. (8) in vector form in order to simultaneously compute the expression for the posterior variance for the joint vector $[\mathbf{f}(t)^{T}\ f(x,t)]^{T}$ given all the available measurements. As for Eq. (8) the expression follows from Lemma 1 in [33] and reads as

[TABLE]

In Proposition 8, observe that $\widehat{\mathbf{f}}(t)$ and $\Sigma^{\mathbf{f}}(t)$ are the estimated process and covariance as returned from Algorithm 1. Hence, if $t=t_{k}$ then $f(x,t)$ makes use of $\widehat{\mathbf{f}}(t)$ computed with the estimated state, output of the Kalman equations. Differently, if $t>t_{k}$, then the open-loop predicted state is used (see Section 4).

Proposition 8 states that the output of the Kalman filter captures all the necessary information, contained in the measurements, to estimate the entire process. Indeed, $\widehat{\mathbf{f}}$ is a sufficient statistic to reconstruct $f$ over any spatio-temporal location $(x,t)\in\mathcal{X}\times{\mathbb{R}}_{+}$. Hence, this result can be regarded as a Kalman Representer Theorem for GPs generated by separable kernels. Precisely it states that the current estimate for the entire process $f$ is captured by $M$ basis functions, i.e., the spatial kernel sections $\bar{K}_{\rm s}(x,{\mathcal{I}})$, being $M$ the number of distinct spatial input locations contained in ${\mathcal{I}}$, whose coefficients are computed by means of a Kalman filter with state dimension $M\times r$. Moreover, Algorithm 1 together with Proposition 8 outline an exact methodology to implement the filtering procedure, illustratively represented in Figure 4.

Remark 9** (On non stationary time kernels)**

*So far, both the results of Propositions 5 and 8 rely on Assumption 3 of stationary time kernel $h(\tau)$. However, it is worth stressing the fact that the results hold for a more general class of time kernel functions. In particular, conditioned to the existence of a (possibly) time-varying state-space representation for the time kernel, the results extend both to kernels whose PSD are not rational as well as to non-stationary kernels. Indeed, given a state-space representation for $h(\cdot,\cdot)$, the estimation procedure outlined in Algorithm 1 is exact and returns the minimum variance estimate $\widehat{\mathbf{f}}$ and $\mathbf{f}$ as defined in (12). Consequently, the result of Proposition 8 seamlessly holds.

According to the statement of Proposition 8, in Appendix 10.3 we report the proof for the case of stationary kernels. However, for the sake of completeness, we refer the interested reader to Appendix 10.3.1 where we briefly outline the proof for the case of non-stationary kernels.*

6 Computational Complexity

One of the underlying reasons to build a recursive filtering procedure is to break down the computational complexity induced by the classical GP approach which grows cubically with the total number of collected measurements. Interestingly, thanks to its recursive implementation, the computational complexity of the proposed scheme scales as

[TABLE]

where $r$ is the order of a single state-space model $S_{i}$ in (9), $M=|{\mathcal{I}}|$, $M_{k}=|\mathcal{M}(k)|$ and $P$ is the number of “prediction” locations $x\in\mathcal{X}$ where we want to extend the estimate.

Conversely, the classical GP approach is characterized by a complexity of order

[TABLE]

Hence, in a real-time implementation, the computational cost per iteration for the proposed scheme scales linearly with the model complexity $r$. Conversely the cost for the classical GP regression implementation grows cubically with the total number of collected measurements. For further analysis of the computational complexity characterizing the proposed regression scheme against that of the classical GP approach we refer the interested reader to our prelimiary work [34].

7 Adaptive Input Location Space ${\mathcal{I}}(k)$

In Proposition 5 we have assumed that the input locations always fall in a fixed set ${\mathcal{I}}$. In other words, even if the measurements collected at $t_{k}$ might come from a time-varying subset $\mathcal{M}(k)$, one has $\mathcal{M}(k)\subseteq{\mathcal{I}}\ \forall k$. In this section we will remove this constraint to allow measurements of $f$ to be collected over an adaptive set of input locations, i.e., the location set becomes now a function of the time ${\mathcal{I}}(k)$. See Figure 5 for an illustrative representation of the time evolution of ${\mathcal{I}}(k)$. This scenario is important in many applications. Consider for instance a group of air vehicles whose aim is to estimate meteorological phenomena from punctual measurements of, e.g., cumulus-type clouds [35, 36]. As the vehicles proceed in time and space, they collect measurements coming from new input locations. Accordingly, starting from the first location, they might want to “expand” their state-space in order to better estimate the process. Nevertheless, at some point, the vehicles could reach their memory capacity and should start to discard previously collected input locations. This calls for the development of suitable strategies both to “expand” and to “contract” the state-space.

It is worth already pointing out that our strategy to reconstruct $f$ from measurements on the time-varying location set ${\mathcal{I}}(k)$ is sub-optimal for two reasons: i) when a new input location is visited, the optimal state estimate extension would require to run a Kalman filter reprocessing all the measurements collected from the beginning of the experiment. This is because the underlying dynamical model corresponding to the enlarged location set is different and so is the state evolution. It is thus necessary to design a sub-optimal, yet sensible, strategy to extend both the state-space model and the estimate without reprocessing any past data; ii) when input locations are discarded for memory constraints some information is inevitably lost.

As in the previous Sections, the sampling instants where measurements are collected (possibly from both old and newly visited input locations) are denoted as $t_{\ell}$ for $\ell\in{\mathbb{Z}}_{+}$. For ease of notation we sometimes refer to $t_{\ell}$ simply as $\ell$. Now, without loss of generality, assume that, at the generic time instant $k-1$, the location set is given by

[TABLE]

and that, at time instant $k$, just one new input location222For clarity of presentation we assume to visit just one new location but the procedure can be seamlessly extended to the case of multiple locations. $x_{M}$ is visited, thus

[TABLE]

According to (8), the measurement taken on $x_{M}$ at $k$ is denoted as

[TABLE]

Observe that, in general, at time $k$, in addition to $y_{M}(k)$, the system might collect measurements taken also on some old input locations of ${\mathcal{I}}(k-1)$; more precisely, for $s\leq M-1$, let

[TABLE]

denote the set of input locations which are visited together with $x_{M}$ at time $k$, and let $\left\{y_{i_{1}}(k),\ldots,y_{i_{s}}(k)\right\}$ be the corresponding measurements taken.

Before describing how the measurements set $\left\{y_{i_{1}}(k),\ldots,y_{i_{s}}(k),y_{M}(k)\right\}$ is used to properly provide an estimate of $f$ over the augmented input location set ${\mathcal{I}}(k)$, we introduce the following notation. First, since, as previously stressed, the proposed procedure is only sub-optimal, we use $\widetilde{\cdot}$ instead of the more common $\widehat{\cdot}$ to denote the returned estimates. Secondly, given a generic set of locations ${\mathcal{I}}^{\prime}=\{x_{1},\ldots,x_{\ell}\}$ we define $\mathbf{f}_{{\mathcal{I}}^{\prime}}(\cdot):=[f(x_{1},\cdot),\ldots,f(x_{\ell},\cdot)]^{T}$; similarly for the state we use the symbol $\mathbf{s}_{{\mathcal{I}}^{\prime}}(\cdot)$.

Now, we assume that at time $k-1$, estimates of $\mathbf{f}_{{\mathcal{I}}(k-1)}$ and $\mathbf{s}_{{\mathcal{I}}(k-1)}$, obtained processing all the measurements up to $k-1$, are available. Consistently with the adopted notations, we denote them, respectively, as

[TABLE]

accordingly, we denote the corresponding covariance matrices by

[TABLE]

Next we show how the new measurements $\left\{y_{i_{1}}(k),\ldots,y_{i_{s}}(k),y_{M}(k)\right\}$ can be properly exploited to update $\widetilde{\mathbf{f}}_{{\mathcal{I}}(k-1)}(k-1)$, $\widetilde{\mathbf{s}}_{{\mathcal{I}}(k-1)}(k-1)$, $\widetilde{\Sigma}^{\mathbf{f}}_{{\mathcal{I}}(k-1)}(k-1)$, $\widetilde{\Sigma}^{\mathbf{s}}_{{\mathcal{I}}(k-1)}(k-1)$ into the enlarged statistics

[TABLE]

Specifically, the procedure we propose consists of the following five major steps which are performed in order:

estimate update on ${\mathcal{I}}(k-1)$: in this step, we exploit Algorithm 1 which uses the measurements collected on the old input locations, i.e., $\left\{y_{i_{1}}(k),\ldots,y_{i_{s}}(k)\right\}$, and the statistics $\widetilde{\mathbf{s}}_{{\mathcal{I}}(k-1)}(k-1)$ and $\widetilde{\Sigma}^{\mathbf{s}}_{{\mathcal{I}}(k-1)}(k-1)$ (which play the role respectively of $\widehat{\mathbf{s}}(k-1|k-1)$ and $\Sigma(k-1|k-1)$) to compute the updated estimate of $f$, sampled over ${\mathcal{I}}(k-1)$, and its corresponding covariance, that is,

[TABLE] 2. 2.

space prediction on $x_{M}$: based on $\widetilde{\mathbf{f}}_{{\mathcal{I}}(k-1)}(k)$ and by leveraging Proposition 8, it is possible to predict $f$ over the new location $x_{M}$; specifically from (8), we have

[TABLE] 3. 3.

estimate update on ${\mathcal{I}}(k)$: First, the estimate $\widetilde{f}(x_{M},k)$ is used to build the augmented vector

[TABLE]

whose corresponding covariance matrix $\bar{\Sigma}^{\mathbf{f}}_{{\mathcal{I}}(k)}(k)\in{\mathbb{R}}^{M\times M}$ is computed according to formula in (5), by replacing $\Sigma^{\mathbf{f}}(t)$ with $\widetilde{\Sigma}^{\mathbf{f}}_{{\mathcal{I}}(k-1)}(k)$, $\Gamma_{x}$ with $\Gamma_{x_{M}}=h(0)\bar{K}_{\rm s}(x_{M},{\mathcal{I}}(k-1))$, $V_{x}$ with $V_{x_{M}}=h(0)\bar{K}_{\rm s}(x_{M},x_{M})$ and where $V_{\mathbf{f}}=h(0)\bar{K}_{\rm s}({\mathcal{I}}(k-1),{\mathcal{I}}(k-1))$.

Secondly, $\bar{\mathbf{f}}_{{\mathcal{I}}(k)}(k)$ and $\bar{\Sigma}^{\bf f}_{{\mathcal{I}}(k)}(k)$ are updated to obtain the estimates $\widetilde{\mathbf{f}}_{{\mathcal{I}}(k)}(k)$ and $\widetilde{\Sigma}^{\mathbf{f}}_{{\mathcal{I}}(k)}(k)$, by a Kalman-like correction step which exploits the measurement $y_{M}(k)$; precisely, observing that, Eq. (21) can be rewritten as

[TABLE]

where $\Upsilon:=[0\ \cdots\ 0\ 1]$, we define the Kalman gain as

[TABLE]

obtaining that

[TABLE] 4. 4.

contraction: if memory requirements are hit, some input location must be discarded. Consequently, one has to update ${\mathcal{I}}(k)$, and contract the vector of function estimate $\widetilde{\mathbf{f}}_{{\mathcal{I}}(k)}(k)$ and the corresponding covariance matrix $\widetilde{\Sigma}^{\mathbf{f}}_{{\mathcal{I}}(k)}(k)$ (line 20-22). It is worth stressing that the particular heuristic to choose the input location to discard depends on the specific application. 5. 5.

state reconstruction: finally, observe that, in the previous steps, we have computed the estimates $\widetilde{\mathbf{f}}_{{\mathcal{I}}(k)}(k)$, $\widetilde{\Sigma}^{\mathbf{f}}_{{\mathcal{I}}(k)}(k)$, without having reconstructed the estimates $\widetilde{\mathbf{s}}_{{\mathcal{I}}(k)}(k)$, $\widetilde{\Sigma}^{\mathbf{s}}_{{\mathcal{I}}(k)}(k)$. However the state statistics are necessary in our machinery since they represent the required inputs to run Algorithm 1, either between two subsequent measurements time instants, i.e., when $t\in(t_{k},t_{k+1})$ or when new measurements are taken, i.e, when $t=t_{k}$.

Observe that, maintaining optimality would require to go back in time and to restart the filtering procedure reprocessing all the measurements, but this becomes unfeasible as the number of measurements increases. The alternative idea we propose is to assume that all the information collected up to $k$ has been actually collected at instant $k$ according to a particular “static” virtual measurement model of the form

[TABLE]

which is consistent with $\widetilde{\mathbf{f}}_{{\mathcal{I}}(k)}(k)$ and $\widetilde{\Sigma}^{\mathbf{f}}_{{\mathcal{I}}(k)}(k)$, where $\check{\mathbf{y}}$ is the vector containing the virtual measurements while $\check{\mathbf{v}}$ is the virtual noise vector. In particular, we show that $\check{\mathbf{y}}$ and $\check{\Sigma}^{\mathbf{v}}$ are uniquely determined. Indeed, recall that in absence of measurements, $\Sigma_{0}$ is the covariance of any state, as defined by the Lyapunov equation in the first lines of Algorithm 1. According to the state-space model of Eq. (10), $\mathbf{f}_{{\mathcal{I}}(k)}(k)=C(k)\mathbf{s}_{{\mathcal{I}}(k)}(k)$, $C(k)=\bar{K}_{\rm s}({\mathcal{I}}(k),{\mathcal{I}}(k))^{1/2}(I\otimes H)$, with prior variances respectively equal to

[TABLE]

Then, to obtain $\widetilde{\mathbf{s}}_{{\mathcal{I}}(k)}(k)$ and its corresponding covariance $\widetilde{\Sigma}^{\mathbf{s}}_{{\mathcal{I}}(k)}(k)$, it is necessary to compute $\check{\mathbf{y}}$ and $\check{\Sigma}^{\mathbf{v}}$. By exploiting the classical formulas of Eqs. (1)$\div$(3) combined with the matrix inversion lemma we have that

[TABLE]

Then, one obtains

[TABLE]

Finally, one has

[TABLE]

From equation above, notice that the explicit computation of $\check{\mathbf{y}}$ is not necessary and only $\check{\Sigma}^{\mathbf{v}}$ is indeed needed.

The previous five steps are summarized in the algorithmic description provided in Algorithm 2. Clearly, if no new input locations are visited, then Algorithm 2 coincides with Algorithm 1.

Next we establish some properties of Algorithm 2. We start by stressing that that steps 1), 2) and 3) may preserve optimality. This fact is made precise in the following Corollary.

Corollary 10

For a generic time instant $j$, let the whole data set up to $j$ be defined as

[TABLE]

In addition, define the optimal minimum variance estimate and corresponding error covariance at $j$ over ${\mathcal{I}}(j)$ as

[TABLE]

and assume that

[TABLE]

Then, it holds that

[TABLE]

Corollary 10 follows directly from Propositions 5 and 8 and states that, conditioned to the optimality of the output of Algorithm 1, steps 1), 2) and 3) described before Algorithm 2 are indeed optimal and preserve the minimum variance estimate. Yet, when a change in the location set ${\mathcal{I}}$ occurs, optimality is inevitably lost in the state reconstruction step 5) which is obviously only sub-optimal. In particular observe that the adaptive strategy we have proposed to handle new input locations is based on suitable initializations of $\widetilde{\mathbf{s}}_{{\mathcal{I}}(k)}(k)$ and $\widetilde{\Sigma}^{\mathbf{s}}_{{\mathcal{I}}(k)}(k)$ which lead to suboptimal filters. One could wonder if the optimal function estimate is inevitably lost or can be reobtained by the algorithm. Clearly, it is not possible to establish any kind of convergence if the input locations vary infinitely often in time. However, the following result shows that, after any change of the input location set, and in absence of other perturbations, fast (indeed exponential) convergence to the entire (infinite-dimensional) optimal estimate holds.

Proposition 11** (Asymptotic optimality of the sub-optimal estimator)**

Consider the system (11) obtained by sampling (9) over time instants $t_{k}$ satisfying $t_{k+1}-t_{k}>\Delta>0$. Assume that, up to the instant $t_{k}$, input locations have been added and/or removed from ${\mathcal{I}}$ as described above. Assume also that all the measurements collected for $t>t_{k}$ fall in the current sampling grid ${\mathcal{I}}(k)$, i.e., ${\mathcal{I}}(j)\equiv{\mathcal{I}}(k)$ for any $j>k$, namely, no other perturbation of the input location set occur. Let be $\widehat{f}(x,t):={\mathbb{E}}[f(x,t)|\mathcal{D}(t)]$ the optimal minimum variance estimate of $f$ given all the measurements up to $t$. Then,

[TABLE]

where the convergence is exponential in mean square sense.

8 Simulations

Here we present some simulations to show the effectiveness of the proposed estimation schemes.

We test both Algorithm 1 and Proposition 8 on synthetic and real field data. Finally, we test the proposed adaptive strategy on synthetic data only. All the simulations are run in MATLAB® on a 2,7 GHz Intel Core i5 processor with 16GB RAM.

8.1 Synthetic Data

In this section we compare, in terms of computational performance, the proposed estimation scheme with the classical iterative GP procedure (which we refer to as “truncated GP” since it is based on a finite memory approach [19]), which assumes perfect knowledge of the kernel for the modeled process. We recall [34] that the truncated GP procedure is characterized by a computational complexity of order

[TABLE]

being $q$ the finite memory buffer length. We work on a 1D space. More specifically, $\mathcal{X}$ consists of a line of length $100$ [p.u.] which has been uniformly sampled every $1$ [p.u.] ($|{\mathcal{I}}|=M=100$). The sampling time is fixed and equal to $0.2$ [s], while the simulation time is $10$[s]. We assume $M_{k}=M$, that is, we collect measurements from all the locations at every time instant and we do not perform predictions, i.e., $P=0$. Thus we collect a total of $10/0.2\times 100=5000$ measurements. Notice that, thanks to this choices, the computational complexities per iteration (see Section 6) reduce to $\mathcal{O}(rM^{3})$ for Kalman and to $\mathcal{O}(q^{3}M^{3})$ for the truncated GP, respectively. Therefore, $r$ and $q$ represent a measure for the complexity of the corresponding approach. Finally, measurements noise is set $\sigma=1$ [p.u.].

We test the proposed approach on a process whose kernel does not satisfy Assumption 3. In particular, the selected process is drawn from a spatio-temporal separable Gaussian kernel $K$ with

[TABLE]

where $\sigma_{s}=5$[p.u.] and $\sigma_{t}=\sqrt{2}$[s]. Observe that, in this case it is necessary to compute a rational approximation $\widehat{S_{\rm r}}(\omega)$ of the true PSD $S(\omega)$, in order to retrieve a state-space representation of the process according to Proposition 4. To do so we compute $\widehat{S_{\rm r}}(\omega)$ as the solution of the following parametric non-linear weighted least-squares problem

[TABLE]

where $r$ is the model order, while $\{a_{i}\}_{i=0}^{r}$ and $\{b_{i}\}_{i=0}^{r-1}$ are the coefficients of the spectral factor $W(\mathbf{i}\omega)$ of $S_{\rm r}(\omega)$.

We compare the estimation methods in terms of

1. (i)

CPU time per iteration; 2. (ii)

estimation fit computed as

[TABLE]

where $\widehat{\mathbf{f}}_{*}$ denotes the estimate obtained at the end of the simulation time, i.e., $T=10$[s], either using the proposed Kalman-based approach or the truncated GP; while $\widehat{\mathbf{f}}_{\rm np}$ denotes the classical GP estimate which uses all the available measurements, i.e., $q=\infty$.

For the truncated GP, Figure 6 shows the fit as a function of the buffer length $q$, while for the proposed Kalman-based Algorithm 1, the fit is plotted as function of the model order $r$. In general, it can be seen that, for the same level of complexity, i.e., $q$ vs. $r$, Algorithm 1 achieves a better fit. We stress the fact that the performance in terms of fit for the truncated GP highly depends on the ratio between the process and the measurements noise. Indeed, for high process noise, the information contained in the measurements collected during the last few iterations already contains all the necessary information to reconstruct the process. Thus, the fit curve would increase more rapidly. Conversely, Kalman is optimal hence it does not depend on the ratio. Table 1 reports the value of the Fit defined in(26) and of the CPU execution time for the proposed Kalman-based approach with $r=6$ against the truncated GP for three different values of buffer length $q$ corresponding to $\sigma$, $2\sigma$ and $3\sigma$ of the time kernel, respectively. Note that the proposed Kalman-based approach behaves almost perfectly as the classical GP approach using all the available measurements (reported in the table last row as $q=\infty$, which needs almost 20[s] to run). The only discrepancy is due to the rational approximation of the kernel needed to implement Algorithm 1. Conversely, the truncated GP needs a computational time of at least one order of magnitude higher to achieve the same level of estimation accuracy in terms of Fit. It is worth stressing that this depends on the time kernel used for estimation. Indeed, in the example above, since we used a Gaussian time kernel we needed a rational approximation $\widehat{S_{\rm r}}$ of at least $r=6$ to achieve 99.4% performance with Algorithm 1. Conversely, Table 2 reports the values obtained using the Laplace time kernel equal to

[TABLE]

which is characterized by a rational PSD, see Example 7. In this case, since the time kernel has a rational PSD, Algorithm 1 achieves 100% of accuracy. Moreover, it requires less CPU time than before since the state-space model corresponding to the Laplace kernel PSD is of order $r=1$. Conversely, to achieve a level of accuracy similar to one of Table 1, the truncated GP needs more memory steps and thus its CPU time keeps increasing.

8.2 Colorado Weather Data

As second application, we consider weather forecasting on real field collected data. We exploit the same data-set used in [37, 32, 23]. This consists of spatio-temporal weather data, i.e., precipitations and maximum and minimum temperature, collected every month during the years 1895-1997 from 367 different weather stations around Colorado, USA333https://www.image.ucar.edu/Data/US.monthly.met/CO.shtml. The data-set is actually a subset of a larger data-set including 11918 weather stations. In particular, the considered subset has been extracted from the larger one considering only stations laying in the rectangular lon/lat region $[-109.5,-101]\times[36.5,41.5]$ and deleting those collecting only one type of measurement. This leads to a data-set consisting of $453612$ measurements. A great amount of data are Not Available (NA). This makes it suitable for prediction and forecasting. Before presenting our simulations in details, a first comparison between the proposed Kalman-based approach and the truncated GP method is offered in Table 3. This reports the memory and the computational requirements (per iteration) for Algorithm 1 and for both the classical and the truncated GP methods. For the latter, the table shows the requirements assuming to use the data from three time windows of different length. First of all it is worth noticing how the classical approach is not feasible. Conversely, even if the truncated approach leads to more reasonable implementations, it is still not comparable to the Kalman-based approach. Observe that the time window length (and thus memory and computational requirements) for the truncated GP method largely depends on the process and, in particular, on its temporal correlation. A common choice is to consider data within the $3\sigma$ confidence interval444With $3\sigma$ confidence interval we mean a time window $T$ such that $\int_{0}^{T}h(\tau)d\tau=0.99h(0)$. Since for the Normal distribution this translates to consider 3 deviations from the mean we use the same nomenclature. only. In our specific case, given the estimated values for the kernel hyper-parameters (see below), this translates in using data from the last $\approx 2$ years.

In the following we focus on precipitations measurements only, assuming a noise standard deviation equal to $5\%$ of the corresponding absolute measured value. To model the spatial covariance we use an exponential kernel

[TABLE]

while, in order to exploit the seasonal periodicity of the precipitations ($f=1/12$), for the time covariance we resort to a stationary, periodically decaying kernel equal to

[TABLE]

which is characterized by a rational PSD equal to

[TABLE]

which leads to a factorization (6) with

[TABLE]

The kernel hyper-parameters are estimated by minimization of the -log-marginal likelihood (see Appendix 11) over all available data in the period $1895\div 1995$. For inference, from the remaining two years data, i.e., $1996\div 1997$, we extract a subset corresponding to $80\%$ of randomly picked weather stations. For test we use the data collected from the remaining $20\%$ of the stations. Figure 7 shows the contour of the estimates (Figure 7a) and of the corresponding posterior variance (Figure 7b) in geographic coordinates over the entire region of interest during October 1997. Black crosses (+) identify weather stations used for inference while black circles ($\bullet$) identify stations used for test. Interestingly, from Figure 7b, observe how the minima of the posterior variance are usually attained in the locations corresponding to the weather stations used for inference. This is not always the case since there might be weather stations where measurements are missing.

Finally, Figure 8 shows the time evolution corresponding to two randomly picked sensors, namely, sensor $\#18$ (Figure 8a) which owns to the set of sensors used for inference, and sensor $\#130$ (Figure 8b) which owns to the test set. Figure 7 explicitly shows their geographical positions marked with red circled crosses. For sensor $\#18$, notice that the estimates are always within the confidence interval and how the procedure automatically extends the result even when actual measurements are not available. Similarly, for sensor $\#130$, whose measurements are not used for inference, both estimate and confidence interval are quite accurate with respect to the corresponding measurements, showing the effectiveness of the proposed procedure.

8.3 Adaptive Input Location Set

Before presenting some concluding remarks, we test here the adaptive grid strategy proposed in Section 7 on a synthetic data-set. For simplicity we work on a 1-D environment, i.e., a line of length $1$[p.u.]. We want to estimate a spatio-temporal GP process with kernel

[TABLE]

where $\sigma_{s}=0.005$[p.u.], $\sigma_{t}=100$[s] and $\lambda=1$. We assume to have at our disposal a mobile robot which is able to collect one single measurements per iteration every $T=1$[s] from $50$ uniformly distributed points along the 1[p.u.] line. More precisely, at each iteration the robot moves along the line according to a random walk with a forcing term in order to avoid it to jump back and forth around the same measurement locations. The motivating application we have in mind is monitoring and estimation of meteorological phenomena from punctual sampling of, e.g., cumulus-type clouds [35, 36] as pursued in the ongoing project “SkyScanner”555https://www.laas.fr/projects/skyscanner/ which comprises different partners among which the RIS research group at LAAS/CNRS in Toulouse, France. Since clouds move in space and time it is necessary to follow them while extracting useful information from the collected spatio-temporal measurements. In such a scenario, the proposed strategy perfectly suits since the estimator, while following the moving clouds, updates the grid locations retaining only the last and more valuable ones.

We compare the adaptive strategy against the Kalman-based Algorithm 1. We recall that Algorithm 1 assumes ${\mathcal{I}}$ to be fixed and equal to the entire set of 50 measurement locations. Thus it produces an estimate for every single measurement location for the entire simulation time, chosen to be equal to $100$[s]. Conversely, for the adaptive strategy we choose a maximum memory capacity of $10$ input locations. Thus, the location set as well as the state space are expanded until $10$ different locations are visited. After that, when a new location is visited, the oldest visited one is discarded. In addition, to show the optimal limit behavior of Proposition 11, for the adaptive strategy, after $50$[s] we block the input locations set to coincide with the last $10$ visited input locations.

Figure 9 shows two different time shots corresponding to two different time instants. In particular, Figure 9a corresponds to $t=20$[s] when the memory capacity must still be hit. Interestingly, the estimate returned by the proposed adaptive strategy nicely reproduce that returned by Algorithm 1. Conversely, Figure 9b corresponds to $t=90$[s] when the memory capacity have been hit, the algorithm have already started to drop old locations and the input location set have been blocked. In this case, observe how, over the common set of locations, thank to the asymptotic optimality results of Proposition 11, the statistics returned by the adaptive strategy perfectly coincide with those of Algorithm 1. Both the figures report the corresponding confidence interval computed simply as $\pm 3\sqrt{\mathrm{diag}(\Sigma^{\mathbf{f}})}$. It is worth noting that for Algorithm 1 $\Sigma^{\mathbf{f}}$ actually correspond to the filter performance. Conversely,because of the suboptimality, for the adaptive strategy $\widetilde{\Sigma}^{\mathbf{f}}$ is not the true error covariance matrix (the true error covariance is retrieved in the limit when the input location set is blocked) but just an approximation.

Finally, to better appreciate the performance of the proposed algorithm, we refer the interested reader to a full video of the simulation666http://automatica.dei.unipd.it/people/todescato/publications.html.

9 Conclusions & Future Works

In this work we focused on building an efficient GP estimator for spatio-temporal dynamical Gaussian processes. The main idea was to couple Kalman-filtering and GP regression. In particular, assuming space/time separability of the covariance (kernel) of the process and rational time spectrum, we built a finite-dimensional discrete-time state-space process representation over a finite dimensional set of input locations. Our major finding is that the Kalman filter state at instant $t_{k}$ represents a sufficient statistic to compute the minimum variance estimate of the process at any $t\geq t_{k}$ over the entire infinite dimensional spatial domain. This result can be interpreted as a novel Kalman representer theorem for dynamic GPs. Then, we extended the study to situations where the sampling locations can vary over time by designing a novel computational scheme. The proposed strategies are tested on both synthetic and real field data, also providing comparisons with standard GP and truncated GP estimation techniques. Future work will consider application of the proposed strategies in real-time applications as cloud monitoring and the extension to distributed GP regression for multi-agent systems.

10 Proofs of Propositions

10.1 Proof of Proposition 4

First of all, notice that the process $\mathbf{f}(t)$ is a Gaussian process since it is the solution of a linear differential equation driven by Gaussian noise $\mathbf{w}(t)$. To conclude the proof we need to show that the covariance of $\mathbf{f}(t)$ is indeed $\bar{K}=\bar{K}_{\rm s}h(\tau)$. As previously shown, the first two equations of model (9) are the state-space representation of the rational power spectral density $S_{\rm r}(\omega)$ thus ${\mathbb{E}}\left[z_{i}(t+\tau)z_{i}(t)\right]=h(\tau)$. It follows that

[TABLE]

10.2 Proof of Proposition 5

First, consider time instants $t=t_{k}$, $k\in{\mathbb{Z}}_{+}$. In this case the result directly follows by applying the standard Kalman Eqs.(5) to the discrete model (11). Finally, for the case when $t\in]t_{k},t_{k+1}[$ note that, given the last state-space estimate $\widehat{\mathbf{s}}(k|k)$, Algorithm 1 returns the best “open-loop” time prediction in accordance with the underlying state-space model.

10.3 Proof of Proposition 8

Since $K$ is separable, we have that $K(x,x^{\prime},t,t^{\prime})=K_{\rm s}(x,x^{\prime})h(\tau)$ with $\tau:=t-t^{\prime}$. We assume to be at time instant $t\geq t_{k}$, while $t_{j}<t_{k}$ represents a generic previous sampling instant, i.e., when measurements have been collected. Now, let us define the following additional symbols

[TABLE]

thus implying that $\Gamma_{x}=\Gamma_{x}(0)$ and $V_{\mathbf{f}}=V_{\mathbf{f}}(0)$.

We first study $p(\boldsymbol{\varphi}(t_{k}),f(x,t)|\mathbf{f}(t))$. For the conditional variance, we have that

[TABLE]

where

[TABLE]

and where

[TABLE]

From (29) it follows that the second term in the right hand side of (28) is equal to

[TABLE]

Hence, by subtracting it to the first term on the right hand side of (28), i.e., $Var([\boldsymbol{\varphi}^{T}(t_{k})\ f(x,t)]^{T})$, the last column and the last row cancel out (except for the diagonal block), meaning that $\boldsymbol{\varphi}(t_{k})$ and $f(x,t)$ are conditionally independent given $\mathbf{f}(t)$. Thus, we have that

[TABLE]

Thank to this we can write

[TABLE]

so $f(x,t)$ is conditionally independent from all the past contained in $\boldsymbol{\varphi}(t_{k})$. Then, we have that

[TABLE]

where the first equality holds because we are conditioning on a larger $\sigma$-algebra; the second holds thanks to (30) and the third comes from Proposition 5. Finally, for the posterior variance we exploit the result of Lemma 1 contained in Appendix A of [33] combined with the conditional independence stated in Eq. (30). For convenience, we recall

Lemma 12** (Lemma 1 in [33])**

Assume that

[TABLE]

where $z^{*}$ is a scalar and $\epsilon$ is independent of $[z^{*}\ \boldsymbol{\eta}^{T}]^{T}$. Then,

[TABLE]

Following the above notation we have that $z^{*}=f(x,t)$, $\sigma_{*}^{2}=V_{x}$, $\boldsymbol{\eta}=[\boldsymbol{\varphi}(t_{k})^{T}\ \mathbf{f}(t)^{T}]^{T}$, $\Gamma=\Gamma_{x}$, $V=V_{\mathbf{f}}$, $\Sigma_{\varepsilon}=R$ and $F=I$. However, in view of Eq. (30), i.e., the conditional independence of $f(x,t)$ from $\boldsymbol{\varphi}(t_{k})$, in our specific case, given $\mathbf{f}(t)$, the past process contained in $\boldsymbol{\varphi}(t_{k})$ does not bring any useful additional information and thus, it is sufficient to consider $\boldsymbol{\eta}=\mathbf{f}(t)$. Finally, it is easy to see that in our specific case $Var[\boldsymbol{\eta}|\mathbf{y}]\equiv\Sigma^{\mathbf{f}}(t)$ as returned by Algorithm 1 and thus the result follows.

10.3.1 Proof of Remark 9 – On non-stationary time kernels

The proof of Remark 9 follows almost straightforward from the above result. It is indeed sufficient, slightly changing the notation, to define

[TABLE]

instead of $\Gamma_{x}(\tau)$ and $V_{\mathbf{f}}(\tau)$, respectively. Then, the same reasoning applies and it is possible to conclude for the conditional independence of $f(x,t)$ from the past $\boldsymbol{\varphi}(t_{k})$ given $\mathbf{f}(t)$. The interesting and delicate step here is to conclude that

[TABLE]

still holds. This is indeed true because, conditioned to an exact state space representation of the kernel $h(\cdot,\cdot)$, which now does not satisfy Assumption 3, Algorithm 1 still is an exact procedure which outputs the minimum variance estimate $\widehat{\mathbf{f}}(t)$ of $\mathbf{f}(t)$.

10.4 Proof of Proposition 11

Consider system (11) and two filters associated to it which differ only in their initializations at $k$ obtained assigning two distinct couples $\widehat{\mathbf{s}}(k|k)$ and $\Sigma^{\mathbf{s}}(k|k)$777In our context, one can think e.g. of the first filter initialized through (24) and of the second one which has used all the past measurements to obtain the correct initial condition. The time-varying system (11) is obtained by sampling the continuous-time system (10) and, hence, is exponentially stable. In fact, the transition matrix $F$ in (10) is stable by construction so that there exist positive scalars $a$ and $b$, with $b<1$, such that

[TABLE]

So, for any integer $k$ and $j$, the transition matrices $A(k)$ in (11) must satify

[TABLE]

where $c=b^{\Delta}$ (recall that $t_{k+1}-t_{k}>\Delta$ by assumption) so that $c<1$.

Now, in view of this result, according to the definitions in Section 2 of [38] it is easy to conclude that (11) is uniformly stabilizable and detectable. Hence, Theorem 5.3 in [38] ensures that the two filters are exponentially stable and, using the same arguments contained in the proof of Theorem 7.5 of [39], that they converge to the same estimator. A common Kalman filter state will be thus asymptotically reached (in the mean square sense). As $t\rightarrow+\infty$, mean square convergence of $\widetilde{f}(x,t)$ to the optimal estimator over the entire spatial domain now comes directly from Proposition 8 (Kalman representer theorem).

11 Hyper-parameter estimation

Given a certain data-set, in order to maximize the estimation performance of the filter, it is necessary to perform a suitable selection of the kernel hyper-parameters. To this end, a natural choice is the minimization of the marginal negative log-likelihood associated to a specific set of “training” measurements. To do so, given a set of possible hyper-parameters combinations, it is possible to run in parallel multiple Kalman regressors, each one associated to a particular combination of parameters. From the Kalman filter Eqs. (5), the marginal negative log-likelihood associated to each regressor, and denoted with the symbol $\ell(k)$, can be recursively updated, at iteration $k$, as

[TABLE]

where $i(k)$ and $I(k)$ represent the innovation vector and its covariance during iteration $k$, respectively, and it is assumed $\ell(0)=0$.