Numerical low-rank approximation of matrix differential equations

TL;DR

This paper introduces a new splitting integrator for large-scale matrix differential equations that efficiently handles stiffness and preserves key properties of solutions, demonstrated through numerical examples including weather modeling.

Contribution

A novel splitting integrator based on dynamical low-rank approximation is proposed for stiff matrix differential equations, improving efficiency and property preservation.

Findings

Handles stiffness efficiently

Preserves symmetry and positive semidefiniteness

Effective in weather simulation scenarios

Abstract

The efficient numerical integration of large-scale matrix differential equations is a topical problem in numerical analysis and of great importance in many applications. Standard numerical methods applied to such problems require an unduly amount of computing time and memory, in general. Based on a dynamical low-rank approximation of the solution, a new splitting integrator is proposed for a quite general class of stiff matrix differential equations. This class comprises differential Lyapunov and differential Riccati equations that arise from spatial discretizations of partial differential equations. The proposed integrator handles stiffness in an efficient way, and it preserves the symmetry and positive semidefiniteness of solutions of differential Lyapunov equations. Numerical examples that illustrate the benefits of this new method are given. In particular, numerical results for the efficient simulation of the weather phenomenon El Ni\~no are presented.

Figures

Tables

Table 1. Table 1: Defects in symmetry for Algorithm 1 applied to the DLE of Section 3.2 for d = 400 𝑑 400 d=400 . The defect ( 15 ) is displayed as a function of the approximation rank (rk) and the number of time steps (ns) at T = 0.1 𝑇 0.1 T=0.1 .

	ns = $2$	ns = $2^4$	ns = $2^7$	ns = $2^{11}$	ns = $2^{13}$
rk = 2	4.7294e-16	1.6989e-16	1.5067e-16	7.3221e-16	5.1379e-16
rk = 4	7.7743e-16	3.8092e-16	4.9242e-16	1.0861e-15	2.7226e-15
rk = 6	1.0498e-15	3.4863e-16	5.0365e-16	3.2494e-16	3.9935e-16
rk = 8	1.0781e-15	4.1657e-16	3.9016e-16	2.6837e-15	1.2762e-14
rk = 10	8.2787e-16	7.7903e-16	1.3127e-15	1.7855e-15	2.8277e-15
rk = 12	1.1677e-15	4.6731e-16	6.3586e-16	2.1097e-15	6.9804e-15
rk = 14	1.8577e-15	5.5461e-16	8.0892e-16	1.9972e-15	2.5059e-15

Table 2. Table 2: Defects in symmetry for Lie splitting of Section 2 applied to the DLE of Section 3.2 for d = 400 𝑑 400 d=400 without symmetry-preserving modification. The defect ( 15 ) is displayed as a function of the approximation rank (rk) and the number of time steps (ns) at T = 0.1 𝑇 0.1 T=0.1 .

	ns = $2$	ns = $2^4$	ns = $2^7$	ns = $2^{11}$	ns = $2^{13}$
rk = 2	9.5178e-03	3.2545e-03	4.2299e-04	4.3455e-05	5.2796e-06
rk = 4	1.2097e-02	3.8332e-03	2.4587e-04	2.0065e-05	2.3394e-06
rk = 6	2.3350e-04	4.1849e-03	2.6199e-04	2.0735e-05	2.3876e-06
rk = 8	1.2875e-05	8.0839e-04	1.9463e-04	1.5981e-05	1.8315e-06
rk = 10	2.5457e-07	1.5020e-04	1.0291e-04	8.1713e-06	9.2641e-07
rk = 12	1.3272e-09	5.9891e-05	4.7831e-05	3.9946e-06	4.5438e-07
rk = 14	1.0354e-12	5.5624e-06	1.8540e-05	1.7506e-06	1.9895e-07

Table 3. Table 3: Defects in positive semidefiniteness for Algorithm 1 applied to the DLE of Section 3.2 for d = 400 𝑑 400 d=400 . The defect ( 16 ) is displayed as a function of the approximation rank (rk) and the number of time steps (ns) at T = 0.1 𝑇 0.1 T=0.1 .

	ns = $2$	ns = $2^4$	ns = $2^7$	ns = $2^{11}$	ns = $2^{13}$
rk = 2	4.9230e-15	2.7584e-15	2.7871e-15	4.3101e-15	2.3554e-15
rk = 4	4.1822e-15	4.6147e-15	3.2650e-15	5.4743e-15	3.2069e-15
rk = 6	5.9181e-15	1.9928e-15	2.5106e-15	3.0357e-15	4.0661e-15
rk = 8	5.1317e-15	2.6077e-15	1.7991e-15	3.0388e-15	6.6778e-15
rk = 10	7.9000e-15	3.3707e-15	2.7210e-15	3.3223e-15	2.3704e-15
rk = 12	5.7540e-15	2.2482e-15	2.0631e-15	2.6968e-15	5.7000e-15
rk = 14	4.2693e-15	7.0395e-15	2.6659e-15	2.0693e-15	3.2132e-15

Table 4. Table 4: Error of the covariance matrix 𝒫 ​ ( T ) 𝒫 𝑇 \mathcal{P}(T) in ( 21 ) of size 3900 × 3900 3900 3900 3900\times 3900 in Frobenius norm ( 14 ) obtained by the low-rank integrator described in Algorithm  1 with step size τ = 0.1 𝜏 0.1 \tau=0.1 days. The reference solution is computed with step size τ = 0.001 𝜏 0.001 \tau=0.001 days.

