SymPKF: a symbolic and computational toolbox for the design of parametric Kalman filter dynamics

Recent researches in data assimilation lead to the introduction of the parametric Kalman filter (PKF): an implementation of the Kalman filter, where the covariance matrices are approximated by a parameterized covariance model. In the PKF, the dynamics of the covariance during the forecast step relies on the prediction of the covariance parameters. Hence, the design of the parameter dynamics is crucial while it can be tedious to do this by hand. This contribution introduces a python package, SymPKF, able to compute PKF dynamics for univariate statistics and when the covariance model is parameterized from the variance and the local anisotropy of the correlations. The ability of SymPKF to produce the PKF dynamics is shown on a non-linear diffusive advection (Burgers equation) over a 1D domain and the linear advection over a 2D domain. The computation of the PKF dynamics is performed at a symbolic level, but an automatic code generator is also introduced to perform numerical simulations. A final multivariate example illustrates the potential of SymPKF to go beyond the univariate case.


Introduction
The Kalman filter (KF) (Kalman 1960) is one of the backbone of data assimilation.This filter represent the dynamics of Gaussian distribution all along the analysis and forecast cycles, and takes the form of two equations representing the evolution of the mean and of the covariance of the Gaussian distribution.
While the equations of the KF are simple linear algebra, the large dimension of linear space encountered in the realm of data assimilation make the KF impossible to handle, and this is particularly true for the forecast step.This limitation has motivate some approximation of covariance matrix to make the KF possible.For instance, in ensemble method (Evensen 2009), the covariance matrix is approximated by a sample estimation, where the time evolution of the covariance matrix is then deduced from the forecast of each individual sample.In the parametric Kalman filter (PKF) (Pannekoucke et al. 2016(Pannekoucke et al. , 2018b,a),a), the covariance matrix is approximated by a parametric covariance model, where the time evolution of the matrix is deduced from the time integration of the parameters' evolution equations.
One of the major limitation for the PKF is the design of the parameter evolution equations While it is not mathematically difficult, this step requires the calculus of many terms which is made difficult by hand and might include mistake during the computation.To facilitate the derivation of the parametric dynamics and certify the correctness of the resulting system a symbolic derivation of the dynamics would be welcome.
The goal of the package SymPKF1 (Pannekoucke 2021b) is to facilitate the computation of the PKF dynamics for a particular class of covariance model, the VLATcov model, which is parameterized by the variance and the anisotropy.The symbolic computation of the PKF dynamics relies on a computer algebra system (CAS) able to handle abstract mathematical expression.A preliminary version has been implemented with Maxima2 .However, in order to create and integrated framework which would include the design of the parametric system as well as its numerical evaluation, the symbolic python package Sympy (Meurer et al. 2017) has been preferred for the present implementation.In particular, SymPKF comes with an automatic code generator so to provide an end-to-end exploration of the PKF approach, from the computation of the PKF dynamics to its numerical integration.
The paper is organized as follows.The next section provides the background on data assimilation and introduces the PKF.The Section 3 focuses on the PKF for univariate VLATcov models, in the perspective of its symbolic computation by a CAS.Then, the package SymPKF is introduced in Section 4 from its use on the non-linear diffusive advection (the Burgers' equation) over a 1D domain.A numerical example illustrates the use of the automatic code generator provided in SymPKF .Then, the example of the linear advection over a 2D domain shows the ability of SymPKF to han-dle 2D and 3D domains.The section ends with a simple illustration of a multivariate situation, which also show that SymPKF applies on a system of prognostic equations.The conclusion is given in Section 5.
2 Description of the PKF

Context of the numerical prediction
Dynamics encountered in geosciences is given as a system of partial differential equation where X (t, x) is the state of the system and denotes either a scalar field or multivariate fields in a coordinate system x = (x i ) i∈ [1,d] where d is the dimension the geographical space ; ∂X are the partial derivatives with respect to the coordinate system at any orders, with the convention that order zero denotes the field X itself ; and M denotes the trend of the dynamics.A spatial discretization (e.g. by using finite differences, finite elements, finite volumes, spectral decomposition, etc) transforms Eq. ( 1) into where this time, X (t) is a vector, and M denotes this time the discretization of the trend in Eq. ( 1).Thereafter, X can be seen either as a collection of continuous fields with dynamics given by Eq. (1) or a discrete vector of dynamics Eq. (2).
Because of the sparsity and the error of the observations, the forecast X f is only an estimation of the true state X t , which is known to within a forecast-error defined by e f = X f − X t .This error is often modelled as an unbiased random variable, E e f = 0.In the discrete formulation of the dynamics Eq. ( 2), the forecasterror covariance matrix is given by P f = E e f (e f ) T where the superscript T denotes the transpose operator.Since this contribution is focused on the forecast step, thereafter the upper script f is removed for the sake of simplicity.
We now detail how the error-covariance matrix evolves during the forecast by considering the formalism of the second-order nonlinear Kalman filter.

Second-order nonlinear Kalman filter
A second-order nonlinear Kalman filter (KF2) is a filter that extends the Kalman filter (KF) to the nonlinear situations where the error-covariance matrix evolves tangent-linearly along the trajectory of the mean state and where the dynamics of this mean is governed by the fluctuation-mean interacting dynamics (Jazwinski 1970;Cohn 1993).Hence, we first state the dynamics of the mean under the fluctuation-mean interaction, then the dynamics of the error covariance.Note that the choice of the following presentation is motivated by the perspective of using a computer algebra system to perform the computation.

Computation of the fluctuation-mean interaction dynamics
Because of the uncertainty on the initial condition, the state X is modelized as a Markov process X (t, x, ω) where ω stands for the stochasticity, while X evolves by Eq. ( 1).Hence, ω lies within a certain probability space (Ω, F, P ) where F is a σ−algebra on Ω (a family of subsets of Ω, which contains Ω and which is stable for the complement and the countable union) and P is a probability measure see e.g.(Oksendal 2003, chap.2) ; and where the integer n is either the dimension of the multivariate field X (t, x) or the dimension of its discretized version X (t).The connexion between the Markov process and the parameter dynamics is obtained using the Reynolds averaging technique.So to perform the Reynolds averaging of Eq. ( 1), the first step is to replace the random field by its Reynolds decomposition In this modelling of the random state, E [X ] is the ensemble average or the mean state; e is an error or a fluctuation to the mean, and it is an unbiased random field, E [e] = 0.Then, Eq. ( 1) reads as where η is a control of magnitude introduced to facilitate Taylor's expansion when using a computer algebra system.At the second order, the Taylor's expansion in η of Eq. (3) reads where M and M are two linear operators, the former (the later) refers to the tangent-linear model (the hessian), both computed with respect to the mean state E [X ].The first order expansion is deduced from Eq. (4a) by setting η 2 = 0, which then reads as By setting η to one, the dynamics of the ensemble average is calculated at the second order from the expectation of Eq. (4a) that reads as where ∂e ⊗ ∂e denotes the tensorial product partial derivatives with respect to the spatial coordinates, i.e. terms as ∂ k e∂ m e for any positive integers (k, m).Here, we have used that the partial derivative commutes with the expectation, E [∂e] = ∂E [e], and that E [e] = 0.Because the expectation is a projector, ) is itself.The second term of the right hand side makes appear the retro action of the error onto the ensemble averaged dynamics.Hence, Eq. ( 5) gives the dynamics of the error-mean interaction (or fluctuation-mean interaction).
Note that, the tangent-linear dynamics along the ensemble averaged dynamics Eq. ( 5) is obtained as the difference between the first order Taylor's expansion Eq. (4b) by its expectation, and reads as Now it is possible to detail the dynamics of the error covariance from the dynamics of the error which tangent-linearly evolves along the mean state E [X ].

