# Geometry-preserving lie group integrators for differential equations on the manifold of symmetric positive definite matrices

Lucas Drumetz, Alexandre Reiffers-Masson, Naoufal El Bekri, Franck Vermet

► To cite this version:

Lucas Drumetz, Alexandre Reiffers-Masson, Naoufal El Bekri, Franck Vermet. Geometry-preserving lie group integrators for differential equations on the manifold of symmetric positive definite matrices. GSI 2023 - 6th International Conference on Geometric Science of Information 23, Aug 2023, Saint Malo, France. pp.1-18. hal-03815325v3

HAL Id: hal-03815325

https://imt-atlantique.hal.science/hal-03815325v3

Submitted on 17 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geometry-preserving Lie group integrators for differential equations on the manifold of symmetric positive definite matrices*

Lucas Drumetz<fnref n="1" />, Alexandre Reiffers-Masson<fnref n="1" />, Naoufal El Bekri<fnref n="1" />,<fnref n="2" />, and Franck Vermet<fnref n="2" />

**Abstract.** In many applications, one encounters time series that lie on manifolds rather than a Euclidean space. In particular, covariance matrices are ubiquitous mathematical objects that have a non Euclidean structure. The application of Euclidean methods to integrate differential equations lying on such objects does not respect the geometry of the manifold, which can cause many numerical issues. In this paper, we propose to use Lie group methods to define geometry-preserving numerical integration schemes on the manifold of symmetric positive definite matrices. These can be applied to a number of differential equations on covariance matrices of practical interest. We show that they are more stable and robust than other classical or naive integration schemes on an example.

**Keywords:** Lie groups, Differential equations, Symmetric positive definite matrices, Stochastic differential equations

## 1 Introduction

Ordinary Differential Equations (ODEs) arise everywhere in science, and are a fundamental tool to describe continuous-time dynamical systems [16]. However, numerical integration is almost always required to obtain approximate solutions. Most of the time, the variable to integrate lives in a Euclidean space, typically $\mathbb{R}^n$. One can then choose from many methods, ranging from simple explicit/implicit Euler or Runge-Kutta methods to adaptive time step schemes [7].

In a number of situations, one may require that the variable to integrate lies on a manifold [14]. Examples include flows on spheres, rotation or covariance matrices (or other matrix manifolds)[15]... For embedded submanifolds of $\mathbb{R}^n$, though the underlying vector space makes it possible to apply classical ODE integration methods, they cannot guarantee that the numerical solution stays on the manifold at each time step. Providing these guarantees is crucial for subsequent uses of

<footnote n="1">¹ IMT Atlantique, Lab-STICC, UMR CNRS 6285, Brest, France</footnote>

<footnote n="2">² Univ Brest, UMR CNRS 6205, Laboratoire de Mathématiques de Bretagne Atlantique, France</footnote>

<footnote>lucas.drumetz@imt-atlantique.fr</footnote>

<footnote>* This work was supported by Agence Nationale de la Recherche under grant ANR-21-CE48-0005 LEMONADE.</footnote>
the solution, e.g. computing geodesic distances for Riemannian manifolds [1], or simply keeping the structural or physical interpretation of a variable.

Formally, the flow of a smooth vector field on a smooth manifold $\mathcal{M}$ generated by an ODE with the initial condition $x(0) = x_0$, writes:

$$ \frac{dx}{dt} = F|_{x(t)}(t), \qquad (1) $$

where $x \in \mathcal{M}$ and $F|_{x(t)}(t)$ is a (time dependent) tangent vector to $\mathcal{M}$ at $x$, and $F: [0, +\infty[ \to \mathcal{X}(\mathcal{M})$, with $\mathcal{X}(\mathcal{M})$ the set of smooth vector fields on $\mathcal{M}$ [20]. To integrate such ODEs when $\mathcal{M}$ has the additional structure of a Lie group, several frameworks were developed under the umbrella term "Lie group integrators" [8]. Interestingly, these methods can be extended to any smooth manifold acted upon transitively by a Lie group [15], that is to any homogeneous space.

In this paper, we focus more specifically on the manifold of $n \times n$ symmetric positive definite (SPD) matrices, denoted as $Sym_n^+$ [5]. It is the manifold of (non-degenerate) covariance matrices, fundamental for multivariate statistics. Flows of covariance matrices arise in many applications, such as Brain Computer Interfaces (BCI) [25], Diffusion Tensor Image processing [12], finance [9], control [4], or data assimilation [13], to represent the evolution of second order moments of random variables. For example, the second order moments of the solution of Stochastic Differential Equations (SDEs) provide simplified and interpretable representations of stochastic processes, though partial in general. In data assimilation, quantifying and propagating the uncertainty of the state variable is crucial and is done in practice using covariance matrices [21]. Solutions of covariance matrix ODEs which are not SPD are meaningless in terms of statistical interpretation. Thus, we focus on equations similar to (1), where the manifold $\mathcal{M}$ is $Sym_n^+$, and the RHS of (1) is a symmetric matrix (an element of the tangent space of $Sym_n^+$). In spite of the ubiquiteness of covariance matrices, to the best of our knowledge, Lie group integrators have not been considered yet for $Sym_n^+$.

Our contributions are multiple. i) We show that in the case of $Sym_n^+$, for small enough time steps, classical methods actually remain in the manifold, but, with moderately big time steps, the iterates may cross the boundary of $Sym_n^+$, (consisting in positive semidefinite matrices), leading to meaningless solutions or even diverging algorithms. ii) We propose to use a Lie group action of invertible matrices on $Sym_n^+$ to turn the latter into a homogeneous space, that can be used to integrate many equations of interest. iii) From there, we design Lie group versions of the Euler and Runge-Kutta 4 (RK4) methods (applicable to many other schemes) on $Sym_n^+$. iv) We conduct experiments an example ODE on $Sym_n^+$ related to a multivariate SDE. They indicate that our integrators perform better than classical or naive schemes, in particular when the integration step is large.

## 2 Theoretical results

In this section, we provide sufficient conditions on the integration time step $\rho$ of the form $P_{i+1} = P_i + \rho T$ to either stay in or leave $Sym_n^+$. The results are stated in the following theorem:
**Theorem 1.** Let $\mathbf{P} \in \text{Sym}_n^+$, and $\mathbf{T} \in \text{Sym}_n$. We denote as $\lambda_1 \le \dots \le \lambda_n$ the eigenvalues of $\mathbf{P}$ and as $\nu_1 \le \dots \le \nu_n$ the (real) eigenvalues of $\mathbf{T}$. We define a set $S = \{(i,j) \in \mathbb{N}_*^2, 1 \le i,j \le n, i+j = n+1, \nu_i \le 0\}$. Then the following holds:

1. If $\mathbf{T}$ is positive semidefinite, then $\forall \rho \in \mathbb{R}^+, \mathbf{P} + \rho \mathbf{T} \in \text{Sym}_n^+$.

2. When $\mathbf{T}$ has at least one negative eigenvalue, $S$ in nonempty and if

$$ \rho \ge \rho_{\min} = \min_{(i,j) \in S} -\frac{\lambda_j}{\nu_i} $$

then $\mathbf{P} + \rho \mathbf{T} \notin \text{Sym}_n^+$.

3. When $\mathbf{T}$ has at least one negative eigenvalue, if

$$ \rho < \rho_{\max} = -\frac{\lambda_1}{\nu_1}, $$

