{'fraction_non_alphanumeric': 0.03348393091736358, 'fraction_numerical': 0.007850299593066103, 'mean_word_length': 4.561297220242338, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical Gaussian process Kalman filtering for spatiotemporal systems

We present a novel Kalman filter for spatiotemporal systems called the numerical Gaussian process Kalman filter (GPKF). Numerical Gaussian processes have recently been introduced as a physics informed machine learning method for simulating time-dependent partial differential equations without the need for spatial discretization. We bring numerical GPs into probabilistic state space form. This model is linear and its states are Gaussian distributed. These properties enable us to embed the numerical GP state space model into the recursive Kalman filter algorithm. We showcase the method using two case studies.


I. INTRODUCTION
State estimators are algorithms for reconstructing a system's state from a stream of noisy online measurements and model predictions. The model is usually in the form of a probabilistic state space model. Given linear finite-dimensional dynamics, the process and measurement equation read as with the state at time t as x t ∈ R dx , and the measurements as y t ∈ R dy . Dynamics are described by matrix A while measurements are obtained through C. Process q t−1 and measurement noise r t are modeled as Gaussian distributed zero-mean white noise with covariance matrix Q and R, respectively. Process and measurement equations (1) - (2) can also be represented by probability density functions p(·) that take the form of Gaussian distributions N(·) due to the nature of the noise processes. In particular we have p(y t |x t ) = N (y t |Cx t , R) .
The following question arises: How should the measurement y t be used to correct the model prediction of x t in real time?
One answer to this question is the Kalman filter. In fact, it gives the optimal estimate m opt t of the current state x t under consideration of the measurement history up to the present y 1:t = {y 1 , . . . , y t } [1]. It is optimal in the sense that the expected squared error between estimated and true state conditioned on the measurement history is minimized. The Kalman filter estimates the state in a recursive manner. In Bayesian statistical terms this can be written as p(x t |y 1:t ) = p(y t |x t )p(x t |y 1:t−1 ) p(y t |y 1:t−1 ) .
Since the recursion starts with a Gaussian distribution and the dynamics are linear, all the above distributions remain Gaussian, and the posterior estimate can be obtained in closed form by calculating its mean and covariance.
Remember that the above model equations described finitedimensional systems. Infinite-dimensional filtering theory certainly exists [2], but direct practical infinite-dimensional Kalman filtering does not to the best knowledge of the authors. Finite-dimensional approximations are therefore needed. Using spatial discretization methods, e.g. the finite volume method [3], one usually obtains a large system of ordinary differential equations which can be formulated as the above state space model (1) - (2).
In this paper, we will take a different route. We will use socalled numerical Gaussian processes [4] to solve spatiotemporal models. Gaussian processes (GPs) are a probabilistic machine learning method that can be used for regression [6]. Numerical GPs are structured by the spatiotemporal model in the form of a (linear) time-dependent partial differential equation where L x : is a linear (integro-)differential operator determining the dynamics, and B x is a linear functional imposing the boundary condition. Due to the probabilistic nature of GPs, they provide us with uncertainty quantification. This will be of much benefit for their use in state estimation.
The contributions of this paper are as follows: • we extend numerical GPs with an output channel for online measurements and formulate this as a probabilistic linear state space model which is by definition Gaussian distributed, see Section III; • this allows us to put numerical GPs through the Kalman filter algorithm and therefore derive the numerical Gaussian process Kalman filter (GPKF), see Section IV; • using numerical case studies of spatiotemporal systems, we showcase the numerical GPKF, see Section V. spatial discretization temporal discretization Fig. 1: Spatial discretization methods lay down a mesh on the spatial domain and describe how neighboring mesh cells should interact with each other dictated by the PDE. This results in a system of ODEs that is simulated with traditional solvers. Numerical Gaussian processes [4] on the other hand start with a temporal discretization. The resulting time discretized PDE is used to build up a structured Gaussian processes that is evaluated on arbitrarily chosen regression points. To simulate, one recursively formulates the so-called posterior Gaussian process. Figure adapted from [5].
Preliminary results were presented in the conference paper [7]. This article significantly extends this previous work. In contrast to [7], the preliminaries now contain a didactical example to streamline the section. Furthermore, numerical GPs are explicitly expressed as probabilistic state space models in this paper. Moreover, boundary conditions are treated as part of the model and not as measurements. This affects the KF algorithm as such that boundary conditions enter through the prediction step, rather than the update step of the Kalman filter. Finally, more sophisticated case studies are presented.
Related works: Infinite-dimensional Kalman filtering with Gaussian processes has previously been done by [8], [9]. Therein, certain covariance functions are shown to be convertible via Fourier-transforms to infinite-dimensional state space models. Using these models, infinite-dimensional Kalman filtering can then be done. This approach does not need a structured PDE to inform the GPs, contrary to numerical GPs. Of course, this also means that the GP has to learn the dynamics from the full state vector which might not be accessible in the first place.
Another fully non-parametric approach is given in [10]. Assuming space-time separability of the covariance function, a finite-dimensional discrete-time state space model is built up that can be used for Kalman filtering. Again, the full states need to be accessible for this method to work.
Finite-dimensional filtering with Gaussian processes has been done by [11]. Different Bayesian filters such as particle filters, extended and unscented Kalman filters are recovered. Closely related is the work of [12], wherein an analytic moment-based Gaussian process filter is presented. A general perspective on finite-dimensional Gaussian filtering is given in [13].
Mathematical notation: Scalars are lowercase non-bold symbols x, while vectors are bold x, and matrices are uppercase and bold X.
A random variable x stemming from probability density function p(·) is symbolically written as x ∼ p(x). Gaussian random variables and their distributions are denoted as x ∼ N(x|m, P ) with mean m and variance P .
A covariance matrix K (X, X) ∈ R M ×M in the Gaussian process regression context is built up by evaluating the underlying covariance function element-wise

II. PRELIMINARIES
In this section we will introduce Gaussian process regression, linear operators in combination with GP regression, and numerical GPs for time-dependent partial differential equations.

Gaussian process regression
A Gaussian process is a stochastic process f (x) that is fully defined by its mean function m(x) and covariance function k(x, x ) What makes GPs usable in practice is their property that any finite dimensional collection of random variables with In regression, we wish to infer the process f (·) based on its inputs x and its (possibly noisy) output y. The output is not restricted to algebraic expressions such as y = f (x) + noise, but can also be generalized to linear transforms of the process y = L x f (x) + noise.
We will illustrate GP regression using the following differential example   with a Dirichlet boundary condition rendering the boundary operator to the identity operator B x = I x , and therefore f b = f (x b = 0). To start, we place a GP prior on f (x) and write where we chose the prior consisting of a zero-mean function and a squared exponential covariance function k(x, x ) = Here, σ 2 and l are hyper-parameters that permit to adjust the shape of the covariance function, and therefore the behavior of f (x), to the data.
Next, we make use of the fact that a linear transform of a GP is another GP that is structured by the linear transform, see [14] and the references therein. We therefore have For our example, k f f and k bb are identical due to the Dirichlet boundary condition. The covariance function k yy can easily be derived under use of the formal definition of a covariance function (9). Note the important distinction between L x and L x , meaning that the operator should be applied with respect to either the first or second argument, that is x or x . We continue by formulating the evaluated joint model p(f , f b , y) at test locations x * = {x * 1 , . . . , x * S }, the boundary x b , and observation locations with covariance matrix Before doing inference, we learn the hyper-parameters of the covariance function by minimizing the negative log marginal likelihood (NLML) − log p(f b , y). We will touch on this more thoroughly later on.
Inference is now done by calculating the conditional distribution p(f |f b , y). For Gaussian distributions, this results in another Gaussian with conditional mean Here, another hyper-parameter σ 2 n has been introduced to model measurement noise as additive zero mean Gaussian white noise. This hyper-parameter is also learned via the NLML. As can be seen in Fig. 2, the structured GP (18) is able to infer f from derivative observations.
To summarize, we give a cooking recipe for probabilistic inference with GPs: 1) Carefully place a GP prior to prevent inversion of linear operators (prior on f in our example to prevent inversion of d dx ). 2) Build up output GPs that are linear transforms of the GP prior (y, f b in our example). 3) Construct the joint model p(f, f b , y) and learn the hyperparameters via the marginal distribution p(y, f b ). 4) Do inference by conditioning p(f |f b , y).
Numerical Gaussian processes for time-dependent partial differential equations So far we looked at static examples. We are however interested in spatiotemporal systems described through timedependent partial differential equations (PDEs). We assume that the PDE model is known and linear. Given an initial condition and boundary conditions, we wish to solve the PDE via GPs. We can do this with the recently introduced numerical GPs [4]. We will present the essential idea here, referring the reader to the original paper for further details.
Numerical GPs are built on the fact that linear transforms of a GP result in another GP. In fact, numerical GPs use a very similar procedure as shown earlier.
The first step is to discretize the time-dependent PDE (7) in time using one of the many existing methods (explicit Euler for simplicity here) Next, a GP prior is placed on and therefore f t is also a GP defined as and so is the boundary condition Then, the evaluated joint model p(f t , f b t , f t−1 ) on arbitrary test points X t , the boundaries X b , and previous points X t−1 is formulated as Here, the inputs to the individual covariance matrices are omitted. They can be inferred from the sub-and superscripts, e.g.
Here the solution from the previous simulation step f t−1 acts as an artificial measurement. Hyper-parameters are optimized in each simulation step by minimizing − log p(f b t , f t−1 ).