Computation of the error-covariance dynamics
In the discretized form, the dynamics of the error Eq. ( 6) reads as where M stands for the tangent-linear model So the dynamics of the error-covariance matrix, P = E ee T , is given by (M T is the adjoint of M), or its integrated version where M t←0 is the propagator associated to the time integration of Eq. ( 7), initiated from the covariance P 0 .

Setting of the KF2
Gathering the dynamics of the ensemble mean given by the fluctuation-mean interaction Eq. ( 5) and the covariance dynamics Eq. ( 8) leads to the second-order closure approximation of the extended KF, that is the forecast step equations of the KF2.
Similarly to the KF, the principal limitation of the KF2 is the numerical cost associated with the covariance dynamics Eq. ( 8): living in a discrete world, the numerical cost of Eq. ( 8) dramatically increases with the size of the problem.As an example, for the dynamics of simple scalar field discretized with n grid points, the dimension of its vector representation is n, while the size of the error-covariance matrix scales as n 2 ; leading to a numerical cost of Eq. (8) between O(n 2 ) and O(n 3 ).
We now introduce the parametric approximation of covariance matrices which aims to reduce the cost of the covariance dynamics Eq. (8).

Formulation of the PKF prediction
The parametric formulation of covariance evolution stands as follows.If P(P) denotes a covariance model featured by a set of parameters P = (p i ) i∈I , then there exists a set P f t featuring the forecast error covariance matrix so that P(P f t ) ≈ P f t .Hence, starting from the initial condition P f = P f 0 , if the dynamics of the parameters P f t is known, then it is possible to approximately determine P f t ≈ P(P f t ) without solving Eq. ( 8) explicitly.This approach constitutes the so-called parametric Kalman filter (PKF) approximation, introduced by Pannekoucke et al. (2016Pannekoucke et al. ( , 2018a) (P16, P18).
We now focus on the PKF applied to a particular family of covariance models.

PKF for VLATcov models
This part introduces a particular family of covariance models, parameterized by the fields of variances and of local anisotropy tensor: the VLATcov models (Pannekoucke 2020).What makes this covariance model interesting is that its parameters are related to the error field and thus, it is possible to determine the dynamics of the parameters.So to introduce VLATcov models, we first present the diagnosis of the variance and of the local anosotropy tensor, then we present two examples of VLATcov models and we end the section by the description of the dynamics of the parameters.

Definition of the fields of variance and of local anisotropy tensor
From now, we will focus on the forecast-error statistics, so the upperscript f is removed for the sake of simplicity.Moreover, for a function f , when there is no confusion, the value of f at a point x is written either as f (x) or as f x .
The forecast error being unbiased, E [e] = 0, its variance at a point x is defined as When the error is a random differentiable field, the anisotropy of the two-points correlation function x)e(y)] is featured, from the second order expansion by the local metric tensor g(x), and defined as where ρ x (y) = ρ(x, y) e.g.
The metric tensor is a symmetric positive definite matrix, and it is a 2 × 2 (3 × 3) matrix in a 2D (3D) domain.
Note that it is useful to introduce the local aspect tensor (Purser et al. 2003) whose the geometry goes as the correlation, and defined as the inverse of the metric tensor where the superscript −1 denotes the matrix inverse.What makes the metric tensor attractive, either at a theoretical or at a practical level, is that it is closely related to the normalized error ε = e √ V by (see e.g.(Pannekoucke 2020) for details).
Hence, a VLATcov model is a covariance model characterized by the variance field, and by the anisotropy field, the latter being defined either by the metric-tensor field g or by the aspect-tensor field s.To put some flesh on the bones, two examples of VLATcov models are now presented.