2013			2015
weeks	rank	error	weeks	rank	error
1	10	3.1117e-05	1	10	5.0599e-05
	100	8.8863e-07		100	1.4408e-06
2	10	3.1229e-05	2	10	5.3995e-05
	100	1.7762e-06		100	2.8798e-06
3	10	3.1687e-05	3	10	5.5525e-05
	100	2.6620e-06		100	4.3158e-06

Equations

$$\dot{X}(t)=AX(t)+X(t)A^{\mathsf{T}}+G(t,X(t)),X(t_0)=X_0,t\in [t_0,T],$$

$$X(t)={\mathrm{e}}^{(t-t_0)A}X(t_0){\mathrm{e}}^{(t-t_0)A^{\mathsf{T}}}+{\int }_{t_0}^t{\mathrm{e}}^{(t-s)A}G(s,X(s)){\mathrm{e}}^{(t-s)A^{\mathsf{T}}}\text{d}s$$

$$\dot{M}(t)$$

$$\dot{N}(t)$$

$${\Phi }_{\tau }^A(Z)=M(t_0+\tau )={\mathrm{e}}^{\tau A}Z{\mathrm{e}}^{\tau A^{\mathsf{T}}}.$$

$${\mathcal{L}}_{\tau }Z={\Phi }_{\tau }^G\circ {\Phi }_{\tau }^A(Z),$$

$${\mathcal{S}}_{\tau }Z={\Phi }_{\tau /2}^A\circ {\Phi }_{\tau }^G\circ {\Phi }_{\tau /2}^A(Z).$$

$$M_1={\mathrm{e}}^{\tau A}US{V}^{\mathsf{T}}{\mathrm{e}}^{\tau A^{\mathsf{T}}}.$$

$$\Vert \dot{Y}(t)-G(t,Y(t)){\Vert }_F=\min ,$$

$$\dot{Y}(t)=P(Y(t))G(t,Y(t)),$$

$$\dot{Y}=G(t,Y)V{V}^{\mathsf{T}}-U{U}^{\mathsf{T}}G(t,Y)V{V}^{\mathsf{T}}+U{U}^{\mathsf{T}}G(t,Y).$$

$$\dot{X}(t)$$

$$X(t_0)$$

$$\dot{M}(t)=AM(t)+M(t)A^{\mathsf{T}},$$

$$\dot{N}(t)=Q,$$

$$U_1{\hat{S}}_1=U_0{S}_0+\tau Q{U}_0.$$

$${\tilde{S}}_0={\hat{S}}_1-\tau U_1^{\mathsf{T}}Q{U}_0.$$

$$S_1={\tilde{S}}_0{U}_0^{\mathsf{T}}{U}_1+\tau U_1^{\mathsf{T}}Q{U}_1.$$

$$S_1=U_1^{\mathsf{T}}{U}_0{S}_0{U}_0^{\mathsf{T}}{U}_1+\tau U_1^{\mathsf{T}}Q{U}_1$$

$${\partial }_{t}w=\Delta w,w{|}_{\partial \Omega }=0$$

$$\Vert Z{\Vert }_F=\frac{1}{d}\sqrt{\sum _{i,j=1}^d{Z}_{ij}^2},Z\in {\mathbb{R}}^{d\times d}.$$

$$d_{\text{sym}}=\frac{\Vert Y-Y^{\mathsf{T}}{\Vert }_F}{\Vert Y_{\text{ref}}{\Vert }_F},$$

$$d_{\text{psd}}=\frac{\Vert Y-\hat{Y}{\Vert }_F}{\Vert Y_{\text{ref}}{\Vert }_F}$$

$$X_{n+1}=X_n+\tau \left(A{X}_{n+1}+X_{n+1}{A}^{\mathsf{T}}+Q\right).$$

$$0=(X_n+\tau Q)+\left(\tau A-\frac{1}{2}I\right)X_{n+1}+X_{n+1}\left(\tau A^{\mathsf{T}}-\frac{1}{2}I\right).$$

$$0=\tilde{A}X+X{\tilde{A}}^{\mathsf{T}}+\tilde{B}{\tilde{B}}^{\mathsf{T}},$$

$$X=\tilde{V}Y{\tilde{V}}^{\mathsf{T}}.$$

$$\frac{\Vert \tilde{A}{X}^m+X^m{\tilde{A}}^{\mathsf{T}}+\tilde{B}{\tilde{B}}^{\mathsf{T}}{\Vert }_2}{2\Vert \tilde{A}{\Vert }_F\Vert Y^m{\Vert }_F+\Vert \tilde{B}{\Vert }_F^2}\le \text{tol},$$

$$\tilde{A}=\tau A-\frac{1}{2}I\text{and}\tilde{B}=[\sqrt{\tau }R{Z}_n],$$

$$\text{d}X(t)=\mathcal{A}X(t)+F(t),t\in [t_0,T],$$

$$\text{d}\hat{X}(t)$$

$$\text{d}\mathcal{P}(t)$$

Full text

1 Introduction