then $\mathbf{P} + \rho \mathbf{T} \in \text{Sym}_n^+$.

The proof can be found in the appendices. This theorem means that when the vector used for the update happens to be positive semidefinite, then a classical Euler or RK step will remain on $\text{Sym}_n^+$. However, in the general case, that vector may have negative eigenvalues, and we have shown that the next iterate will leave (resp. remain on) $\text{Sym}_n^+$ for time steps that are too large (resp. small enough). In practice, the lower bound is obtained by searching the best value in $S$, whose cardinal is the number of negative eigenvalues of $\mathbf{T}$. Time steps whose values lie in between both bounds may or may not stay on the manifold. Thus, algorithms at least guaranteeing that the trajectory remains on $\text{Sym}_n^+$ are necessary. We stress that even when iterates of classical methods remain on $\text{Sym}_n^+$, the geometry of the manifold is not accounted for, which may lead to low quality solutions nonetheless.

## 3 Background on Lie group integrators

The general idea behind basic Lie group methods relies on the fact that the flow of a simple class of vector fields on the Lie group is easy to compute via the Lie exponential map. The corresponding equations are analogous to the linear ODE $dx/dt = \mathbf{A}x$ in Euclidean spaces. In the case of general ODEs, discretizing by temporarily fixing the vector field of a general equation with a nonconstant $\mathbf{A}$ (depending on $x$ and $t$), an approximate Euler-like scheme can be computed step by step. Higher order schemes such as RK4 require a vector space structure to be able to manipulate vector fields at different locations and times. Then, we have to further translate the ODE from the Lie group to the Lie algebra.

Interestingly, all these methods can be effortlessly extended to smooth manifolds on which we can find a *transitive Lie group action* [20,8] (i.e. homogeneous spaces). Throughout this section, we follow [15] (Chap.2). Here, we limit ourselves to matrix Lie groups for simplicity. Let $\mathcal{M}$ be the smooth manifold, $G$ the matrix Lie group. A smooth map $\Lambda: G \times \mathcal{M} \to \mathcal{M}$ is a Lie group action if and only if

$$ \forall x \in \mathcal{M}, \Lambda(\mathbf{I}, x) = x, \qquad (2) $$

$$ \forall x \in \mathcal{M}, \forall \mathbf{A}, \mathbf{B} \in G, \Lambda(\mathbf{A}, \Lambda(\mathbf{B}, x)) = \Lambda(\mathbf{A}\mathbf{B}, x). \qquad (3) $$
$\Lambda$ is further said to be transitive if

$$ \forall x, y \in \mathcal{M}, \exists \mathbf{A} \in G, \Lambda(\mathbf{A}, x) = y. \qquad (4) $$

This means that any point of $\mathcal{M}$ can be reached from any other using the group action with an element of $G$. Equivalently, the group action has only one orbit. To every Lie group $G$ is associated a Lie algebra $\mathfrak{g}$, which is a vector space (the tangent space to the Lie group at the identity). For matrix Lie groups, the Lie algebra is also a set of matrices. Associated to a transitive Lie group action on a smooth manifold is a *algebra action* that translates the group action into an infinitesimal action giving an element of tangent space to the manifold at every point. It determines the type of equations that can be dealt with. It is defined [15] (Lemma 2.6) as a map $\lambda_* : \mathfrak{g} \times \mathcal{M} \to \mathfrak{X}(\mathcal{M})$ such that, for a given point $x \in \mathcal{M}$:

$$ \lambda_*(\mathbf{A})(x) = \frac{d}{ds} \Lambda(\rho(s), x)|_{s=0}, \qquad (5) $$