Examples of VLATcov models
The covariance model based on the the heterogeneous diffusion operator of Weaver and Courtier ( 2001) is often introduced in numerical weather or ocean prediction to model heterogeneous correlation functions.This model own the property that, under the local homogenous assumption (that is when the spatial derivatives are negligible) the local aspect tensors of the correlation functions are twice the local diffusion tensors.(Pannekoucke and Massart 2008;Mirouze and Weaver 2010).Hence, by defining the local diffusion tensors as half the local aspect tensors, the covariance model based on the heterogeneous diffusion equation is a VLATcov model.
Another example of heterogeneous covariance model is the heterogenous Gaussian covariance model where ν is a field of symmetric positive definite matrices, and |ν| denotes the matrix determinant.P he.g (V, ν) is a particular case of the class of covariance models deduced from Theorem 1 of Paciorek and Schervish (2004).Again, this covariance model own the property that, under local homogeneous assumptions, the local aspect tensor is approximately given by ν, i.e. for any point x, Hence, as for the covariance model based on the diffusion equation, by defining the field ν as the aspect tensor field, the heterogeneous Gaussian covariance model is a VLATcov model (Pannekoucke 2020).
At this stage, all the pieces of the puzzle are put together to build the PKF dynamics.We have covariance models parameterized from the variance and the local anisotropy, which are both related to the error field: knowing the dynamics of the error leads to the dynamics of the VLATcov parameters.This is now detailed.

PKF prediction step for VLATcov models
When the dynamics of the error e is well approximated from the tangent-linear evolution Eq. ( 6), the connection between the covariance parameters and the error, Eq. ( 9) and Eq. ( 13), makes possible to establish the prediction step of the PKF (Pannekoucke et al. 2018a), which reads as the dynamics of the ensemble average (at the secondorder closure) coupled with the dynamics of the variance and the metric where it remains to replace the dynamics of the error (and its normalized version ε = e/ √ V ) from Eq. ( 6), and where property that the expectation operator and the temporal derivative commutes, has been used to obtain Eq. (16b) and Eq.(16c).
The set of equations ( 16) is at the heart of the numerical sobriety of the parametric approach.In contrast to the matrix dynamics of the KF, the PKF approach is designed for the continuous world, leading to PDEs for the parameters dynamics in place of ODEs Eq. ( 8) for the full matrix dynamics.For the above mentioned scalar fields, introduced is the computation of the algorithmic complexity in section 2.1, the cost of Eq. ( 16) is O(n).Moreover, the dynamics of the parameters sheds light on the nature of the processes governing the dynamics of covariances ; and it does not require any adjoint of the dynamics (Pannekoucke et al. 2016(Pannekoucke et al. , 2018a)).
Note that Eq. ( 16) can be formulated in terms of aspect tensors, thanks to the definition Eq. ( 12): since, sg = I, its time derivative and then where it remains to replace occurrences of g by s −1 in the resulting dynamics of the mean, the variance and the aspect tensor.Hence, the PKF forecast step for VLATcov model is given by either the system Eq.( 16) (in metric), or by its aspect tensor formulation thanks to Eq. ( 17).Whatever the formulation considered, it is possible to carry out the calculations using a formal calculation language.However, even for simple physical processes, the number of terms in formal expressions can become very large, e.g. it is common to have to manipulate expressions with more than a hundred terms.Thus, any strategy that simplifies the assessment of PKF systems in advance can quickly become a significant advantage.
In the following section, we present the splitting method that allows the PKF dynamics to be expressed by bringing together the dynamics of each of the physical processes, calculated individually.

The splitting strategy
When there are several processes in the dynamics Eq. ( 1), the calculation of the parametric dynamics can be tedious even when using a computer algebra system.To better use digital resources, a splitting strategy can be introduced (Pannekoucke et al. 2016(Pannekoucke et al. , 2018a)).
While the theoretical background is provided by the Lie-Trotter formula for Lie derivatives, the well-known idea of time-splitting is easily catched from a first order Taylor expansion of an Euler numerical scheme: The computation of a dynamics over a single time step δt, can be done in two times following the numerical scheme where at order δt, this scheme is equivalent to , that is the Euler step of Eq. ( 18).Because f 1 and f 2 can be viewed as vector fields, the fractional scheme, joining the starting point (at t) to the end point (at t + δt), remains to going through the parallelogram, formed by the sum of the two vectors, along its sides.Since there are two paths joining the extreme points, starting the computation by f 2 is equivalent to starting by f 1 (at order δt), this corresponds to the commutativity of the diagram formed by the parallelogram.
Appendix A shows that a dynamics given by Eq. ( 18) implies a dynamics of the error, the variance, the metric and the aspect, written as a sum of trends.Hence, it is possible to apply a splitting for all these dynamics.
As a consequence for the calculation of the parametric dynamics: calculating the parametric dynamics of Eq. ( 18) is equivalent to calculating separately the parametric dynamics of ∂ t X = f 1 (X ) and ∂ t X = f 2 (X ), then bringing together the two parametric dynamics into a single one by summing the trends for the mean, the variance, the metric or the aspect dynamics.This splitting applies when there are more than two processes and appears as a general method to reduce the complexity of the calculation.