Matrix differential equations arise quite naturally in many applications in science and engineering. Perhaps the most studied matrix differential equations are differential Riccati (DREs) and differential Lyapunov equations (DLEs) as they play a crucial role in many applications. For instance, they arise in optimal control problems like linear quadratic regulator and linear quadratic Gaussian design, $H_{\infty}$ control of linear time varying systems, optimal filtering, differential games, model reduction of linear time varying systems, damping optimization in mechanical systems, and control of shear flows subject to stochastic excitations. For an overview and more details, we refer to [1, 3, 28].

In the literature, there is a large variety of approaches to compute the solution of DREs and DLEs, see, e.g., [9, 11]. However, due to efficiency reasons, standard methods are not suitable for large-scale problems, in general. A major source of large-scale problems are partial differential equations (PDEs). In two and three space dimensions, their spatial discretization leads to systems with a significant number of degrees of freedom. The system matrices have often a particular structure, and they are usually sparse. In general, the arising differential equation are stiff which in turn requires integrators that handle the stiffness in an efficient way.

The problem of solving large-scale DREs/DLEs has recently received considerable attention. Various integrators based on low-rank approximations have been developed. In particular, in [5, 6] the authors proposed efficient numerical methods which are based on a matrix valued version of backward differentiation formulas (BDFs) and Rosenbrock methods, respectively. Further, a low-rank splitting method for large-scale DREs has been introduced in [32]. A different approach to compute low-rank approximations to large-scale matrix differential equations was proposed in [21]. It is the so-called dynamical low-rank approximation. Its strength relies on the fact that a low-rank approximation to the solution is computed by solving differential equations only for its low-rank factors, see also [20, 22, 25].

The performance of low-rank based methods relies on the decay of the singular values of the solutions. This phenomenon has been deeply studied and is frequent in applications, see, e.g., [4, 30] and references therein.

In this paper, we focus on a new approach for solving DREs and DLEs that arise in PDE based models. Following the ideas of [7] and [22] we develop an efficient algorithm based on a splitting integrator. We split the vector field into a linear stiff and a non-linear non-stiff part. The stiff problem is efficiently integrated by an exponential integrator. We employ here the Leja method which proves to be very competitive, see [7]. We combine this integrator with the dynamical low-rank approximation for the non-stiff part. In particular, we employ the projector-splitting scheme, proposed in [22]. This choice makes our low-rank integrator very robust with respect to small singular values. In this work we keep the rank fixed during time integration. Changing dynamically the rank would require a strategy for adapting the rank, based on a control of the corresponding error. Such a rank-adaptive approach is feasible but beyond the scope of our work. However, we illustrate the behaviour of our integrator with respect to different approximation ranks by numerical experiments.

The paper is organized as follows. In Section 2 we introduce the class of matrix differential equations considered and we propose our low-rank numerical integrator. In Sections 3 and 4, we specialize our integrator for DLEs and DREs, and we discuss a possible extension to generalized DREs. We give details on the implementation and present numerical results that support our chosen approach. Finally, some conclusions are given.

2 Matrix differential equations and numerical integrators

In this work we are interested in the numerical solution of the following class of matrix differential equations

[TABLE]

where $X(t)$, $A\in\mathbb{R}^{d\times d}$, $G:\mathbb{R}\times\mathbb{R}^{d\times d}\rightarrow\mathbb{R}^{d\times d}$, and $(\cdot)^{\mathsf{T}}$ denotes the transpose. This class includes, among others, differential Lyapunov and Riccati equations. In particular, we consider here problems stemming from PDEs. Therefore, the matrix $A$ is typically the spatial discretization of a differential operator. The function $G$ on the other hand is assumed to be non-stiff. By means of the variation-of-constants formula, the solution of (1) can be written as

[TABLE]

for $t_{0}\leq t\leq T$.

The stiffness of (1) is a limitation for most of the standard numerical integrators. Explicit methods have to use tiny step sizes whereas the use of implicit methods can result in high computational cost. Both approaches are thus not efficient. As a remedy, we design here a numerical integrator which is able to handle the stiff part of the solution (2) in an exact and efficient way. The key idea is the use of splitting methods. For an introduction to this class of numerical integrators we refer to [16]. The structure of (1) suggests a splitting into two terms, where we split the stiff linear part from the non-stiff non-linearity. Following this approach, two subproblems arise:

[TABLE]

where $t\in[t_{0},T]$. A clear advantage of this splitting approach is that the stiffness is only present in the first subproblem (3a), whereas the second equation (3b) is a non-stiff ordinary differential equation. Therefore, we have more freedom for the numerical integration of (3b). Depending on the way of recombining the partial flows of (3), we obtain splitting methods with different orders of convergence.

Let us denote with $\Phi_{\tau}^{A}(Z)$ the exact solution of (3a) at $t_{0}+\tau$ with initial value $M(t_{0})=Z$, i.e.,

[TABLE]

Further, let $\Phi^{G}_{\tau}(Z)$ denote the exact solution of (3b) with initial value $N(t_{0})=Z$, i.e., $\Phi_{\tau}^{G}(Z)=N(t_{0}+\tau)$. The simplest splitting integrator for solving (1) is then given by

[TABLE]

where $Z=X_{0}$. The result $\mathcal{L}_{\tau}X_{0}$ is the numerical approximation to $X(t_{0}+\tau)$. This is a first-order method and called Lie splitting in the literature. A second-order method, the so-called Strang splitting, is given by the symmetric formula