III. NUMERICAL GAUSSIAN PROCESSES AS PROBABILISTIC STATE SPACE MODELS
In this section we will derive a probabilistic state space model from numerical Gaussian processes. This probabilistic state space model is linear and its random states are, by definition of GPs, Gaussian distributed. These model properties allow us in Section IV to go through the Kalman filter algorithm, thereby deriving the numerical Gaussian process Kalman filter.
We start by introducing an output channel for noisy online measurements that are linear transforms of f (t, x). Mathematically this reads as The measurement operator is H x : Obtained measurements are point evaluations of f y t at locations y i . In this work, we consider measurement noise r t as zeromean spatiotemporal white noise with constant variance over y. To model this, we place a GP prior on r t that is independent of all the other GP priors. We write where the covariance function is k rr t,t = σ 2 r δ(y − y ) and E[r t r t+τ ] = 0. Other additive noise models are conceivable. Particularly the spatial white noise assumption could be relaxed to something more suitable under use of an appropriate covariance function [6].
Regarding the PDE itself, we add spatiotemporal process noise to the PDE to obtain Discretization in time yields with process noise modeled as an independent zero-mean GP prior with a white noise kernel . Of course, the above discretization is mathematically not rigorous because (30) is not continuous due to the spatiotemporal white noise. A more rigorous treatment would require Itô-Calculus [15]. We could formally circumvent this by first discretizing and then adding process noise q t−1 (x) (without ∆t), as done in the finite-dimensional state space models. However, we found that exclusion of ∆t in the process noise term drastically reduced estimation performance. Moreover, above we recovered a time-discrete form that is similar to the Euler-Maruyama method. This brief discussion lays bare the gap between proper theoretical treatment and actual implementation, as lamented as early as in [16].
The updated GP for f t (x) is The complete multi output GP now reads as (33) All covariance functions inside K GP are known as they can be built up from the prior covariance function k f f t−1,t−1 (for the explicit Euler) and the respective linear operators A x , H x or functionals B x , see Appendix C.
We can recover a probabilistic state space model from (33). For this purpose we assume that (33) has been evaluated at measurement locations Y , boundaries X b and arbitrarily chosen test locations X. The covariance functions therefore become covariance matrices, e.g.
The process equation is with the covariance matrix taking the role of the dynamic matrix in a state space model. Here, A t is time-varying, indicated by the subscript, due to changing hyper-parameter values and possibly changing regression points. The process noise covariance matrix of the state space model, usually denoted Q, is given by The measurement equation reads as with taking the role of the measurement matrix in a state space model. The measurement noise covariance matrix of the state space model, usually denoted R, is given by The noise covariance matrices Q and R in a classical state space model are usually static matrices. The GP model however naturally allows for an adaptive noise description, due to online adaptation of the hyper-parameters, thereby drastically reducing the usually required tedious fine-tuning of the Kalman filter. This will be showcased with case studies in Section V and discussed in Section VI.
For the recursive Kalman filter algorithm we will also need the following with Equipped with a probabilistic state space model, we are now ready to write down the recursive Kalman filter equations.