Discussion/intermediate conclusion
However, while the computation of the system Eq.( 16) is straightforward, since it is similar to the computation of the Reynolds equations (Pannekoucke et al. 2018a), it is painful because of the numerous terms it implies, and there is a risk to introduce errors during the computation by hand.
Then, once the dynamics of the parameters is established, it remains to design a numerical code to test whether the uncertainty is effectively well represented by the PKF dynamics.Again, the design of a numerical code is not necessarily difficult but with the numerous terms, the risk of introducing an error is important.
To facilitate the design of the PKF dynamics as well as the numerical evaluation, the package SymPKF has been introduced to perform the VLATcov parameter dynamics and to generate a numerical code used for the investigations (Pannekoucke 2021b).The next section introduces and details this tool.

Symbolic computation of the PKF for VLATcov
In order to introduce the symbolic computation of the PKF for VLATcov model, we consider an example: the diffusive non-linear advection, the Burgers equation, which reads as where u stands for the velocity field and corresponds to a function of the time t and the space of coordinate x, and where κ is a diffusion coefficient (constant here).This example illustrates the workflow leading to the PKF dynamics.It consists in defining the system of equation in Sympy, then to compute the dynamics Eq. ( 16), we now detail these two steps.

Definition of the dynamics
The definition of the dynamics relies on the formalism of sympy as shown in Fig. 1.The coordinate system is first defined as instances of the class Symbols.Note that the time is defined as sympkf.twhile the spatial coordinate is let to the choice of the user, here x.
Then, the function u is defined as an instance of the class Function, as a function of (t, x).
In this example, the dynamics consists in a single equation defined as an instance of the class Eq, but in the general situation where the dynamics is given as a system of equation, the dynamics has to be represented as a python list of equations.
A preprocessing of the dynamics is then performed to determine several important quantities to handle the dynamics: the prognostic fields (functions for which a time derivative is present), the diagnostic fields (functions for which there is no time derivative in the dynamics), the constant functions (functions that only depend on the spatial coordinates), and the constants (pure scalar terms, that are not function of any coordinate).This preprocessing is performed when transforming the dynamics as an instance of the class PDESystem, and whose default string output delivers a summary of the dynamics: for the Burgers' equation, there is only one prognostic function, u(t, x), and one constant, κ.
The prognostic quantities being known, it is then possible to perform the computation of the PKF dynamics, as discussed now.

Computation of the VLATcov PKF dynamics
Thanks to the preprocessing, we are able to determine what are the VLATcov parameters needed to compute the PKF dynamics, that is the variance and the anisotropy tensor associated to the prognostic fields.For the Burgers' equation, the VLATcov parameters are the variance V u and the metric tensor g u = (g u,xx ) or its associated aspect tensor s u = (s u,xx ).Note that, in sympkf, the VLATcov parameters are labeled by there corresponding prognostic fields so to facilitate their identification.This labelling is achieved when the dynamics is transformed as an instance of the class SymbolicPKF.This class is at the core of the computation of the PKF dynamics from Eq. ( 16).
As discussed in Section 2.2.1, the PKF dynamics relies on the second-order flucutation-mean interaction dynamics where each prognostic function is replaced by a stochastic counter-part.Hence, the constructor of SymbolicPKF converts each prognostic functions as a function of an additional coordinate, ω ∈ Ω.For the Burgers' equation, u(t, x) becomes u(t, x, ω).
Since the computation of the second-order  fluctuation-mean interaction dynamics relies on the expectation operator, an implementation of this expectation operator has been introduced in sympkf: it is defined as the class Expectation build by inheritance from the class sympy.Function so to leverage on the computational facilities of sympy.
The implementation of the class Expectation is based on the linearity of the mathematical expectation operator with respect to deterministic quantities, and its commutativity with partial derivatives and integrals with respect to coordinates different from ω e.g. for the Burgers' equation Note that E [u(t, x, ω)] is a function of (t, x) only: the expectation operator converts a random variable into a deterministic variable.
Then, the symbolic computation of the second-order fluctuation-mean interaction dynamics Eq. ( 16a) is performed, thanks to sympy, by following the steps as described in Section 2.2.1.In particular, the computation also leads to the tangent-linear dynamics of the error Eq. ( 6), from which it is possible to compute the dynamics of the variance Eq. ( 16b) and of the metric tensor Eq. (16c) (or its associated aspect-tensor version).Applying these steps, and the appropriate substitutions, is achieved in the back-office when calling the in_metric or in_aspect python's property of an instance of the class SymbolicPKF.This is shown for the Burgers' equation in Fig. 2, where the background computation of the PKF dynamics leads to a list of the three coupled equations corresponding to the mean, the variance and the aspect tensor, similar to the system Eq.( 22) first obtained by Pannekoucke et al. (2018a).
Hence, from SymPKF , for the Burgers' equation, the VLATcov PKF dynamics given in aspect tensor reads as where here s u,xx is the single component of the aspect tensor s u in 1D domains.Note that in the output of the PKF equations, as reproduced in Eq. ( 21), the expectation in the dynamics of the mean is replaced by the prognostic field, that is for the Burgers' equation: E [u] (t, x) is simply denoted by u(t, x).
While the Burgers' equation only contains two physical processes i.e. the non-linear advection and the diffusion, the resulting PKF dynamics Eq. ( 21) makes appear numerous terms, which justifies the use of symbolic computation, as above mentioned.The computation of the PKF dynamics leading to the metric and to the aspect tensor formulation takes about 1s of computation (Intel Core i7-7820HQ CPU at 2.90GHz x 8), while it can take more than one hour by hand.
In this example, the splitting strategy has not been considered to simplify the computation of the PKF dynamics.However, it can be done by considering the PKF dynamics for the advection ∂ t u = −u∂ x u and the diffusion ∂ t u = κ∂ 2 x u, and computed separately, then merged to find the PKF dynamics of the full Burgers' equation.For instance, the Fig. 3 show the PKF dynamics for the advection (first cell) and for the diffusion (second cell), where the output can be trace back in Eq. ( 2) e.g. by the terms in κ for the diffusion.
Thanks to the symbolic computation using the expectation operator, as implemented by the class Expectation, it is possible to handle terms as E ε u ∂ 4 x ε u during the computation of the PKF dynamics.The next section details how these terms are handeled during the computation and the closure issue they bring.