[TABLE]

The numerical evaluation of (4) can be computationally expensive if the dimension $d$ is large. In this case we propose to approximate the action of the flow $\Phi^{A}_{\tau}$ by means of a low-rank approximation. Note that the differential equation (3a) is rank preserving, i.e, the rank of the solution $M(t)$ does not depend on time (see [18, Lemma 1.22]). A rank-$r$ decomposition ($r\leq d$) of the initial data $Z$ is typically given by a truncated singular value decomposition (SVD) of the form $USV^{\mathsf{T}}$, where $U,V\in\mathbb{R}^{d\times r}$ have orthonormal columns and $S\in\mathbb{R}^{r\times r}$ is diagonal. Using this approximation for the initial data $Z$, we obtain a rank-$r$ approximation $M_{1}$ to the exact solution $M(t_{0}+\tau)=\Phi_{\tau}^{A}(Z)$ by

[TABLE]

Note that the action of the matrix exponential $\mathrm{e}^{\tau A}$ to the skinny matrices $U$ and $V$ can be efficiently computed, for instance, by Taylor interpolation [2], interpolation at Leja points [7], and Krylov subspace methods [15, 29]. We will compute these actions here with the Leja method. The resulting matrices $\mathrm{e}^{\tau A}U$ and $\mathrm{e}^{\tau A}V$ can be orthogonalized by means of a QR decomposition such that $\mathrm{e}^{\tau A}U=\widetilde{U}R$, where $\widetilde{U}\in\mathbb{R}^{d\times r}$ has orthonormal columns and $R~{}\in~{}\mathbb{R}^{r\times r}$. Similarly we have $\mathrm{e}^{\tau A}V=\widetilde{V}P$. After defining $\widetilde{S}=RSP^{\mathsf{T}}$, a rank-$r$ solution of (3a) is given as $M_{1}=\widetilde{U}\widetilde{S}\widetilde{V}^{\mathsf{T}}$.

For the solution of (3b) we make use of the dynamical low-rank approach proposed in [21] (see also [20, 22, 25]). This is a differential equation based approach to efficiently compute low-rank approximations to time dependent large matrices or to solutions of large-scale matrix differential equations. Let us denote by $\mathcal{M}_{r}=\mathcal{M}_{r}^{d\times d}$ the manifold of rank-$r$ matrices of dimension $d\times d$. At any time $t$, a rank-$r$ approximation to the solution $N(t)\in\mathbb{R}^{d\times d}$ of (3b) is a matrix $Y(t)\in\mathcal{M}_{r}$ defined by the following condition: for every $t$, the matrix $\dot{Y}(t)\in\mathcal{T}_{Y(t)}\mathcal{M}_{r}$ is such that

[TABLE]

where $\mathcal{T}_{Y(t)}\mathcal{M}_{r}$ is the tangent space of $\mathcal{M}_{r}$ at the current state $Y(t)$ and $\lVert\cdot\rVert_{F}$ denotes the Frobenius norm. The minimization condition leads to a differential equation for $Y$ on the low-rank manifold, which has to be solved numerically. Standard integrators fail due to the presence of small singular values. Therefore, a particular integrator, the so-called projector-splitting integrator was developed in [22]. Its convergence for non-stiff $G$ and its robustness with respect to the small singular values was proven in [20].

Condition (8) states that $\dot{Y}(t)$ is obtained by an orthogonal projection of $G(t,Y(t))$ onto the tangent space $\mathcal{T}_{Y(t)}\mathcal{M}_{r}$, i.e.

[TABLE]

where $P(Y)$ is this orthogonal projection. The matrix $Y(t)$ is represented as $U(t)S(t)V(t)^{\mathsf{T}}$, where $U,V\in\mathbb{R}^{d\times r}$ with orthonormal columns and $S\in\mathbb{R}^{r\times r}$. Note that $S$ is not assumed to be diagonal. Using the explicit form of the projector, which was computed in [21], the evolution equation for the approximation $Y$ is seen to be

[TABLE]

The projector-splitting integrator itself makes use of splitting methods. The right-hand side of (9) is split up into three subproblems. The integrator only deals with differential equations for the low-rank factors $U,S$ and $V$. The practical algorithm for a general nonlinearity $G(t,Y)$ is given in [20, Sect. 2.2]. Starting from $U_{0}S_{0}V_{0}^{\mathsf{T}}$, a rank-$r$ approximation to the initial data $N(t_{0})$, one step of the integrator is as follows.

1. a.

Solve $\dot{K}(t)=G(t,K(t)V_{0}^{\mathsf{T}})V_{0}$ with initial value $K(t_{0})=U_{0}S_{0}$. Then orthonormalize $K(t_{1})$ by a QR decomposition to get $K(t_{1})=U_{1}\widehat{S}_{1}$, where $U_{1}\in\mathbb{R}^{d\times r}$ has orthonormal columns and $\widehat{S}_{1}\in\mathbb{R}^{r\times r}$. 2. b.

Solve $\dot{S}(t)=-U_{1}^{\mathsf{T}}G(t,U_{1}S(t)V_{0}^{\mathsf{T}})V_{0}$ with initial value $S(t_{0})=\widehat{S}_{1}$. Set $\widetilde{S}_{0}=S(t_{1})$. 3. c.