IV. NUMERICAL GAUSSIAN PROCESS KALMAN FILTER DERIVATION
In this section we will embed the numerical GP state space model into the Kalman filter algorithm. Afterwards, we will briefly discuss the Kalman filter marginal likelihood for the use of hyper-parameter estimation.
There are two ways to derive the Kalman filter. The first derivation is rooted in recursive least squares regression, see e.g. [17]. The second derivation takes a Bayesian filtering perspective, see e.g. [18]. In this article, we take the latter perspective to derive the numerical Gaussian process Kalman filter since the Bayesian perspective neatly fits into the Gaussian process framework. Recall that in Kalman filtering, we want to calculate the probability distribution of a dynamic state given all the measurement history up to the present. For spatiotemporal systems, we additionally condition on the boundary data history. Boundary conditions are treated as part of the model and therefore enter through the prediction step of the Kalman filter. According to Bayes' rule, the posterior distribution of the state estimate is given by We make the usual assumptions for the Kalman filter [18] and an additional one for the boundary data: Assumption 1) States are assumed to be Markovian, i.e. the current state f t is conditionally independent of anything that happened before t − 1 given the previous state f t−1 . Assumption 2) Given the current state f t , the current measurement f y t is conditionally independent of the measurement f y 1:t−1 and state histories f 1:t−1 , as well as the boundary history up to the present f b 1:t . Assumption 3) Given the previous state f t−1 , the current boundary condition f b t is conditionally independent of the measurement f y 1:t−1 and boundary history data f b 1:t−1 . In this article, only one-step discretization schemes such as the Euler method are considered. Using multi-step methods would require to extend the Markovian property to an accordingly higher order.
The individual terms of (43) can be calculated in closed form for Gaussian distributions as The road map ahead to derive these expressions is shown in Fig. 3. To calculate the prior distribution p(f t |f y 1:t−1 , f b 1:t ), we will first formulate the joint distribution between states f t , f t−1 , and current boundary condition f b t , conditioned on measurement and boundary histories. Under use of Assumptions 1 and 3 1 , the prior joint is calculated as The covariance matrices A t and P GP,f t have been introduced in (34), while A f b t and P GP,f b t have been introduced in (42). Using Lemma A two times, the joint mean m of (47) is and the covariance is Here we introduced The prior distribution is now obtained from (47) by first marginalizing over f t−1 and then conditioning on f b t . Prior mean m − t and covariance Equipped with this, we can calculate the joint distribution of the current state (prior distribution) and measurements (likelihood). This will allow us to write down the posterior distribution later on. Using Assumption 2, we have The covariance matrices C t and P GP,f y t were introduced in (37). We use Lemma A one more time to get the joint mean as and the covariance as (57) The posterior distribution is obtained by conditioning the joint distribution (55) on the current measurement using Lemma B. The posterior mean is and the corresponding posterior variance is To summarize the numerical GPKF, we have the prediction step as and the update step as Due to having to include boundary conditions, the prediction step doesn't perfectly align with the traditional KF. It is not possible to simply write the prior mean prediction as the previous posterior mean mapped through the dynamic matrix A t . Treating boundary conditions as measurements results in much simpler expressions that smoothly align with the traditional KF, see [7]. Of course this means that prediction steps do not abide boundary conditions which can be problematic when multiple prediction steps have to be computed between update steps.