Comments on the computation of the VLATcov
PKF dynamics and the closure issue 4.3.1 Computation of terms E ∂ α ε∂ β ε and their connection to the correlation function 21), are directly connected to the correlation function ρ(x, y) = E [ε(x)ε(y)] whose Taylor expansion is written as However, during its computation, the VLATcov PKF dynamics makes appear terms E ∂ α ε∂ β ε with |α| ≤ |β|, where for any α, ∂ α denotes the derivative with respect to the multi-index α = (α i ) i∈[1,n] , α i denoting the derivative order with respect to the i th coordinate x i of the coordinate system ; and where the sum of all derivative order is denoted by |α| = i α i .The issue it that these terms in E ∂ α ε∂ β ε are not directly connected to the Taylor expansion Eq. ( 22).The interesting property of these terms is that they can be rewords as spatial derivatives of terms in the form E [ε∂ γ ε].More precisely, any term E ∂ α ε∂ β ε can be written from derivative of terms in E [ε∂ γ ε] where |γ| < |α| + |β|, and the term E ε∂ α+β ε (see Appendix B for the proof).So, to replace any term in E ∂ α ε∂ β ε by terms in E [ε∂ γ ε] where |γ| < |α| + |β|, a substitution dictionnary is computed in SymPKF , and stored as the variable subs_tree.The computation of this substitution dictionnary is performed thanks to a dynamical programming strategy.Latter, the integer |α| + |β| is called the order of the term E ∂ α ε∂ β ε .Fig. 4 shows the substitution dictionnary computed for the Burgers' equation.It appears that terms of order lower than 3 can be explicitly written from the metric (or its derivatives) while terms of order larger than 4 cannot: this is known as the closure issue (Pannekoucke et al. 2018a).The term E ε∂ 4 x ε , which features long-range correlations, cannot be related neither to the variance or the metric, and has to be closed.We detail this point in the next section.

Analytical and data-driven closure
A naïve closure for the PKF dynamics Eq. ( 21) would be to replace the unknown term E ε u ∂ 4 x ε u by zero.However, in the third equation that corresponds to the aspec tensor dynamics, the coefficient −3κ of the diffusion term ∂ 2 x s u being negative, it follows that the dynamics of s u numerically explodes at an exponential rate.Of course, because the system represent the uncertainty dynamics of the Burgers' equation Eq. ( 20) that is well-posed, the parametric dynamics should not explodes.Hence, the unknown term E ε u ∂ 4 x ε u is crucial: it can balance the negative diffusion so to stabilize the parametric dynamics.
For the Burgers' equation, a closure for E ε u ∂ 4 x ε u has been previously proposed (Pannekoucke et al. 2018a), given by where the symbols ∼ is used to indicate that this is not an equality but a proposal of closure for the term in the left hand side, and which leads to the closed system (24) The closure Eq. ( 23) results from a local Gaussian approximation of the correlation function.Previous numerical experiments have shown that this closure is well adapted to the Burgers' equation (Pannekoucke et al. 2018a).But the approach that has been followed to find this closure is quite specific, and it would be interesting to design a general way to find such a closure.
In particular, it would be interesting to search a generic way for designing closures that leverages on the symbolic computation, which could be plugged with the PKF dynamics computed from SymPKF at a symbolic level.To do so, we propose an empirical closure which leverages on a data-driven strategy so to hybride the machine learning with the physics, as proposed by Pannekoucke and Fablet (2020) with their neural network generator PDE-NetGen.
The construction of the proposal relies on the symbolic computation shown in Fig. 5: The first step is to consider an analytical approximation for the correlation function.For the illustration, we consider that the local correlation function is well approximated by the quasi-Gaussian function Then, the second step is to perform the computation of the Taylor's expansion of Eq. ( 22) at a symbolic level.This is done thanks to sympy with the method series applied to Eq. ( 25) for δx near the value 0 and at a given order, e.g. for the illustration expansion is computed as the sixth order in Fig. 5.
Then, the identification with the Taylor's expansion Eq. ( 22), leads to the closure .
(26) While it looks like the closure Eq. ( 23), the coefficient are not the same.But this suggests that the closure of E ε u ∂ 4 x ε u can be expanded as ) where a 4 = (a 4 0 , a 4 1 , a 4 2 ) are three unknown reals.A data-driven strategy can be considered to find an appropriate value of a 4 from experiments.This has been investigated by using the automatic generator of neural network PDE-NetGen which bridges the gap between the physics and the machine-learning (Pannekoucke and Fablet 2020), and where the training has lead to the value a 4 ≈ (0.93, 0.75, −1.80) ± (5.1 10 −5 , 3.6 10 −4 , 2.7 10 −4 ) (estimation obtained from 10 runs).Since this proposal is deduced from symbolic computation, it is easy to build some proposals for higher-order unknown terms as it is shown in Fig. 5 for the term E ε u ∂ 5 x ε u .Whatever if the closure has been obtained from an analytical or an empirical way, it remains to compute the closed PKF dynamics to assess its performance.To do so a numerical implementation of the system of partial differential equation has to be introduced.As for the computation of the PKF dynamics, the design of a numerical code can be tedious, with a risk to introduce errors in the implementation due to the numerous terms occurring in the PKF dynamics.This task can be tedious especially in reasearch when To facilitate the research on the PKF, SymPKF comes with a python numerical code generator, which provides an end-to-end investigation of the PKF dynamics.This code generator is now detailed.