Solve $\dot{L}(t)=G(t,U_{1}L(t)^{\mathsf{T}})^{\mathsf{T}}U_{1}$ with initial value $L(t_{0})=V_{0}\widetilde{S}_{0}^{\mathsf{T}}$. Then orthonormalize $L(t_{1})$ by a QR decomposition to get $L(t_{1})=V_{1}S_{1}^{\mathsf{T}}$, where $V_{1}\in\mathbb{R}^{d\times r}$ has orthonormal columns and $S_{1}\in\mathbb{R}^{r\times r}$.

The above scheme can be tailored to the particular form of $G$. More details are given in the sections below. A detailed convergence analysis of the proposed integrator will be presented elsewhere.

3 Differential Lyapunov equations

Differential Lyapunov equations arise for $G(t,X)=Q$ with $Q$ being a constant matrix. Thus, we consider the following initial value problem

[TABLE]

where $X(t)$, $A$, $Q$ $\in\mathbb{R}^{d\times d}$, and $t\in[t_{0},T]$. The matrix $A$ typically arises from the discretization of a differential operator. Further, $Q$ and the initial data $X_{0}$ are symmetric and positive semidefinite. Since (10) is a linear differential equation with constant coefficients the solution exists for all times. Moreover, the solution is also symmetric and positive semidefinite. This is a straightforward consequence of (2), see also [1].

3.1 A low-rank split-step integrator

As explained in Section 2 we split (10) into the following two subproblems:

[TABLE]

where $t\in[t_{0},T]$. Note that the exact solution of (11b) simply is $N(t)=N_{0}+~{}tQ$. Employing this representation directly in our splitting method, however, would require one additional SVD per time step. To avoid this computational overhead, we follow the approach chosen in [20] and use the projector-splitting method.

Since the solution of (10) is symmetric and positive semidefinite we represent it as a rank-$r$ matrix of the form $USU^{\mathsf{T}}$, where $U\in\mathbb{R}^{d\times r}$ has orthonormal columns and $S\in\mathbb{R}^{r\times r}$ is symmetric. This is an SVD-like decomposition with the additional constraint of symmetry.

The linear problem (11a) can be treated as explained in Section 2. The projector-splitting integrator for the solution of (11b) in its standard formulation, however, does not preserve the symmetry and positive semidefiniteness of the solution. Let $U_{0}S_{0}U_{0}^{\mathsf{T}}$ be a low-rank decomposition of the initial data $N_{0}$. In order to preserve the symmetry of the problem we propose a modification of the algorithm given at the end of Section 2, see also [20, Sect. 2.2]. The modified projector-splitting integrator is described in step 5 of Algorithm 1. In the first substep we update the value of the left-sided orthonormal factor $U_{0}$ to $U_{1}$, whereas we keep the right-sided one till the last step. There we get $L(t_{1})=S_{1}\widetilde{U}_{1}^{\mathsf{T}}$. Imposing $\widetilde{U}_{1}=U_{1}$ and taking advantage of the orthonormality of $U_{1}$ we obtain $S_{1}=L(t_{1})U_{1}$. The low-rank approximation $N_{1}$ at $t_{0}+\tau$ is then simply $U_{1}S_{1}U_{1}^{\mathsf{T}}$.

In Algorithm 1 we summarize the above proposed Lie splitting for the solution of (10). Note that the flows of all these substeps can be computed exactly. Moreover, the whole time integration is performed by working only with the low-rank factors of the solution, which has several computational advantages such as less computing time and less memory requirements. By using the explicit structure of the subproblems, it is easy to show that $S_{1}$ and consequently $N_{1}$ are symmetric and positive semidefinite.

Lemma 1**.**

The matrix $S_{n}$ arising in step 5c of Algorithm 1 is symmetric and positive semidefinite.

Proof.

Without loss of generality we prove the lemma for $n=0$. For simplicity, we will also omit the upper index $A$ and write $U_{0}$ instead of $U_{0}^{A}$, etc. Starting from the rank-$r$ approximation $U_{0}S_{0}U_{0}^{\mathsf{T}}$ of $N_{0}$, the first step of the projector-splitting integrator gives

[TABLE]

The solution of the second substep is

[TABLE]

From step 5c we get

[TABLE]

Using the orthonormality of $U_{1}$, we compute $\widehat{S}_{1}=U_{1}^{\mathsf{T}}U_{0}S_{0}+\tau U_{1}^{\mathsf{T}}QU_{0}$ from (12). Inserting this formula into (13) we get the symmetric formula

[TABLE]

which proves the assertion. ∎

3.2 Numerical results for a parabolic problem

In order to illustrate the behavior of the integrator, we take a simple test problem. We consider the heat equation

[TABLE]

on $\Omega=(0,1)^{2}$ with homogeneous Dirichlet boundary conditions. The associated DLE is of the form (10), where $A$ arises from the spatial discretization of the Laplacian. Thus, we have $A=\widetilde{A}\otimes I+I\otimes\widetilde{A}$, where $\widetilde{A}$ is the 1D standard centered finite differences matrix with $\widetilde{d}$ uniformly spaced grid points in the interior. In the following we take $\widetilde{d}=20$ and obtain $d=\widetilde{d}^{2}=400$. The matrix $Q$ in our example has random coefficients and is of rank 5. The initial value $X_{0}$ is chosen to be a random matrix of rank 10. This choice ensures the low-rank behaviour of the solution of the DLE. This can be observed in Figure 1, left where we plot the first 90 singular values of a numerical solution computed with DOPRI5 [17] at the final integration time $T=0.1$. This solution is taken as the reference solution for all tests in this section. The code DOPRI5 is an explicit Runge–Kutta method of order 5 with adaptive step size strategy.