Hyper-parameter estimation via the marginal likelihood
Hyper-parameters are estimated in each update step. In GP regression the hyper-parameters are usually estimated by minimizing the negative log marginal likelihood (NLML). The marginal likelihood of the Kalman filter, marginal with respect to the current state f t , is found in the recursive Bayes denominator (43) as The corresponding NLML is with N y being the number of measurements. The partial derivatives of this NLML required for minimization can be calculated analytically and are shown in Appendix D.

V. SIMULATION CASE STUDIES
In this section we will showcase the numerical Gaussian process Kalman filter (GPKF) using different case studies. For all case studies the relative error between posterior mean estimate m t and analytical solution f t,ref evaluated on the same points X is defined as with euclidean distance · . All case studies have been implemented in MATLAB. A.
Step shaped one-dimensional advection equation We study the one-dimensional advection equation with initial and boundary condition as Here, H(x) is the Heaviside step function. We solve (72) with the implicit Euler and a step size of ∆t = 0.005. The GP prior is f t ∼ GP(0, k f f t,t (x, x )) with a neural network covariance function .

(75)
This covariance function is capable of capturing discontinuities. The hyper-parameters are σ 2 0 and σ 2 . Measurements are noisy point evaluations so that the measurement equation reads as with r t ∼ N(0, σ 2 r,true ) being white noise. Measurements are received every third model step size so that between each KF update step there are three KF prediction steps. Hyperparameters are learned by minimizing the NLML (70) in each update step.
Estimation results are illustrated in Fig. 4. The initial estimate (blue curve with uncertainty band in Fig. 4a) is far off from the true solution (black line). However, rather quickly the estimate converges closer, as indicated by a dropping relative error (black curve in Fig. 4b). The relative error drops to a value of around 0.1 before it slowly starts to increase. This divergence is not uncommon in Kalman filtering [19]. One possible cause, that is also present in Fig. 4a, is a low posterior covariance matrix P t . As a result, the KF puts less importance on measurements and more on model predictions, which in this case aren't perfect due to the temporal discretization required for numerical GPs. Remedies such as the Kalman filter with fading memory exist [17].
The measurement noise level σ r is correctly learned online to be σ r,true . Estimated process noise σ q is non-zero, likely reflecting the aforementioned model imperfection. For a few time steps at the beginning its estimation gets exceedingly large. This is cut off from the plot.