Automatic code generation for numerical simulations
While compiled language with appropriate optimization should be important for industrial applications, we chose to implement a pure python code generator which offers a simple research framework for exploring the design of PKF dynamics.It would have been possible code generator already based on sympy (see e.g.Louboutin et al. ( 2019)) but such code generators being domain specific, it was less adapted to the investigation of the PKF for arbitrary dynamics.In place, we consider a finite difference implementation of partial derivatives with respect to spatial coordinates.The default domain to perform the computation is the periodic unit square of dimension the number of spatial coordinates.The length of the domain can be specified along each direction.The domain is regularly discretized along each direction while the number of grid-points can be specified   for each direction.
The finite difference takes the form of an operator F that approximates any partial derivate at a second order of consistency: for any multi-index α, δx 2 is finite.The operator F computed with respect to independent coordinates commute, e.g.
The finite difference of partial derivative with respect to multi-index is computed sequentially e.g.
The finite difference of order α with respect to a single spatial coordinate is the centered finite difference based on α + 1 points.
For instance, Fig. 6 shows how to close the PKF dynamics for the Burgers' equation following P18, and how to build a code from an instance of the class sympkf.FDModelBuilder: it creates the class ClosedPKFBurgers.In this example, the code is executed at run time, but it can also be written in an appropriate python's module for adapting the code to a particular situation or to check the correctness of the generated code.At the end, the instance closed_pkf_burgers of the class ClosedPKFBurgers is created, raising a warning to indicate that the value of constant κ has to be specified before to perform a numerical simulation.Note that it is possible to set the value of kappa as a keyword argument in the class ClosedPKFBurgers.Fig. 6 also show a sample of the generated code with the implementation of the computation of the first order partial derivative ∂ x V u , which appears as a centered finite difference.Then, the sample of code show how the partial derivatives are used to compute the trend of the system of partial differential equations Eq. ( 24).
The numerical integration is handle through the inheritance mechanism: the class ClosedPKFBurgers inherits the integration time loop from the class sympkf.Model as described by the UML dia-gram shown in Fig. 7.
In particular, the class Model contains several time schemes e.g. a fourth-order Runge-Kutta scheme.The details of the instance closed_pkf_burgers of the class ClosedPKFBurgers make appears that the closed system Eq.( 24) will be integrated by using a RK4 time scheme, on the segment [0, D] (here D = 1) with periodic boundaries, and discretized by 241 points.
Thanks to the end-to-end framework proposed in SymPKF , it is possible to perform a numerical simulation based on the PKF dynamics Eq. ( 23).To do so, we set κ = 0.0025 and consider the simulation starting from the Gaussian distribution N (u 0 , P f h ) of mean u 0 (x) = U max [1 + cos(2π(x − D/4)/D)]/2 with U max = 0.5, and of covariance matrix where V h = 0.01U max and l h = 0.02D ≈ 5dx.The time step of the fourth-order Runge-Kutta scheme is dt = 0.002.The evolution predicted from the PKF is shown in Fig. 8 (solid lines).This simulation illustrates the time evolution of the mean (panel a) and of the variance (panel b) ; the panel (c) represents the evolution of the correlation length-scale defined from the aspect tensor as L(x) = s u,xx (x).Note that at time 0, the length-scale field is L(x) = l h .For the illustrations, the variance (the length-scale) is normalized by its initial value V h (l h ).
In order to show the skill of the PKF applied on the Burgers, when using the closure of P18, an ensemble validation is now performed.Note that the code generator of SymPKF can be used for an arbitrary dynamics e.g. the Burgers' equation itself.Hence, a numerical code solving the Burgers' equation is rendered from its symbolic definition.Then an ensemble of 1600 forecasts is computed starting from an ensemble of initial errors at time 0. The ensemble of initial errors is sampled from the Gaussian distribution N 0, P f h of zero mean and covariance matrix P f h .Note that the ensemble forecasting implemented in SymPKF as the method Model.ensemble_forecast(see Fig. 7) leverages on the multiprocessing tools of python, so to use the multiple cores of the CPU, when present.On the computer used for the simulation, the forecasts are performed in parallel on the 8 cores.The ensemble estimation of the mean, the variance and the length-scale are shown in Fig. 8 (dashed lines).Since the ensemble is finite, a sampling noise is visible e.g. on the variance at the initial time that is not strictly equals to V h .In this simulation, it appears that the PKF (solid line) coincide with the ensemble estimation (dashed lines) which shows the ability of the PKF to predict the forecasterror covariance dynamics.Note that the notebook corresponding to the Burgers' experiment is available in the example directory of SymPKF .
While this example shows an illustration of SymPKF in 1D domain, the package also applies in 2D and 3D domains, as presented now.