In Figure 1, right we show the error behaviour of the Lie splitting (5) as described in Algorithm 1. The errors are measured in an appropriately scaled Frobenius norm

[TABLE]

We observe that the error of our method is composed by two different contributions, the error due to the outer splitting into (11a) and (11b) and the error due to the low-rank approximation. As long as the error due to the low-rank approximation is not dominant, we observe the expected order of convergence one for Lie splitting. On the other hand, if the low-rank approximation is poor, decreasing the step size will not improve the quality of the solution. We observe a stagnation of the error around certain values depending on the rank. These errors can be related with the magnitude of the first singular value discarded for each choice of the approximation rank, as can be observed by comparing Figure 1, left and right. In Figure 2 we show the corresponding results for Strang splitting for $d=400$ and $d=3600$, respectively. On the left we observe the expected order of convergence two for the outer splitting. Again, when the approximation rank is too low, the outer error becomes independent of any step size refinement. An order reduction shows up in the figure on the right. It is due to the fact that the inhomogeneity $Q$ is not in the domain of the Laplacian. Therefore, $AQ$ cannot be bounded independently of $d$, see [13].

In Table 1 we list the defects in symmetry

[TABLE]

where $Y$ is the solution obtained with Lie splitting as defined in Algorithm 1 and $Y_{\text{ref}}$ is the reference solution. We observe that the symmetry is numerically preserved, as expected from Lemma 1. In Table 2 we show that the standard non-symmetric integrator as presented in Section 2 is not symmetry-preserving. It is therefore necessary to modify the projector-splitting as explained in detail in Algorithm 1.

Due to Lemma 1, our integrator also preserves positive semidefiniteness of the solution. To show this, we denote with $\widehat{Y}$ the nearest symmetric positive semidefinite matrix to $Y$ obtained by the method proposed in [19]. In Table 3 we list the defects in positive semidefiniteness

[TABLE]

for different approximation ranks and step sizes.

3.3 Comparison with a standard approach

In order to give more insight into the performance of the low-rank splitting described in Algorithm 1 we carry out a comparison with a standard method. The basic idea is to discretize the DLE in time to convert the differential problem into an algebraic one. Since our aim is to compare the first-order low-rank splitting, we apply here the backward Euler method with time step $\tau$ to obtain the numerical approximation $X_{n+1}$ to $X(t_{n+1})$. This gives

[TABLE]

After recombining the terms, we end up with the following algebraic Lyapunov equation (ALE)

[TABLE]

Several approaches are available in the literature for solving this type of equations. We employ the one proposed in [31]. It consists of projecting the ALE onto an approximation space, generated as combination of Krylov subspaces in $A$ and $A^{-1}$. Such an approach is also called extended Krylov method. The reduced equation is then solved by means of a direct solver.

In particular, consider an ALE of the form

[TABLE]

with $\widetilde{A}\in\mathbb{R}^{d\times d}$ and $\widetilde{B}\in\mathbb{R}^{d\times s}$, $s\ll d$. Let us denote with $\widetilde{V}$ the matrix whose columns span the approximation space and with $Y$ the solution of the projected ALE. Then the solution of the original ALE is reconstructed as

[TABLE]

A low-rank solution is generated discarding all the eigenvalues of $Y$ which are smaller than a certain tolerance tolY. Then the solution $X$ is given in low-rank form as $ZZ^{\mathsf{T}}$. A second user-supplied parameter, denoted by tol, is used for the stopping criterion of the algorithm. The iteration is stopped if

[TABLE]

where the superscript $m$ denotes the number of the iteration. The algorithm associated with these procedure is called K-PIK (Krylov-plus-inverted Krylov). For further details we refer to [31].

In the following we report some numerical experiments for the example we considered in the beginning of Section 3.2. We solve it until $T=0.1$. An algebraic equation of the form (17) has to be solved for each time step. In order to use the K-PIK procedure we set

[TABLE]

where $Z_{n}$ is the low-rank factor of the solution $X_{n}$ at time $t_{n}$, i.e., $X_{n}=Z_{n}Z_{n}^{\mathsf{T}}$ and $R$ is the low-rank factor of $Q$, i.e., $Q=RR^{\mathsf{T}}$.

In Figure 3 we illustrate the behaviour of this integrator for different values of tolY. The results are all obtained in Matlab with $\texttt{tol}=10^{-12}$. The projected ALE is solved with Matlab’s lyap function. Moreover, as a comparison, we compute the solution of (10) by directly solving (17) by the routine lyap. Since in both the cases the problem is discretized in time in the same way, this comparison reduces to a comparison between the K-PIK procedure and the lyap method. In Figure 3, left we show the error behaviour of the two methods. The K-PIK procedure is tested for different values of tolY. The error we observe is basically just the first order error of the backward Euler method. Only when the time step becomes very small ($\tau\approx 10^{-6}$) we observe an additional error if tolY is not chosen sufficiently small. This is basically the error due to the neglect of some eigenvalues of the solution. In Figure 3, right we observe how the choice of tolY influences the rank of the solution. The full rank given by the solution computed with lyap is recovered by the K-PIK procedure with tolY= $10^{-12}$.