B. Liouville's equation
We study Liouville's equation wherein the vector fieldẋ is governed bẏ Initial and boundary conditions are   We discretize (78) in time using the implicit Euler method with a step size of ∆t = 0.005, and place a GP prior on f t ∼ GP(0, k f f t,t (x, x )) with the following squared exponential covariance function Here, the hyper-parameters are σ 2 n , l 1 , and l 2 . Through noisy measurements we receive the marginal distribution of x 1 as with r t ∼ N(0, σ 2 r,true ) being white noise. The vector fieldẋ describes harmonic oscillations, plotting the solution f (t, x) results therefore in circular motions. The simulation runs for one full rotation T sim = 2π. Although measurements are formally the x 1 -marginal, over time measurements will be marginals of the initial distribution from all angles due to the circular motion of the joint distribution f (t, x). Interested readers are referred to [20] for further details. As discovered and explained by [21], this case study illustrates the relatedness to tomography. Measurements are received every Tsim 64 , which translates to 19 KF predictions steps for every KF update step.
State estimation results are displayed in Fig. 5 while noise level estimations are displayed in Fig. 6. Through measurements of the x 1 −marginal f (t, x 1 ) the numerical GPKF is able to reconstruct the latent x 2 −marginal distribution f (t, x 2 ) and the joint distribution f (t, x 1 , x 2 ). Noise level estimation σ r oscillates around its true level σ r,true , while process noise σ q is set to zero by the numerical GPKF. The relative estimation error drops significantly before settling in at around 0.16.
Hyper-parameter optimization requires an ad-hoc approach for this case study. Simply using the marginal likelihood p(f y t |f y 1:t−1 , f b 1:t ) of the Kalman filter as before will result in ever larger hyper-parameter values of the marginalized dimension, i.e. l 2 . We conjecture that due to the marginalization of x 2 , its impact on the marginal likelihood has been largely nullified, rendering l 2 unidentifiable. A remedy to this would be to use the Radon-transform and thereby treat this case study as a problem of tomographic type. Through a change of variables, known as input warping in the GP community, the Radontransform would convolute the dimensions x 1 and x 2 , and therefore also both length-scales l 1 and l 2 , thereby preventing cancellation of l 2 through the measurement marginalization. See [22] for the use of GPs in tomographic inverse problems. Instead, we successfully employed the function as the objective function to minimize. Here, we used the previous posterior mean m t−1 as fixed targets for f t−1 . The number of data points is N . This objective function now contains both length-scales l 1 and l 2 explicitly and therefore works.

VI. DISCUSSION
In this section we elaborate on computational load, optimality, observability, and distinguishing features of the numerical Gaussian process Kalman filter. Right plots also show noisy measurements of f (t, x 1 ) as blue circles. Middle plots show the posterior mean estimate for the joint distribution f (t, x 1 , x 2 ), as well as the analytical solution (transparent). Uncertainty bands have been omitted. The initial estimate is based on the black dots (first row, middle plot).
Online measurements of f (t, x 1 ) shown as blue dots (right plots).