where $\rho(s)$ is a smooth curve on $G$, parameterized by a scalar $s$, with initial value $\rho(0) = \mathbf{I}$ and initial speed $\mathbf{A} \in \mathfrak{g}$ ($\rho'(0) = \mathbf{A}$). On matrix groups, $\rho$ can be written as a Taylor expansion:

$$ \rho(s) = \mathbf{I} + s\mathbf{A} + o(s). \qquad (6) $$

Intuitively, this curve represents a direction $\mathbf{A}$ on the Lie algebra towards which we can move infinitesimally from any point in $G$.

The first step is to write the differential equation in terms of the algebra action associated with an adequately chosen group action ($\Lambda$, with the associated $\lambda_*$):

$$ \frac{dx}{dt} = F|_{x(t)}(t) = \lambda_*(\xi(x(t), t))(x(t)), \qquad (7) $$

with the initial condition $x(t_i) = x_i$. $\xi: \mathcal{M} \times \mathbb{R}^+ \to \mathfrak{g}$ is a smooth function.

Thanks to [15] (Lemma 2.7), we know there exists $\mathbf{Y}(t) \in G$ following the differential equation (with $\mathbf{Y}(0) = \mathbf{I}$):

$$ \frac{d\mathbf{Y}}{dt} = \xi(x, t)\mathbf{Y}(t). \qquad (8) $$

such that $\Lambda(\mathbf{Y}(t), x_0)$ is the solution of Eq (7). At this step, we can already design a simple Lie-Euler method by temporarily fixing $\xi$ to its current value $\xi(x_i, t_i)$. Then, the solution of the "frozen" ODE (8) on $G$ is [15] (Theorem 2.8):

$$ \tilde{\mathbf{Y}}(t) = \exp(t\xi(x_i, t_i)), \qquad (9) $$

with exp the matrix (Lie group) exponential. By setting $t = t_{i+1} = t_i + h$, we obtain the next iterate on $G$. We come back to $\mathcal{M}$ using the group action:

$$ x_{i+1} = x(t_i + h) = \Lambda(\tilde{\mathbf{Y}}(t+h), x_i). \qquad (10) $$

The Lie-Euler method is summarized in Algorithm 1. For other schemes requiring
**Algorithm 1 Lie-Euler algorithm**

**Require:** Transitive group action $\Lambda$ of $G$ on $M$, associated algebra action $\lambda_*$, Smooth function $\xi$ from $M$ to $\mathfrak{g}$ such that the ODE on $M$ writes as in 7. Initial condition $x_0 \in M$. Step size $h$. Number of time steps $N$.

$$i = 0$$

$$x_i \leftarrow x_0$$

while $i \le N$ do

$$\mathbf{A} \leftarrow h\xi(x_i, t_i)$$

$$x_{i+1} \leftarrow \Lambda(\exp(\mathbf{A}), x_i)$$

$$i \leftarrow i + 1$$

end while

to combine the current iterate on the manifold and intermediary values of the flow, e.g. RK4, we cannot form linear combinations of evaluations of the vector field for different times $t$ and locations $x$, since those live on different tangent spaces, nor can we add them to points on $M$. Thus, we need an ODE defined on a Euclidean space: the Lie algebra $\mathfrak{g}$. Fortunately, [15] (Lemma 3.1) tells us that the solution of Eq. (8) writes as:

$$\mathbf{Y}(t) = \exp(\boldsymbol{\theta}(t)) \qquad (11)$$

where $\boldsymbol{\theta}(t) \in \mathfrak{g}$ solves an ODE in the Lie algebra:

$$\frac{d\boldsymbol{\theta}}{dt} = \exp_{\boldsymbol{\theta}(t)}^{-1}(\boldsymbol{\xi}(x, t)) \qquad (12)$$

with $\boldsymbol{\theta}(0) = \mathbf{0}$. $\exp_{\boldsymbol{\theta}(t)}^{-1}$ is the inverse of the derivative of the matrix exponential at $\boldsymbol{\theta}(t)$. Its approximation at order 4 (see [15], (Def. 2.18)), for any $\mathbf{A}, \boldsymbol{\theta} \in \mathfrak{g}$ is

$$d \exp_{\boldsymbol{\theta}}^{-1}(\mathbf{A}) = \mathbf{A} + \frac{1}{2}[\boldsymbol{\theta}, \mathbf{A}] + \frac{1}{12}[\boldsymbol{\theta}, [\boldsymbol{\theta}, \mathbf{A}]] - \frac{1}{720}[\boldsymbol{\theta}, [\boldsymbol{\theta}, [\boldsymbol{\theta}, \mathbf{A}]]] + o(||\boldsymbol{\theta}||^4) \quad (13)$$

using the matrix commutator $[\mathbf{A}, \mathbf{B}] = \mathbf{A}\mathbf{B} - \mathbf{B}\mathbf{A}$. On Eq. (12), we can apply a classical RK4 scheme. Each necessary evaluation $\mathbf{K}_i \in \mathfrak{g}$ of the vector field can be done using the algebra action $\lambda_*$ at the right location and time. In the Lie algebra, computing weighted combinations of these evaluations is possible, giving a final vector field $\boldsymbol{\Theta} \in \mathfrak{g}$. Finally, we obtain $x_{i+1}$ on the manifold via:

$$x_{i+1} = x(t_i + h) = \Lambda(\exp(\boldsymbol{\Theta}), x_i). \qquad (14)$$

Error order is guaranteed to be the same as the Euclidean scheme as long as there are enough terms (we gave enough for RK4) in the approximation (13). The so-called Runge-Kutta-Munthe-Kaas (RKMS) 4 method is an adaptation of the general RKMS algorithm [15], and can be extended to any choice of coefficients defining a classical RK method. This method is summarized in Algorithm 2.

## 4 Application to the SPD manifold

Here, we examine a suitable Lie group action to build Lie group integration schemes on $\mathrm{Sym}_n^+$. First, $\mathrm{Sym}_n^+$ is indeed a smooth manifold, whose tangent space
**Algorithm 2** Runge-Kutta-Munthe-Kaas 4 algorithm

**Require:** Transitive group action $\Lambda$ of $G$ on $M$, associated algebra action $\lambda_*$, Smooth function $\xi$ from $M$ to $\mathfrak{g}$ such that the ODE on $M$ writes as in 5. Initial condition $x_0 \in M$. Step size $h$. Number of time steps $N$.

$i = 0$

$x_i \leftarrow x_0$

while $i \le N$ do

$$A_1 \leftarrow h\xi(x_i, t_i)$$

$$K_1 \leftarrow A_1$$

$$A_2 \leftarrow h\xi(\Lambda(\exp(K_1/2), x_i), t_i + h/2)$$

$$K_2 \leftarrow \exp_{K_1/2}^{-1}(A_2)$$

$$A_3 \leftarrow h\xi(\Lambda(\exp(K_2/2), x_i), t_i + h/2)$$

$$K_3 \leftarrow \exp_{K_2/2}^{-1}(A_3)$$

$$A_4 \leftarrow h\xi(\Lambda(\exp(K_3), x_i), t_i)$$

$$K_4 \leftarrow \exp_{K_3}^{-1}(A_4)$$

$$\Theta = K_1/6 + K_2/3 + K_3/3 + K_4/6$$

▷ This operation is possible and stable

because all summands are in $\mathfrak{g}$

$$x_{i+1} \leftarrow \Lambda(\exp(\Theta), x_i)$$

$$i \leftarrow i + 1$$

end while

at each point can be identified with the set of symmetric matrices $\text{Sym}_n$. Thus, any differential equation on $\text{Sym}_n^+$ has a symmetric matrix as a RHS.

We choose the Lie group to be the general linear group $GL_n(\mathbb{R})$. Its Lie algebra is simply $\mathbb{R}^{n \times n}$. The group action we consider is:

$$\Lambda: GL_n(\mathbb{R}) \times \text{Sym}_n^+ \rightarrow \text{Sym}_n^+ \\ (M, P) \mapsto MPM^T. \qquad (15)$$

It is a well known group action, and we can easily check it is indeed transitive. For covariance matrices, it corresponds to the effect of an invertible linear transformation of a random vector on its covariance matrix.

Using (5), we can derive the algebra action $\lambda_* : \mathbb{R}^{n \times n} \times \text{Sym}_n^+ \rightarrow \mathfrak{X}(\text{Sym}_n^+)$:

$$\lambda_*(M)(P) = \frac{d}{ds}\Lambda(Y(s), P)|_{s=0} = MP + PM^T, \qquad (16)$$

where $Y(s) = I + sM + o(s)$ is a smooth curve on the Lie group with $Y(0) = I$ and initial speed $M \in \mathbb{R}^{n \times n}$. Following Eq. (7), we can tackle equations of the form

$$\frac{dP}{dt} = \xi(P, t)P + P\xi(P, t)^T, \qquad (17)$$

with $\xi: \text{Sym}_n^+ \times \mathbb{R}^+ \rightarrow \mathbb{R}^{n \times n}$ any smooth function. Eq. (17) is not very restrictive and many equations of interest can be written this way. For instance, with a constant $\xi$, Eq. (17) governs the dynamics of the covariance of a random variable that propagates via a deterministic linear dynamical system [21]. More complex functions $\xi$ can model more complex situations.
# 5 Case Study

## 5.1 Multivariate Geometric Brownian Motion

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Trajectories of the covariance matrices in $Sym_2^+$ obtained by the different RK4 methods for Eq. (20) using the embedding $(a^c, b^c) \mapsto (a, b, c) \in \mathbb{R}^3$. The boundary of $Sym_2^+$ is given by the surface (in gray) defined by $ab - c^2 = 0$. Left: first case, where the iterates always stay in the manifold but only the Lie-RK4 method is accurate. Right: second case, where the Euclidean method produces iterates outside $Sym_2^+$.</figcaption></figure>

We are interested here in a multivariate generalization of a Geometric Brownian Motion (GBM), given by the (Itô) SDE [3]:

$$dX = \left( \mathbf{A} + \frac{1}{2}\mathbf{B}^2 \right) X dt + \mathbf{B} X dW_t, \qquad (18)$$

with $X$ a random vector of size $n$, $\mathbf{A}, \mathbf{B} \in \mathbb{R}^{n \times n}$ two commuting matrices, such that the eigenvalues of $\mathbf{A} + \frac{1}{2}\mathbf{B}^2$ have a strictly negative real part. $W_t \in \mathbb{R}$ is a Brownian motion. A closed form expression of the trajectories for a deterministic initial condition $\mathbf{x}_0$ is given by $X(t) = \exp(t\mathbf{A} + \mathbf{B}W_t)\mathbf{x}_0$. We can derive ODEs followed by the mean $\mathbf{m}$ (taking expectations in (18)) and covariance matrix $\mathbf{P}$ (using Itô's Lemma [19] on $\mathbf{X}\mathbf{X}^T$, taking expectations, and a few algebraic manipulations, see [11]). They provide a broad summary of the statistics of the process up to order two (much simpler than e.g. the Fokker-Planck equation):

$$\frac{d\mathbf{m}}{dt} = \left( \mathbf{A} + \frac{1}{2}\mathbf{B}^2 \right) \mathbf{m} \qquad (19)$$

$$\frac{d\mathbf{P}}{dt} = \left( \mathbf{A} + \frac{1}{2}\mathbf{B}^2 \right) \mathbf{P} + \mathbf{P} \left( \mathbf{A} + \frac{1}{2}\mathbf{B}^2 \right)^T + \mathbf{B}(\mathbf{P} + \mathbf{m}\mathbf{m}^T)\mathbf{B}^T. \qquad (20)$$

Eq. (20) can indeed be put in the form of Eq. (17), using

$$\xi(\mathbf{P}) = \left( \mathbf{A} + \frac{1}{2}\mathbf{B}^2 \right) + \frac{1}{2}\mathbf{B}(\mathbf{P} + \mathbf{m}\mathbf{m}^T)\mathbf{B}^T \mathbf{P}^{-1}. \qquad (21)$$
We can only detail this example (with $n = 2$) here, but our method also applies e.g. to the covariance of a multivariate Ornstein-Uhlenbeck process [17], or to various Riccati equations encountered in control [4], with suitable choices of $\xi$. Still, we will consider two cases, only differing in the way the matrix $\mathbf{A}$ is defined, and the chosen time step, all other parameters remaining the same.

## 5.2 Choice of numerical values for the SDE (18)

For the case study of Section 5, we chose $n = 2$ to be able to visualize the trajectories in the natural embedding of $\mathrm{Sym}_2^+$ in $\mathbb{R}^3$, see Fig 1. $\mathbf{P}_0$ was generated by using scikit-learn's function `make_spd_matrix` [22], with the same random state in both cases, and multiplying the result by a factor of 0.15. In practice, we have $\mathbf{P}_0 \approx \begin{pmatrix} 0.3383 & -0.0716 \\ -0.0716 & 0.0743 \end{pmatrix}$. We chose commuting $\mathbf{A}$ and $\mathbf{B}$ matrices so as to be able to use the closed form solution given in [3]. We further chose them such that $\mathbf{A} + \frac{1}{2}\mathbf{B}^2$ is negative definite. More in detail, we chose a symmetric

$\mathbf{B} = \begin{pmatrix} -0.4 & 0.1 \\ 0.1 & -0.2 \end{pmatrix}$ To make sure $\mathbf{A}$ and $\mathbf{B}$ commute, we chose $\mathbf{A}$ by diagonalizing $\mathbf{B} = \mathbf{O}\mathbf{D}\mathbf{O}^T$, changing the eigenvalues in $\mathbf{D}$, and then computing $\mathbf{A} = \mathbf{O}\mathbf{D}'\mathbf{O}^T$.

In the first case, we chose $\mathbf{D}' = \begin{pmatrix} -10 & 0 \\ 0 & -2 \end{pmatrix}$, and in the second case we chose

$\mathbf{D}' = \begin{pmatrix} -4 & 0 \\ 0 & -8 \end{pmatrix}$.

We set $t \in [0, 2]$ (resp. $t \in [0, 1.5]$), and show results with 30 (resp. 11) evenly spaced time steps in the first (resp. second) case. In the second case, we note that this yields a time step $h = 0.15$. Using the eigenvalues of $\mathbf{P}_0$ and of the RHS of eq. (20), we obtain $\rho_{\max} = 0.02121$, and $\rho_{\min} = 0.1345$. For the RK4 algorithm, similarly, we obtain $\rho_{\max} = 0.0223$ and $\rho_{\min} = 0.1414$. With our choice of time step, the classical Euler and RK4 methods will then fail in the second case. In practice, the upper bound $\rho_{\max}$ seems quite loose on that example and it seems experimentally that the iterate stays on the manifold as long as $\rho \le \rho_{\min}$.

## 5.3 Results

We compare three explicit RK4 schemes for Eq. (20): a Euclidean scheme, a variant where each step is brought back to the manifold using a Riemannian exponential map, with the affine invariant metric of [5] (Riemmanian RK4), and our method (Lie RK4). We compare them to a reference obtained from a classical RK4 method with a small enough time step so we know the trajectory remains on $\mathrm{Sym}_n^+$ and the integration error is small. We start from a Gaussian initial condition $\mathbf{X}_0 \sim \mathcal{N}(\mathbf{m}_0, \mathbf{P}_0)$. Then, the process converges to a distribution given by a Dirac centered at $\mathbf{0}$. We show the trajectories obtained for the three competing methods embedded in $\mathbb{R}^3$ in Fig. 1 (left: first case, right: second case), and plot two distances between the trajectories and the reference in Fig. 2: the Euclidean distance in $\mathbb{R}^{n \times n}$, as well as the affine invariant distance [5] on $\mathrm{Sym}_2^+$. We see from both Figs 1 and 2 that in the first case, only the Lie-RK4 method is
able to accurately follow the reference dynamics for this value of the integration step. The Euclidean RK4 method is unable to produce a good solution while still producing iterates remaining on the manifold, which may be quite deceptive in applications. The Riemannian RK4 method here actually performs worse than its Euclidean counterpart. Indeed, thus method is not truly intrinsic to the manifold, so a bad step in the Euclidean domain often cannot be completely made up for. In the second case, the configuration is different and the integration step larger, so the result is that Euclidean RK4 produces iterates that immediately leave the manifold, producing a statistically meaningless trajectory. This is why the Riemannian distance cannot be computed for this method in the right panel of Fig 2. In this case, Riemannian RK4 is able to guarantee that the iterates remain in the manifold, but the solution does not follow accurately the true dynamics. Finally, Lie-RK4 produces a very accurate solution guaranteed to stay in the manifold, even in such a challenging situation.

For smaller time steps, Lie RK4 remains the most precise up to a point, then classical RK4 may becomes marginally better than our approach in some cases, probably due to accumulating errors during the additional computations. However, in many applications, the step size is imposed by the problem, e.g. when observation data have a low sampling rate.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Euclidean and Riemannian distances between each integrated trajectory and the reference, on the left and right of each panel, respectively. The case corresponding to the left of Fig. 1 is on the left panel, and the other case is in the right panel. In the latter, the Riemannian distance cannot be computed for the Euclidean method since the iterates are outside the manifold.</figcaption></figure>

# 6 Conclusion

We have presented a Lie group framework to define structure-preserving integration schemes for flows of SPD matrices. Our fully intrinsic integrators keep iterates on the manifold, and provides smaller integration error than classical or naive methods, especially for large time steps. This will be useful in our future work to learn and represent uncertainty in data assimilation [10,18] or controls from observation data when governing equations are unknown (by learning a function $\xi$ matching the data). In such cases, the time step is imposed by data and the training process may lead to ill-conditioned equations.
# 7 Appendices

## 7.1 Proof of Theorem 1

In this section, we prove Theorem 1, which provides a sufficient condition on the integration step $\rho$ for the classical explicit Euler method (or any explicit Runge-Kutta method, for that matter) to yield an iterate that does not belong to $\text{Sym}_n^+$, in the case where the (necessarily symmetric) RHS of the ODE evaluated at the current point has at least one negative eigenvalue. $\mathbf{T}$ has to be understood to be the right hand side of the ODE, computed for the current value of the matrix to be integrated. Besides, we also provide a condition on the time step to guarantee that the next iterate of a classical Euler method remains on the manifold. If explicit RK methods are used instead of explicit Euler methods, then $\mathbf{T}$ denotes the (symmetric) vector field used in the RK update instead and the results still apply.

*Proof.* 1. if $\mathbf{T}$ is positive semidefinite, $\forall \mathbf{x} \in \mathbb{R}^n$ we have

$$ \mathbf{x}^T(\mathbf{P} + \rho\mathbf{T})\mathbf{x} = \mathbf{x}^T\mathbf{P}\mathbf{x} + \rho\mathbf{x}^T\mathbf{T}\mathbf{x} > 0, $$

since $\mathbf{x}^T\mathbf{T}\mathbf{x} \ge 0$. Then $\mathbf{P} + \rho\mathbf{T} \in \text{Sym}_n^+$ for any $\rho \ge 0$.

2. if $\mathbf{T}$ has at least one negative eigenvalue, we use inequalities due to Weyl (1912), [24,23] bounding the eigenvalues of a sum of symmetric matrices by sums of suitable pairs of eigenvalues of each matrix:

**Lemma 1 (Weyl inequalities [24]).** Let $\mathbf{N}, \mathbf{R} \in \text{Sym}_n$, and $\mathbf{M} = \mathbf{R} + \mathbf{N}$. Then denoting as $\lambda_1 \le \dots \le \lambda_n$ the eigenvalues of $\mathbf{R}$, $\nu_1 \le \dots \le \nu_n$ the eigenvalues of $\mathbf{N}$, and as $\eta_1 \le \dots \le \eta_n$ the eigenvalues of $\mathbf{M}$, we have

(a) For any $1 \le i, j, k \le n$ such that $i + j = n + k$:

$$ \eta_k \le \lambda_i + \nu_j. \qquad (22) $$

(b) For any $1 \le i, j, k \le n$ such that $i + j = 1 + k$:

$$ \lambda_i + \nu_j \le \eta_k. \qquad (23) $$

In our case, we are interested in knowing whether the smallest eigenvalue of $\mathbf{P} + \rho\mathbf{T}$ is negative or not (i.e. when at least one eigenvalue is negative). Thus, we apply the Weyl inequalities (22) with $\mathbf{R} = \mathbf{P}$ and $\mathbf{N} = \rho\mathbf{T}$, with $k = 1$, denoting here as $\nu_1 \le \dots \le \nu_n$ the eigenvalues of $\mathbf{T}$, and we obtain

$$ \eta_1 \le \lambda_i + \rho\nu_j. $$

for any $(i, j)$ such that $i + j \le n + 1$. If the RHS is smaller than 0, then $\mathbf{P} + \rho\mathbf{T} \notin \text{Sym}_n^+$. Thus, the only case for which the bound is interesting are those for which $\nu_j \le 0$, otherwise the RHS is positive. From this, requiring the RHS to be negative yields the following lower bound on $\rho$:

$$ \rho \ge -\frac{\lambda_j}{\nu_i}, $$
In particular, we are interested in the tightest bound, so in the tuple $(i, j) \in S$ resulting in the smallest RHS, hence the minimum in the statement of the theorem.

3. We prove this point by contraposition. Let us assume that $\mathbf{P} + \rho \mathbf{T} \notin \text{Sym}_n^+$. Then $\eta_1 \le 0$ (at least). This time, we use the second set of Weyl inequalities (23) with $k=1$. Necessarily, we need $i = j = 1$ (the smallest eigenvalues of $\mathbf{P}$ and $\mathbf{T}$) so that $i + j = 1 + k = 2$. We obtain

$$ \lambda_1 + \rho\nu_1 \le \eta_1 \le 0. $$

Again, isolating $\rho$ yields: $\rho \ge -\frac{\lambda_1}{\nu_1}$.

## 7.2 Proof of the transitivity of the group action (15)

*Proof.* Let $\mathbf{X}, \mathbf{Y} \in \text{Sym}_n^+$, and define $\mathbf{G} = \mathbf{Y}^{\frac{1}{2}}\mathbf{X}^{-\frac{1}{2}} \in GL_n(\mathbb{R})$. Then

$$ \begin{aligned} \Lambda(\mathbf{G}, \mathbf{X}) &= \mathbf{G}\mathbf{X}\mathbf{G}^T = \mathbf{Y}^{\frac{1}{2}}\mathbf{X}^{-\frac{1}{2}}\mathbf{X}(\mathbf{Y}^{\frac{1}{2}}\mathbf{X}^{-\frac{1}{2}})^T \\ &= \mathbf{Y}^{\frac{1}{2}}\mathbf{X}^{-\frac{1}{2}}\mathbf{X}^{\frac{1}{2}}\mathbf{X}^{\frac{1}{2}}\mathbf{X}^{-\frac{1}{2}}\mathbf{Y}^{\frac{1}{2}} \\ &= \mathbf{Y}^{\frac{1}{2}}\mathbf{I}\mathbf{Y}^{\frac{1}{2}} \\ &= \mathbf{Y}. \end{aligned} $$

## 7.3 Computation of the Lie algebra homomorphism (16)

We start from the definition given in Eq. (5), and use the fact that on a matrix group such as $GL_n(\mathbb{R})$, a curve $\rho$ for which $\rho(0) = \mathbf{I}$ has an initial speed that can be defined in the usual Euclidean sense, i.e. $\rho'(0) = \mathbf{M}$, with $\mathbf{M} \in \mathcal{M}_n(\mathbb{R})$. Then, the curve can be written as a Taylor expansion at zero:

$$ \rho(s) = \mathbf{I} + s\mathbf{M} + o(s). $$

Then, working from the definition we have

$$ \begin{aligned} \lambda_*(\mathbf{M})(\mathbf{P}) &\triangleq \frac{d}{ds} \Lambda(\rho(s), \mathbf{P})|_{s=0} \\ &= \frac{d}{ds} (\rho(s) \mathbf{P} \rho(s)^T)|_{s=0} \\ &= \frac{d}{ds} ((\mathbf{I} + s\mathbf{M} + o(s)) \mathbf{P} (\mathbf{I} + s\mathbf{M} + o(s))^T)|_{s=0} \\ &= (\mathbf{M} + o(1)) \mathbf{P} (\mathbf{I} + s\mathbf{M} + o(s)) + (\mathbf{I} + s\mathbf{M} + o(s))^T \mathbf{P} (\mathbf{M} + o(1))^T|_{s=0} \\ &= \mathbf{M}\mathbf{P} + \mathbf{P}\mathbf{M}^T. \end{aligned} $$

## 7.4 Derivation of the covariance equations and corresponding $\xi$ for several equations of interest

Here we derive, for several examples of interest, equations followed by covariance matrices. The first three cases concern SDE on a stochastic process $\mathbf{X}_t \in \mathbb{R}^n$ for
which we are interested in the evolution of the covariance matrix (we drop the time index for brevity)

$$ \mathbf{P} = \mathbb{E}[(\mathbf{X} - \mathbb{E}[\mathbf{X}]) (\mathbf{X} - \mathbb{E}[\mathbf{X}])^T] = \mathbb{E}[\mathbf{X}\mathbf{X}^T] - \mathbb{E}[\mathbf{X}]\mathbb{E}[\mathbf{X}]^T. $$

The last case is different in nature and concerns an optimal control problem in continuous time whose solution involves solving an equation on an SPD matrix.

**Linear deterministic system with stochastic initial condition** The first equation of interest is simply a linear deterministic dynamical system, given by

$$ \frac{d\mathbf{X}}{dt} = \mathbf{A}\mathbf{X}. $$

This describes a stochastic process if the initial condition is given by a probability distribution instead of a deterministic value. In that case, taking the expectation of the solutions for any possible $\mathbf{X}_0$, we get the same ODE on $\mathbf{m}(t) = \mathbb{E}[\mathbf{X}(t)]$.

$$ \frac{d\mathbf{m}}{dt} = \mathbf{A}\mathbf{m}. $$

Then, to obtain a differential equation on $\mathbf{P}$, we first compute the derivative of $\mathbf{X}\mathbf{X}^T$:

$$ \frac{d(\mathbf{X}\mathbf{X}^T)}{dt} = \frac{d\mathbf{X}}{dt}\mathbf{X}^T + \mathbf{X}\frac{d\mathbf{X}^T}{dt} = \mathbf{A}\mathbf{X}\mathbf{X}^T + \mathbf{X}\mathbf{X}^T\mathbf{A}^T. $$

Then we can obtain, taking expectations in the previous equation, subtracting $\mathbf{m}\mathbf{m}^T$ to form the covariance (we can swap derivation and expectation by the Leibniz integral rule for Itô diffusions, since $\mathbf{X}_t$ is absolutely continuous for any $t > 0$ and its density is assumed to be smooth enough in both $\mathbf{x}$ and $t$):

$$ \begin{aligned} \frac{d\mathbf{P}}{dt} &= \frac{d}{dt}(\mathbb{E}[\mathbf{X}\mathbf{X}^T] - \mathbf{m}\mathbf{m}^T) \\ &= \mathbf{A}\mathbb{E}[\mathbf{X}\mathbf{X}^T] + \mathbb{E}[\mathbf{X}\mathbf{X}^T]\mathbf{A}^T - \mathbf{A}\mathbf{m}\mathbf{m}^T - \mathbf{m}\mathbf{m}^T\mathbf{A}^T \\ &= \mathbf{A}(\mathbb{E}[\mathbf{X}\mathbf{X}^T] - \mathbf{m}\mathbf{m}^T) + (\mathbb{E}[\mathbf{X}\mathbf{X}^T] - \mathbf{m}\mathbf{m}^T)\mathbf{A}^T \\ &= \mathbf{A}\mathbf{P} + \mathbf{P}\mathbf{A}^T. \end{aligned} $$

From this, it is clear that taking $\xi(\mathbf{P}, t) \equiv \mathbf{A}$ in (17) corresponds exactly to this equation. In other words, the simplest ODE the group action (15) can handle models the evolution of the covariance matrix of a process passing through a linear dynamical system.

**Multivariate Ornstein-Uhlenbeck process** We now switch to an actual (Itô) SDE, in this case a linear one with a constant diffusion:

$$ d\mathbf{X}_t = \mathbf{A}\mathbf{X}_t dt + \mathbf{B}d\mathbf{W}_t, $$
where $X \in \mathbb{R}^n$, $A \in \mathbb{R}^{n \times n}$, $B \in \mathbb{R}^{n \times n}$, $W_t \in \mathbb{R}^n$ is a multivariate Brownian motion (with independent entries). This is a multivariate generalization of the well-known Ornstein-Uhlenbeck process. Since the expectation of the Brownian motion is zero, the ODE followed by the mean of the process is the same as in the previous example. This entails that the mean of the process will converge to zero as well as long as $A$ has no eigenvalues with a positive real part. There is a closed form solution for both the trajectories and the covariance matrix for a deterministic initial condition (see. e.g. [17]). Even if the initial condition is stochastic, we can derive the ODE governing the evolution of the covariance of the process. To do this, we must first derive the SDE followed by $X X^T$, using Itô's Lemma [19]. We have $(X X^T)_{ij} = X_i X_j$. Applying Itô's Lemma to this function yields:

$$d(X_i X_j) = \sum_{k=1}^{n} \frac{\partial(X_i X_j)}{\partial X_k} dX_k + \frac{1}{2} \sum_{k=1}^{n} \sum_{l=1}^{n} \frac{\partial^2(X_i X_j)}{\partial X_k \partial X_l} dX_k dX_l \\ = X_j dX_i + X_i dX_j + dX_i dX_j.$$

Hence, gathering all these terms for $1 \le i, j \le n$ in a matrix form, we get:

$$d(XX^T) = (dX)X^T + XdX^T + dXdX^T.$$

Replacing $dX$ with its expression, and expanding, we obtain:

$$d(XX^T) = (AXdt + BdW_t)X^T + X(AXdt + BdW_t)^T \\ + (AXdt + BdW_t)(AXdt + BdW_t)^T \\ = AXX^T dt + BdW_tX^T + XX^T A^T dt \\ + XdW_t^T B^T + AXX^T A^T dt^2 + AXdtdW_t^T B^T \\ + BdW_tAXdt + BdW_t dW_t B^T.$$

Using the usual multiplication rules, $dt^2 = dtdW_t = 0$ and $dW_t dW_t^T = dtI$, we obtain

$$d(XX^T) = (AXX^T + XX^T A^T + BB^T)dt + XdW_t^T B^T + BdW_tX^T.$$

Note that this is still an Itô SDE, since by vectorizing $XX^T$ into a $n^2$ dimensional vector, the last SDE can be rewritten in the usual form. To do this, we use the well known property of the Kronecker product (denoted as $\otimes$): $\mathrm{vec}(BVA) = (A \otimes B)\mathrm{vec}(V)$, with $\mathrm{vec}$ the vectorization operator. This yields

$$d(\mathrm{vec}(XX^T)) = ((X \otimes A + A^T \otimes X)X + \mathrm{vec}(BB^T))dt + (I + K)(X \otimes B)dW_t,$$

where $K \in \mathbb{R}^{n^2 \times n^2}$ is the commutator matrix, i.e. the matrix such that $\mathrm{vec}(M^T) = K\mathrm{vec}(M)$. Writing out the solution to this SDE, and taking expectations, we get:

$$\frac{d\mathbb{E}[XX^T]}{dt} = A\mathbb{E}[XX^T] + \mathbb{E}[XX^T]A^T + BB^T.$$
Similarly as in the previous section, forming the covariance matrix yields:

$$ 
\begin{aligned}
\frac{d\mathbf{P}}{dt} &= \frac{d}{dt} (\mathbb{E}[\mathbf{X}\mathbf{X}^T] - \mathbf{m}\mathbf{m}^T) \\
&= \mathbf{A}\mathbb{E}[\mathbf{X}\mathbf{X}^T] + \mathbb{E}[\mathbf{X}\mathbf{X}^T]\mathbf{A}^T + \mathbf{B}\mathbf{B}^T - \mathbf{A}\mathbf{m}\mathbf{m}^T - \mathbf{m}\mathbf{m}^T\mathbf{A}^T \\
&= \mathbf{A}(\mathbb{E}[\mathbf{X}\mathbf{X}^T] - \mathbf{m}\mathbf{m}^T) + (\mathbb{E}[\mathbf{X}\mathbf{X}^T] - \mathbf{m}\mathbf{m}^T)\mathbf{A}^T + \mathbf{B}\mathbf{B}^T \\
&= \mathbf{A}\mathbf{P} + \mathbf{P}\mathbf{A}^T + \mathbf{B}\mathbf{B}^T.
\end{aligned}
$$

Finally, by defining $\xi(\mathbf{P}, t) = \mathbf{A} + \frac{1}{2}\mathbf{B}\mathbf{B}^T\mathbf{P}^{-1}$ in (17) we can use our framework to integrate this equation.

**Multivariate Geometric Brownian Motion** The multivariate GBM SDE is given by Eq. (18). When $\mathbf{A}$ and $\mathbf{B}$ commute, and $\mathbf{A} + \frac{1}{2}\mathbf{B}^2$ has no eigenvalue with a negative real part, [3] provides a closed form solution for an initial value of $\mathbf{x}_0$ under the form:

$$ \mathbf{X}(t) = \exp(t\mathbf{A} + \mathbf{B}W_t)\mathbf{x}_0. $$

With our choice of parameters (normal commuting matrices for $\mathbf{A}$ and $\mathbf{B}$, see Sec. 5), we can obtain [3] a closed form solution on the variance of the process at each time step, and for any initial (deterministic) $\mathbf{x}_0$: $\mathbb{E}[||\mathbf{X}(t)||^2]$:

$$ \mathbb{E}[||\mathbf{X}(t)||^2] = \mathbb{E}[\mathbf{X}(t)^T\mathbf{X}(t)] = ||\exp(t\mathbf{Q})\mathbf{x}_0||^2, $$

where $\mathbf{Q} = \frac{\mathbf{A}+(\mathbf{B}+\mathbf{B}^T)^2}{4}$. Then, regardless of the initial distribution, the process converges to $\mathbf{0}$ (with zero covariance).

Letting $\theta \triangleq \mathbf{A} + \frac{1}{2}\mathbf{B}^2$, the ODE followed by $\mathbf{m}(t)$ is

$$ \frac{d\mathbf{m}}{dt} = \theta\mathbf{m}. $$

To obtain the equation on the covariance matrix, we follow the same method as in the previous section, and from the expression of $d(\mathbf{X}\mathbf{X}^T)$ given by Itô's Lemma, we can compute:

$$ 
\begin{aligned}
d(\mathbf{X}\mathbf{X}^T) &= (\boldsymbol{\theta}\mathbf{X}dt + \mathbf{B}\mathbf{X}dW_t)\mathbf{X}^T + \mathbf{X}(\boldsymbol{\theta}\mathbf{X}dt + \mathbf{B}\mathbf{X}dW_t)^T \\
&\quad + (\boldsymbol{\theta}\mathbf{X}dt + \mathbf{B}\mathbf{X}dW_t)(\boldsymbol{\theta}\mathbf{X}dt + \mathbf{B}\mathbf{X}dW_t)^T \\
&= \boldsymbol{\theta}\mathbf{X}\mathbf{X}^Tdt + \mathbf{B}\mathbf{X}dW_t\mathbf{X}^T + \mathbf{X}\mathbf{X}^T\boldsymbol{\theta}^Tdt \\
&\quad + \mathbf{X}dW_t\mathbf{X}^T\mathbf{B}^T + \boldsymbol{\theta}\mathbf{X}\mathbf{X}^T\boldsymbol{\theta}^Tdt^2 + \boldsymbol{\theta}\mathbf{X}dtdW_t\mathbf{X}^T\mathbf{B}^T \\
&\quad + \mathbf{B}\mathbf{X}dW_tdt\mathbf{X}^T\boldsymbol{\theta}^T + \mathbf{B}\mathbf{X}dW_t^2\mathbf{X}^T\mathbf{B}^T,
\end{aligned}
$$

which yields, after applying the conventions mentioned above:

$$ d(\mathbf{X}\mathbf{X}^T) = (\boldsymbol{\theta}\mathbf{X}\mathbf{X}^T + \mathbf{X}\mathbf{X}^T\boldsymbol{\theta}^T + \mathbf{B}\mathbf{X}\mathbf{X}^T\mathbf{B}^T)dt + (\mathbf{B}\mathbf{X}\mathbf{X}^T + \mathbf{X}\mathbf{X}^T\mathbf{B}^T)dW_t. $$
Writing out the solution and taking the expectation:

$$ \frac{d\mathbb{E}[\mathbf{X}\mathbf{X}^T]}{dt} = \theta\mathbb{E}[\mathbf{X}\mathbf{X}^T] + \mathbb{E}[\mathbf{X}\mathbf{X}^T]\theta^T + \mathbf{B}\mathbb{E}[\mathbf{X}\mathbf{X}^T]\mathbf{B}^T. $$

Then, forming the covariance matrix yields:

$$ \begin{aligned} \frac{d\mathbf{P}}{dt} &= \frac{d}{dt} (\mathbb{E}[\mathbf{X}\mathbf{X}^T] - \mathbf{m}\mathbf{m}^T) \\ &= \theta\mathbb{E}[\mathbf{X}\mathbf{X}^T] + \mathbb{E}[\mathbf{X}\mathbf{X}^T]\theta^T + \mathbf{B}\mathbb{E}[\mathbf{X}\mathbf{X}^T]\mathbf{B}^T - \theta\mathbf{m}\mathbf{m}^T - \mathbf{m}\mathbf{m}^T\theta^T \\ &= \theta(\mathbb{E}[\mathbf{X}\mathbf{X}^T] - \mathbf{m}\mathbf{m}^T) + (\mathbb{E}[\mathbf{X}\mathbf{X}^T] - \mathbf{m}\mathbf{m}^T)\theta^T\mathbf{B}\mathbb{E}[\mathbf{X}\mathbf{X}^T]\mathbf{B}^T \\ &= \theta\mathbf{P} + \mathbf{P}\theta^T + \mathbf{B}(\mathbf{P} + \mathbf{m}\mathbf{m}^T)\mathbf{B}^T, \end{aligned} $$

which is indeed equal to Eq. 20. We notice that for the multivariate GBM, the ODEs on the mean and covariance are coupled, since the equation on the covariance involves the mean. In practice, we first integrate the equation on the mean and plug the solution into the covariance ODE.

Finally, choosing $\xi$ as in Eq. (21) puts this equation within our integration framework.

**Riccati differential equation in continuous time finite horizon linear quadratic optimal control** In this last application, we switch to a linear quadratic optimal control problem with finite horizon and continuous time [4]. We consider a state variable $\mathbf{x} \in \mathbb{R}^n$, initially at $\mathbf{x}(0) = \mathbf{x}_0$, that is subject to a linearly controlled linear dynamical system:

$$ \frac{d\mathbf{x}}{dt} = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u}, $$

where $\mathbf{u} \in \mathbb{R}^n$ is the control variable, and $\mathbf{A}, \mathbf{B} \in \mathbb{R}^{n \times n}$. We want to control the system, e.g. to stabilize the trajectories with $t \in [0, t_f]$, if $\mathbf{A}$ is such that without control, the system diverges. This will happen if $\mathbf{A}$ has at least one eigenvalue with a positive real part. In this case, the goal is then to minimize the quadratic cost function wrt to $\mathbf{u}$, leading to a so called Linear Quadratic Regulator (LQR):

$$ \mathcal{L}(\mathbf{u}) = \frac{1}{2} \left( \|\mathbf{x}(t_f)\|_{\mathbf{Q}}^2 + \int_0^{t_f} (\|\mathbf{x}\|_{\mathbf{Q}}^2 + \|\mathbf{u}\|_{\mathbf{R}}^2) dt \right), $$

depending on the choices of the SPD matrices $\mathbf{Q}_{t_f}, \mathbf{Q}$ and $\mathbf{R}$ to set a tradeoff between the energy of the control and how much we want to push the trajectory of $\mathbf{x}$ towards zero. We let $\|\mathbf{x}\|_{\Sigma}^2 = \mathbf{x}^T\mathbf{\Sigma}\mathbf{x}$ for any $\mathbf{x} \in \mathbb{R}^n$, and $\mathbf{\Sigma} \in \text{Sym}_n^+$.

The closed-loop control solution is given in closed form by $\mathbf{u}(t) = \mathbf{K}(t)\mathbf{x}(t)$, where $\mathbf{K}(t)$ is a time dependent gain, given by

$$ \mathbf{K}(t) = \mathbf{R}^{-1}\mathbf{B}^T\mathbf{P}(t). $$
In the expression of the gain, $\mathbf{P}(t) \in \text{Sym}_n^+$ is an SPD matrix that can be obtained by solving the so-called Riccati differential equation:

$$ \frac{d\mathbf{P}}{dt} = -(\mathbf{A}\mathbf{P} + \mathbf{P}\mathbf{A}^T - \mathbf{P}\mathbf{B}\mathbf{R}^{-1}\mathbf{B}^T\mathbf{P} + \mathbf{Q}), $$

with terminal condition $\mathbf{P}(t_f) = \mathbf{Q}_{t_f}$. This time varying algebraic Riccati equation can be derived from the cost function, either using Pontryagin's maximum principle or dynamic programming, through the Hamilton-Jacobi-Bellman equations [4].

Once again, our proposed framework can handle this differential equation (and its many variants) by choosing $\xi(\mathbf{P}, t) = -\mathbf{A} + \frac{1}{2}\mathbf{P}\mathbf{B}\mathbf{R}^{-1}\mathbf{B}^T - \frac{1}{2}\mathbf{Q}\mathbf{P}^{-1}$.

## 7.5 Riemannian metrics on $\text{Sym}_n^+$

In this section, we simply give the expression of the Riemannian distance on $\text{Sym}_n^+$ used in this paper, i.e. the affine invariant distance. We refer to [2],[5],[6] for more details on this and other distances on $\text{Sym}_n^+$ and the associated Riemannian metrics.

The affine-invariant (sometimes called Fisher-Rao) distance [5,6] is defined as

$$ d_A(\mathbf{P}_1, \mathbf{P}_2) = \| \log(\mathbf{P}_1^{-1/2}\mathbf{P}_2\mathbf{P}_1^{-1/2}) \|_F. $$

The corresponding Riemann exponential map $\text{Exp}_P^A$ at $\mathbf{P} \in \text{Sym}_n^+$, used in the Riemannian-RK4 method can be obtained as as [5]:

$$ \text{Exp}_P^A(\boldsymbol{\Sigma}) = \mathbf{P}^{1/2} \exp(\mathbf{P}^{-1/2} \boldsymbol{\Sigma} \mathbf{P}^{-1/2}) \mathbf{P}^{1/2}. $$
References

1. Absil, P.A., Mahony, R., Sepulchre, R.: Optimization algorithms on matrix manifolds. In: Optimization Algorithms on Matrix Manifolds. Princeton University Press (2009)

2. Arsigny, V., Fillard, P., Pennec, X., Ayache, N.: Geometric means in a novel vector space structure on symmetric positive-definite matrices. SIAM journal on matrix analysis and applications 29(1), 328–347 (2007)

3. Barrera, G., Högele, M., Pardo, J.: Cutoff stability of multivariate geometric brownian motion. arXiv preprint arXiv:2207.01666 (2022)

4. Başar, T., Olsder, G.J.: Dynamic noncooperative game theory. SIAM (1998)

5. Bhatia, R.: Positive definite matrices. Princeton university press (2009)

6. Bonnabel, S., Sepulchre, R.: Riemannian metric and geometric mean for positive semidefinite matrices of fixed rank. SIAM Journal on Matrix Analysis and Applications 31(3), 1055–1070 (2010)

7. Butcher, J.C.: Numerical methods for ODEs. John Wiley & Sons (2016)

8. Celledoni, E., Marthinsen, H., Owren, B.: An introduction to lie group integrators–basics, new developments and applications. Journal of Computational Physics 257, 1040–1061 (2014)

9. Dellaportas, P., Pourahmadi, M.: Large time-varying covariance matrices with applications to finance. Technical Report, Athens University of Economics, Department of Statistics (2004)

10. Dridi, N., Drumetz, L., Fablet, R.: Learning stochastic dynamical systems with neural networks mimicking the euler-maruyama scheme. In: 2021 29th European Signal Processing Conference (EUSIPCO). pp. 1990–1994. IEEE (2021)

11. Drumetz, L., Reiffers-Masson, A., Bekri, N.E., Vermet, F.: Geometry-preserving lie group integrators for differential equations on the manifold of symmetric positive definite matrices (2022) https://arxiv.org/abs/2210.08842

12. Dryden, I.L., Koloydenko, A., Zhou, D.: Non-euclidean statistics for covariance matrices, with applications to diffusion tensor imaging. The Annals of Applied Statistics 3(3), 1102–1123 (2009)

13. Evensen, G., Vossepoel, F.C., van Leeuwen, P.J.: Data Assimilation Fundamentals: A Unified Formulation of the State and Parameter Estimation Problem. Springer Nature (2022)

14. Hairer, E.: Solving differential equations on manifolds. Lecture Notes, Université de Geneve (2011)

15. Iserles, A., Munthe-Kaas, H.Z., Nørsett, S.P., Zanna, A.: Lie-group methods. Acta numerica 9, 215–365 (2000)

16. Jordan, D., Smith, P.: Nonlinear ordinary differential equations: an introduction for scientists and engineers. OUP Oxford (2007)

17. Meucci, A.: Review of statistical arbitrage, cointegration, and multivariate ornstein-uhlenbeck. Cointegration, and Multivariate Ornstein-Uhlenbeck (2009)

18. Nguyen, D., Ouala, S., Drumetz, L., Fablet, R.: Variational deep learning for the identification and reconstruction of chaotic and stochastic dynamical systems from noisy and partial observations. arXiv preprint arXiv:2009.02296 (2020)

19. Oksendal, B.: Stochastic differential equations: an introduction with applications. Springer Science & Business Media (2013)

20. Owren, B.: Lie group integrators. In: Discrete Mechanics, Geometric Integration and Lie–Butcher Series, pp. 29–69. Springer (2018)

21. Pannekoucke, O., Ricci, S., Barthelemy, S., Ménard, R., Thual, O.: Parametric kalman filter for chemical transport models. Tellus A: Dynamic Meteorology and Oceanography 68(1), 31547 (2016)
22. Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al.: Scikit-learn: Machine learning in python. the Journal of machine Learning research 12, 2825–2830 (2011)

23. Tao, T.: 254A, Notes 3a: Eigenvalues and sums of Hermitian matrices. https://terrytao.wordpress.com/2010/01/12/254a-notes-3a-eigenvalues-and-sums-of-hermitian-matrices/ (2010), [Online; last accessed 10-May-2023]

24. Weyl, H.: Das asymptotische verteilungsgesetz der eigenwerte linearer partieller differentialgleichungen (mit einer anwendung auf die theorie der hohlraumstrahlung). Mathematische Annalen **71**(4), 441–479 (1912)

25. Zanini, P., Congedo, M., Jutten, C., Said, S., Berthoumieu, Y.: Transfer learning: A riemannian geometry framework with applications to brain–computer interfaces. IEEE Transactions on Biomedical Engineering **65**(5), 1107–1116 (2017)