In Figure 4 we compare the results obtained with the Lie splitting approach described in Algorithm 1 with the ones obtained from the combination of the backward Euler method with the K-PIK procedure. In Figure 4, left the relative errors for the two methods are shown. We observe that the first-order splitting method is more accurate than the backward Euler method for this example. This is due to the fact that we solve the linear subproblem exactly.

In Figure 4, right the computing time at $T=0.1$ is displayed. The solution of the ALE (17) is computed by the K-PIK procedure, where the Choleski decomposition of the matrix $\widetilde{A}$ is computed once and for all. The computation of the matrix exponential required for the Lie splitting in Algorithm 1 was carried out by Leja interpolation [7] with single precision. We conclude that the splitting proposed here is considerably faster than the backward Euler method combined with the K-PIK approach.

3.4 Numerical results for the simulation of El Niño

The quasiperiodic weather phenomenon El Niño which is characterized by an unusual warming of the sea surface in the Indo-Pacific ocean, has a huge impact on the climate worldwide and is also responsible for many natural disasters. Different models (both deterministic and stochastic) are used to describe the variation in the sea surface temperature. We consider in the following a stochastic advection equation driven by additive noise

[TABLE]

where the vector $X$ contains the sea surface temperature (SST) anomalies, the operator $\mathcal{A}=u\cdot\nabla$ describes the advection by the ocean currents $u$, and $F$ is a random vector which aggregates external forces like wind stress or evaporation [8, 26]. Following [26] we construct the SST anomalies, using the dataset OISST [24], whereas for the ocean currents we use the dataset OSCAR [14], both from the National Oceanic and Atmospheric Administration (NOAA). We model the random force $F$ by a Gaussian white noise process with zero mean and covariance operator $\mathcal{Q}$ and compute the first two moments by

[TABLE]

where $t\in[t_{0},T]$ and we denote by $\widehat{X}$ the expectation and by $\mathcal{P}$ the covariance of $X$ [23]. As the solution $X$ is a Gaussian random field, it is completely defined by its second-order statistics.

We discretize the operator $\mathcal{A}$ via centered finite differences and obtain the matrix $A$. The first subproblem (20) is solved by a method based on Leja interpolation. Discretizing the Indo-Pacific ocean we get $3900$ grid points, hence solving the large-scale differential Lyapunov equation (21) with full rank is rather expensive. For such computations, our splitting algorithm has substantial advantages. It requires considerably less computing time and memory. In order to illustrate its accuracy, we make a comparison of the proposed algorithm with a second order scheme. Similarly as in Section 3.3 we apply the backward Euler method to the DLE (21) and solve the resulting ALE by the K-PIK procedure. In order to obtain a second order scheme we then use Richardson extrapolation. Given $Z_{n}$ which is a low-rank factor of $P_{n}$, we make one step with step size $\tau$ and get an approximation $\tilde{Z}_{n+1}$ and two steps with step size $\frac{\tau}{2}$ to obtain $\hat{Z}_{n+1}$. The numerical approximation $Z_{n+1}$ is then given by

[TABLE]

As discussed in Section 3.3 one has to choose parameters tol and tolY, which can have an impact in the solution. In the following we use $\texttt{tol}=10^{-10}$ and $\texttt{tolY}=10^{-12}$. Figure 5 shows the relative error of this approach. The reference solution is obtained by the same method but with step size $2^{-8}$, which is 8 times smaller than the smallest step size used in Figure 5.

In Table 4 we compare the error in Frobenius norm (14) of the low-rank integrator with respect to the second order method (22) for the years $2013$ (no event) and $2015$ (a strong event), respectively. In both experiments, we took as starting time June 15. In the table we display the absolute errors of the low-rank approximation after one, two and three weeks for two different approximation ranks. We observe that the error gets smaller for rank 100. For higher ranks, however, the splitting error dominates in this example.

We further make a realization of the stochastic random field. In Figures 6 and 7 the SST anomalies are given in degree Celsius. We show a section of the Indo-Pacific ocean with Australia being on the lower left part and America on the right part. The black part indicates the land, bright colors indicate higher temperatures and dark ones lower temperatures than usual. In Figure 6 one can see that the temperature is nearly uniformly distributed, whereas we observe from Figure 7 the typical face of an El Niño event with the unusual warming of the sea surface.

4 Differential Riccati equations

We now consider the case where $G(t,X)=Q-XPX$ with constant matrices $Q$ and $P$. The resulting differential equation

[TABLE]

is called differential Riccati equation. Again $X(t)$, $A$, $Q$, $P$ $\in\mathbb{R}^{d\times d}$, and $t\in[t_{0},T]$. The matrices $Q$ and $P$, and the initial value $X_{0}$ are symmetric and positive semidefinite. The global existence and positive semidefiniteness of the solution is guaranteed under these conditions, see [12].

4.1 A low-rank split-step integrator

Following the general idea of the paper, we split (23) into the following two subproblems:

[TABLE]