On computational cost
The main computational cost lies in the hyper-parameter optimization which requires the inversion of the measurement covariance matrix. This is inherent to GPs and scales cubically with the number of training points. There are however different approximations to reduce this to less than cubic scaling, see e.g. [23], [24].
Actual simulation of the model with numerical GPs, i.e. propagation through time, does not require numerical integration and instead only involves algebraic manipulations as given by (21). The two dimensional simulation ran for 111 min, while the one dimensional simulation ran for 9.5 min on a 2020 13" MacBook Pro 2, 3 GHz quad-core i7.
In addition to the number of dimensions x ∈ R dx , computational load is also influenced by the prior covariance function and the PDE structure. Both can lead to less or more difficult to handle contours of the NLML, thereby influencing the computational load via the hyper-parameter optimization. For example, we found simulations of the same PDE using a squared exponential covariance function to run faster and more robustly than using a neural network covariance function.

On optimality
For a linear and Gaussian system, the Kalman filter gives the optimal estimate of the state m opt t given the measurement history n y 1:t . For spatiotemporal systems, we also condition on the boundary history n b 1:t . Optimal is defined in this context as m opt t (n y 1:t , n b 1:t , θ) = arg min mt E (n t − m t ) T (n t − m t )|n y 1:t , n b 1:t , θ . (86) Here, the estimated state is m t and the true (stochastic) state is n t . The hyper-parameters of the Gaussian process are denoted as θ. They influence the state estimate through i) the covariance function itself and ii) their estimated values. Although working with spatiotemporal systems, the input space of numerical GPs is solely the spatial domain. Choosing an appropriate covariance function should therefore be based upon the solution shape with respect to spatial dimensions. For the step shaped one-dimensional advection equation for example, we chose the neural network covariance function due to its ability to handle discontinuities. In the original numerical GP work [4] the Burgers equation is solved with the neural network covariance function as well because its solution develops a shock front.
As showcased and argued in [6], multiple optima of the marginal likelihood can certainly exist, but are usually not a problem. Indeed, the different optima represent different interpretations of the data.

On observability
We showed that numerical Gaussian processes can be brought into probabilistic state space model form. By design of numerical GPs, this state space model is structured by the underlying time-dependent partial differential equation. Showing observability of the PDE might therefore be transferable to the numerical GP state space model (34), (37).
Closely connected to observability is the concept of identifiability. This relates to the estimation of (hyper-)parameters. As seen and explained in the second case study involving Liouville's equation, an ad-hoc replacement for the marginal likelihood might be necessary to make hyper-parameter estimation work.

On distinguishing features
Kalman filters require fine-tuning before they deliver acceptable performance. This fine-tuning involves adjusting the process Q and measurement noise covariances R, as well as the initial error covariance matrix P 0 .
Fine-tuning of numerical GPKFs is either much less tedious or not required at all. Once initialized with reasonable values, process and measurement noise level hyper-parameters σ q and σ r are adapted online to the measurement stream.
As for the initial error covariance matrix P 0 , no manual fine-tuning was required for the case studies due to the probabilistic nature of GPs. P 0 was simply obtained by calculating the posterior GP variance (21) for given initial data. This was done without optimizing the hyper-parameters, as can be seen by the bad initial fit to the step function in Fig. 4a.

VII. CONCLUSIONS
In this article, we presented a new type of Kalman filter for spatiotemporal systems. Building on numerical Gaussian processes, we developed a probabilistic state space model that is linear and has Gaussian distributed states. These properties made it possible to derive the Kalman filter algorithm based on numerical Gaussian process state space models.
The resulting numerical Gaussian process Kalman filter was showcased with a step shaped one-dimensional advection equation where measurements were noisy point evaluations of the solution. In another case study Liouville's equation with integral measurements was studied. In both case studies the numerical GPKF was able to reconstruct the latent state and also estimate the true measurement noise levels through its hyper-parameters.

APPENDIX
Lemmata A and B are taken from [18].

A. Joint distribution of Gaussian variables
If random variables x ∈ R dx and y ∈ R dy have the Gaussian probability distributions x ∼ N (m, P ) (87) y|x ∼ N (Hx + u, R) , than the joint distribution of x, y and the marginal distribution of y are given as x y ∼ N m Hm + u , P P H T HP HP H T + R , y ∼ N Hm + u, HP H T + R .

B. Conditional distribution of Gaussian variables
If the random variables x and y have the joint Gaussian probability distribution than the conditional distribution is For the measurement matrix we have The covariance matrix of the innovation is a bit more elaborate, we have