Illustration of a dynamics in a 2D domain
In order to illustrate the ability of SymPKF to apply in a 2D or a 3D domain, we consider the linear advection of a scalar field c(t, x, y) by a stationary velocity field u = (u(x, y), v(x, y)), which reads as the partial differential equation As for the Burgers' equation, the definition of the dynamics relies on sympy (not shown but similar to the definition of the Burgers' equation as given in Fig. 1).This leads to preprocess the dynamics by creating the instance advection of the class PDESystem, which transform the equation into a system of partial differential equations.In particular, the procedure will diagnose the prognostic functions of a dynamics, here the function c ; the constant functions, here u, v the component of the velocity (u, v) ; for this example there is no constant nor exogenous function.
The calculation of the parametric dynamics is handled by the class SymbolicPKF as shown in the first cell in Fig. 9.The parametric dynamics is a property of the instance pkf_advection of the class SymbolicPKF, and when it is called, the parametric dynamics is computed once time.The parametric dynamics formulated in terms of metric is first computed, see the second cell.For the 2D linear advection, the parametric dynamics is a system of five partial differential equations, as it is shown in the output of the second cell: the dynamics of the ensemble average E [c] which outputs as c for the sake of simplicity (first equation), the dynamics of the variance (second equation) and the dynamics of the local metric tensor (last three equations).In compact form, the dynamics is given by the system which corresponds to the 2D extension of the 1D dynamics first found by Cohn (1993) (Pannekoucke et al. 2016), and validates the computation performed in SymPKF .Due to the linearity of the linear advection Eq. ( 29), the ensemble average Eq.( 30a) is governed by the same dynamics Eq. ( 29).The variance is advected by the flow While both the variance, Eq. (30b), and the metric are advected by the flow, the metric is also deformed by the shear Eq. (30c).This deformation more commonly appears on the dynamics written in aspect tensor form, which is given by where Eq. ( 31c) is similar to the dynamics of the conformation tensor in viscoelastic flow (Bird and Wiest 1995;Hameduddin et al. 2018).This example illustrates a 2D situation, but it runs as well in 3D too.Similarly to the simulation conducted for the Burgers' equation, it is possible to automatically generate a numerical able to perform numerical simulations of the dynamics Eq. ( 31) (not shown here).Hence, this 2D domain example showed the ability of SymPKF to apply in dimensions lager than the 1D.
Before to conclude, we would like to present a preliminary application of SymPKF in a multivariate situation.