where $t\in[t_{0},T]$. The linear problem (24a) is handled as in Section 3.1. The projector-splitting integrator for the solution of (24b) can be specified taking into account the form of the non-linearity. Moreover, we propose the same kind of modification as for DLEs in order to get a low-rank decomposition of the form $USU^{\mathsf{T}}$. The resulting algorithm for the Lie splitting is given in Algorithm 2.

4.2 Numerical results for an optimal control problem

DREs arise, e.g., in optimal control for linear quadratic regulator problems with finite time horizon $T$ for parabolic partial differential equations. Thus we consider the linear control system

[TABLE]

where $A\in\mathbb{R}^{d\times d}$ and $B\in\mathbb{R}^{d\times m}$ are the system matrices, $x\in\mathbb{R}^{d}$ are the state variables and $u\in\mathbb{R}^{m}$ is the control. The output $y\in\mathbb{R}^{q}$ is defined as $y=Cx$, where $C\in\mathbb{R}^{q\times d}$. Both $m$ and $q$ are much smaller than the number $d$ of degrees of freedom. The functional that has to be minimized is

[TABLE]

where $Q\in\mathbb{R}^{q\times q}$ is symmetric and positive semidefinite, and $R\in\mathbb{R}^{m\times m}$ is symmetric and positive definite. The optimal control is given in feedback form by $u_{\text{opt}}(t)=-R^{-1}B^{\mathsf{T}}X(t)x(t)$, where $X(t)$ is the solution of the following DRE

[TABLE]

Note that the solution of (25) converges for $T\rightarrow\infty$ to a steady state. This limit is given as the solution of the algebraic equation

[TABLE]

In order to illustrate the behaviour of Algorithm 2, we consider a test example proposed in [27]. We consider the following diffusion-advection equation

[TABLE]

on $\Omega=(0,1)^{2}$ with homogeneous Dirichlet boundary conditions. The matrix $A$ arises from the spatial discretization of (27) using standard central finite differences, with $\widetilde{d}$ uniformly spaced grid points in each dimension. We denote the discretization points in the interior of $\Omega$ in $x$ direction with $x_{i}=i\delta$, $\delta=(\widetilde{d}+1)^{-1}$ for $i=1,\dots,\widetilde{d}$. We take $B\in\mathbb{R}^{d\times 1}$ with

[TABLE]

and $C\in\mathbb{R}^{1\times d}$ with

[TABLE]

Further, we choose $R=I$ and $Q=100I$.

In Figure 8 we show the behaviour of a reference solution of (25) computed with DOPRI5. The left figure shows the convergence of the solution of (25) to the solution of (26) for two different initial values $X_{0}=0$ and $X_{0}=I$, respectively. As expected this limit is independent of the choice of the initial data. The relative difference between (25) and (26) in Frobenius norm is shown as function of time. Further, in Figure 8, right we display the rank of the reference solution as a function of time.

In the following numerical example we take the initial value $X_{0}=0$, the final time $T=0.1$ and the above reference solution computed with high accuracy. We show the results for $d=400$. The equations in step 5 of Algorithm 2 are quadratic matrix differential equations. We solve them by means of the classical explicit Runge–Kutta method of order 4 with the same step size.

In Figure 9, left we show the error behaviour of the Lie splitting given in Algorithm 2. As pointed out above for DLEs, we observe that the error is composed by two different contributions. The choice of a small approximation rank results in stagnation of the error around certain rank-dependent values. On the other hand, if the approximation error becomes small enough, one observes the usual order of convergence one for the outer Lie splitting. In Figure 10, left we present the corresponding results for Strang splitting. It shows the expected order of convergence two for sufficiently high rank.

In Figures 9, right and 10, right we show the defects in symmetry (15) and positive semidefiniteness (16) in continuous and dotted lines, respectively. Although our method does not preserve these two features of the solutions, it is remarkable that the defects are always negligible with respect to the overall error of the method.

4.3 Generalized differential Riccati equations

In a stochastic version of the linear quadratic regulator problem, a generalised differential Riccati equation (GDRE) arises

[TABLE]

We point out that (28) has the same structure as the DRE (23) except for the term $C^{\mathsf{T}}XC$. The matrix $C$ is a weighting term in the noise perturbation, see [33]. Recently, a numerical method for GRDEs of the form (28) was proposed in [10].

The above GDRE shares the structure of (1) with $G=Q+C^{\mathsf{T}}XC-XBR^{-1}B^{\mathsf{T}}X$. It is therefore possible to employ the approach explained in Section 2. Two subproblems of the form (3a) and (3b) arise and they can be solved in low-rank form as presented in this work. The projector-splitting integrator can be adapted to the form of the nonlinearity, however, a modification to preserve symmetry and positive semidefiniteness might be required.

5 Conclusions

We proposed a new low-rank integrator for a class of matrix differential equations which includes, among others, differential Lyapunov and differential Riccati equations. A low-rank approximation of the solution is computed in a dynamical way, working only with the low-rank factors of the solution. This approach gives substantial advantages in terms of computing time and memory requirements. Moreover, the integrator can handle stiffness in an efficient way. Splitting methods form the core ingredient of the new method. They make it possible to treat the stiff part of the equation separately from the non-stiff one. Numerical results for differential Lyapunov and differential Riccati equations are discussed, and a simulation of the weather phenomenon El Niño is presented.

Acknowledgements

This work was supported by the Austrian Science Fund (FWF) – project id: P27926.