Toward the PKF for multivariate dynamics
SymPKF can be used to compute the prediction of the variance and the anisotropy in a multivariate situation.
Note that one of the difficulty with the multivariate situation is that the number of equations increases linearly with the number of fields and the dimension of the domain e.g. for a 1D (2D) domain and two multivariate physical fields, there are two ensemble averaged fields, two variance fields and two (six) metric fields.Of course this is no not a problem when using a computer algebra system as done in SymPKF .
So to illustrate the multivariate situation, only a very simple example is introduced.Inspired from chemical transport models encountered in air quality, we consider the transport, over a 1D domain, of two chemical species, whose the concentrations are denoted by A(t, x) and (B(t, x), and advected by the wind u(x).For the sake of simplicity, the two species interact following a periodic dynamics, leading to the coupled system Thanks to the splitting strategy, the PKF dynamics due to the advection has already be detailed in the previous section (see Section 4.5), so we can focus on the chemical part of the dynamics which is given by the processes Figure 9: Sample of code and Jupyter notebook outputs: system of partial differential equations produced by sympkf when applied to the linear advection Eq. ( 29).
on the right hand side of Eq. ( 32).The PKF of the chemical part is computed thanks to SymPKF , and shown in Fig. 10.This time, and as it is expected, multivariate statistics appear in the dynamics.Here, the dynamics of the cross-covariance V AB = E [e A e B ] is given by the fifth equation.The coupling makes appear unknown terms e.g. the term which appears in the sixth equation.To go further, some research is still needed to explore the dynamics and the modelling of the multivariate cross covariances.A possible direction is to take advantage of the multivariate covariance model based on balance operator as often introduced in variational data assimilation (Derber and Bouttier 1999;Ricci et al. 2005).Note that such multivariate covariance models has been recently considered for the design of the multivariate PKF analysis step (Pannekoucke 2021a).Another way is to consider a data-driven stragy to learn the physics of the unknown terms from a training based on ensembles of forecasts (Pannekoucke and Fablet 2020).To conclude, this example shows the potential of interest of SymPKF to tackle the multivariate situation.Moreover, the example also shows that SymPKF is able to perform the PKF computation for a system of partial differential equations.However all the equations should be prognostic, SymPKF is not able to handle diagnostic equations.

Conclusion
This contribution introduced the package SymPKF that can be used to conduct the research on the parametric Kalman filter prediction step for covariance models parameterized by the variance and the anisotropy (VLATcov), by providing an end-to-end framework: from the equations of a dynamics to the development of a numerical code.
The package has been first introduced by considering the non-linear diffusive advection dynamics, the Burg- ers' equation.In particular this example shown the ability of SymPKF to handle abstract terms e.g. the unclosed terms formulated with the expectation operator.The expectation operator implemented in SymPKF is a key tool for the computation of the PKF dynamics.Moreover, we showed how to handle closure and how to automatically render numerical codes.
For univariate situations, SymPKF applies in 1D domain as well as in 2D and 3D domains.This has been shown by considering the computation of the PKF dynamics for the linear advection equation on a 2D domain.
A preliminary illustration on a multivariate dynamics showed the potential of SymPKF to handle the dynamics of multivariate covariance.But this point has to be further investigated, and this constitutes the main perspective of development.Moreover, to perform a multivariate assimilation cycle with the PKF, the multivariate formulation of the PKF analysis state is needed.A first investigation of the multivariate PKF has been proposed by Pannekoucke (2021a).
In its present implementation, SymPKF is limited to the computation with prognostic equations.It is not possible to consider dynamics based on diagnostic equations while these are often encountered in atmospheric fluid dynamics e.g. the geostrophic balance.This constitutes another topic of research development for the PKF, facilitated by the use of symbolic exploration.
Note that the expectation operator as introduced here can be used to compute Reynolds equations encountered in turbulence.This open new perspectives of use of SymPKF for other applications that could be interesting especially for automatic code generation.

A Splitting for the computation of the parametric dynamics
In this section we show that using a splitting strategy is possible for the design of the parametric dynamics.For this, it is enough to show that given a dynamics written as the dynamics of the error, the variance, the metric and the aspect all write as a sum of trends depending on each processes f 1 and f 2 .We show this starting from the dynamics of the error.Due to the linearity of the derivative operator, the TL dynamics resulting from Eq. (33) writes ∂ t e = f 1 (e) + f 2 (e). (34) which can be written as the sum of two trends ∂ t e 1 = f 1 (e) and ∂ t e 2 = f 2 (e), depending exclusively on f 1 and f 2 respectively.For the variance's dynamics, ∂ t V = 2E [e∂ t e], substitution by Eq. ( 34) leads to where ∂ t V 1 = 2E [ef 1 (e)] and ∂ t V 2 = 2E [ef 2 (e)], depends exclusively on f 1 and f 2 respectively.Then the standard deviation dynamics, obtained by differenciating σ 2 = V as 2σ∂ t σ = ∂ t V , writes as the sum of two trends ∂ t σ 1 = 1 σ ∂ t V 1 and ∂ t σ 2 = 1 σ ∂ t V 2 , depending exclusively on f 1 and f 2 respectively.It results that the dynamics of the normalized error ε = 1 σ e, deduced from the time derivative of e = σε, ∂ t e = ε∂ t σ + σ∂ t ε, writes 37) and also expands as the sum of two trends , again depending exclusively on f 1 and f 2 respectively.For the metric terms g ij = E [∂ i ε∂ j ε], we deduce that the dynamics with where each partial trend depends exclusively on f 1 and f 2 respectively.To end, dynamics of the aspect tensor s is deduced from Eq. ( 17) which expands as where ∂ t s 1 = −s(∂ t g 1 )s and ∂ t ν 2 = −s(∂ t g 2 )s only depend of on f 1 and f 2 respectively.To conclude, the computation of the parametric dynamics for Eq. ( 33), is deduced from the parametric dynamics of ∂ t X = f 1 (X ) and ∂ t X = f 2 (X ) calculated separately, then merged together to obtain the dynamics of the variance Eq. ( 35), of the metric Eq. ( 38) and of the diffusion Eq. ( 39 The derivative with respect to a zero α i is the identity operator.Note that the multi-index forms a semigroup since for two multi-index α and β we can form the multi-index α + β = (α i + β i ) i∈ [1,n] .Now the property Eq.(B.1) can be proven considering the following recurrent process: Asuming that the property is true for all patterns of degree strictly lower to the degree |α| + |β| Without loss of generality we assume α i > 0 and denote δ i = (δ ij ) j∈ [1,n] where δ ij is the Kroenecker symbol (δ ii = 1, δ ij = 0 for j = i).From the formula and from the commutativity of the expectation operator and the partial derivative with respect to the coordinate system, it results that (41) Considereing the terms of the left hand side.In one hand, we observe that the degree of the first term is decreasing to |α|+|β|−1, from the reccurence assumption E ∂ α−δi ε∂ β ε can be expanded as terms of the form E [ε∂ γ ε].On the other hand, the degree of the second term remains the same, |α| + |β|, but with a shift of the derivative order.This shift of the order can be do it again following the same process, leading after iterations to the term E ε∂ α+β ε .

Figure 1 :
Figure 1: Sample of code and Jupyter notebook outputs for the definition of the Burgers dynamics using sympkf.

Figure 2 :
Figure 2: Sample of code and Jupyter notebook outputs: systems of partial differential equations given in metric and in aspect forms, produced by sympkf when applied to the Burgers' equation Eq. (20).

Figure 3 :
Figure 3: Illustration of the splitting strategy which can be used to compute the PKF dynamics and applied here for the Burgers' equation: PKF dynamics of the Burgers' equation can be obtained from the PKF dynamics of the advection (first cell) and of the diffusion (second cell).

Figure 5 :
Figure 5: Example of a symbolic computation leading to a proposal for the closure of the unknown terms of order 4 and 5.

Figure 6 :
Figure 6: Introduction of a closure and automatic generation of a numerical code in SymPKF .

Figure 7 :
Figure 7: UML diagram showing the inheritance mechanism implemented in SymPKF : the class ClosedPKFBurgers inherits from the class Model which implements several time schemes.Here, closed_pkf_burgers is an instance of the class ClosedPKFBurgers.

Figure 8 :
Figure 8: Illustration of a numerical simulation of the PKF dynamics Eq. (23) (solid line), with the mean (panel a), the variance (panel b) and the correlation length-scale (panel c) which is defined, from the component s u,xx of the aspect tensor, by L(x) = s u,xx (x).An ensemble-based validation of the PKF dynamics is shown in dashed line.

Figure 10 :
Figure 10: Illustration of the computation of the PKF dynamics for a simple multivariate situation by using SymPKF .