Dynamical low-rank integrators for second-order matrix differential equations

In this paper, we construct and analyze a new dynamical low-rank integrator for second-order matrix differential equations. The method is based on a combination of the projector-splitting integrator introduced in Lubich and Oseledets (BIT 54(1):171–188, 2014. https://doi.org/10.1007/s10543-013-0454-0) and a Strang splitting. We also present a variant of the new integrator which is tailored to semilinear second-order problems.

for the factors of a low-rank decomposition resembling a singular value decomposition. Compared to the approximation of the full matrix solution, working only with the factors of the low-rank decomposition significantly reduces the computational costs and the required storage. Unfortunately, the integrator of [10] suffers from illconditioning in the presence of small singular values, a situation which is called over-approximation, i.e., the rank chosen within the method exceeds the rank of the actual solution. A projector-splitting integrator which is robust in the case of over-approximation was introduced in [12] by Lubich and Oseledets. It is based on a clever Lie-Trotter splitting of the projected right-hand side, by which the subproblems can be solved exactly. Another dynamical low-rank integrator for first-order matrix differential equations is the unconventional robust integrator recently presented in [2]. It is especially suited for strongly dissipative problems. Although its derivation is rather different from the construction of the projector-splitting integrator, both methods share their robustness in the case of over-approximation and they satisfy the same error bounds.
Both the projector-splitting integrator and the unconventional integrator have been applied to a variety of first-order matrix differential equations, e.g., for Schrödinger equations in [2], for the Vlasov-Poison equation in [4], for Vlasov-Maxwell equations in [5], and for Burgers' equation with uncertainty in [11]. However, to the best of our knowledge, for second-order matrix differential equations of the form with large m, n, such integrators have not been considered so far. The obvious technique of reformulating (1) into a first-order system and applying the projector-splitting integrator [12] for first-order matrix differential equations behaved poorly in our numerical experiments. The reason might be that the inherent structure of the second-order problem is ignored by this procedure, causing the approximation quality to deteriorate. Therefore, we propose to combine the projector-splitting integrator with a Strang splitting. The resulting method is closely related to the leapfrog scheme and called LRLF method in the following. It is a robust and reliable dynamical low-rank integrator for second-order equations of type (1), provided that the exact solution A(t) and its derivative A (t) can be well approximated by matrices of low rank. When the exact flows used within the Strang splitting preserve the low rank of the previous approximation, our integrator reduces to the leapfrog scheme if the rank chosen in the method is sufficiently large. We also develop a variant of the scheme which is tailored to semilinear second-order equations. For this, we combine our newly developed scheme with the ideas in [15], where a dynamical low-rank integrator for stiff semilinear firstorder equations was derived based on the projector-splitting integrator.
For the projector-splitting integrator, a detailed error analysis was provided in [9]. It relies on an exactness property of the integrator, namely that it provides the exact solution if this solution preserves the (low) rank of the initial value for all times and the exact initial value is used to start the integrator. Unfortunately, this is no longer true for the LRLF scheme. Nevertheless, we will provide error bounds under similar assumptions as in [9]. The paper is organized as follows: In Sect. 2, we briefly recall the projector-splitting integrator introduced in [12] by Lubich and Oseledets. The construction of the LRLF scheme is presented in Sect. 3 and its error analysis in Sect. 4. A modification of the LRLF scheme which is tailored to semilinear second-order problems is developed in Sect. 5.
Throughout this paper, m, n, and r are natural numbers, where w.l.o.g. m ≥ n r . If n > m, we consider the equivalent differential equation for the transpose. By M r we denote the manifold of complex m × n matrices with rank r , The Stiefel manifold of m × r unitary matrices is denoted by where I r is the identity matrix of dimension r and U H is the conjugate transpose of U .
The singular value decomposition of a matrix Y ∈ C m×n is given by where σ 1 ≥ . . . ≥ σ n ≥ 0 are its singular values. It is well known that for r < n, the rank-r best-approximation to Y w.r.t. the Frobenius norm · is given by where Σ = diag(σ 1 , . . . , σ r , 0, . . . , 0) and For a given step size τ we use the notation t k = kτ for any k with 2k ∈ N 0 .

The projector-splitting integrator
In this section, we briefly review dynamical low-rank approximations as introduced in [10,12]. We start with the following problem: Given some time-dependent matrix A(t), find a low-rank approximation A(t) ∈ M r such that In [10], this is done by imposing that A (t), which is contained in the tangent space For an initial value A(0) = A 0 ∈ M r , condition (2) is equivalent to a Galerkin condition. In fact, then A solves the evolution equation where P A(t) denotes the orthogonal projection onto the tangent space T A(t) M r . A natural choice for the initial value A 0 of (3) is the rank-r best approximation to A(0). For all A ∈ M r there is a non-unique low-rank factorization which allows us to write the orthogonal projector P( A) onto T A(t) M r as cf., [10,Lemma 4.1]. The dynamical low-rank integrator constructed in [12] is based on a Lie-Trotter splitting, applied to the differential equation (3) with P( A) given in (4). Given a step size τ > 0, the first integration step consists of solving the three subproblems As shown in [12, Lemma 3.1], all subproblems (5) can be solved exactly on M r when the additional conditions can be written in terms of the increment The integration process is continued with initial value Y γ (τ ) ≈ A(τ ) in the next time step. A single time step of the resultant projector-splitting integrator is presented in Algorithm 1, version (α).
The above approach is also suitable for the construction of low-rank approximations to the unknown solution A(t) of the first-order differential equation As explained in [12,Section 3.4], the only change affects the replacement of the increment A in Algorithm 1 in a way resembling the explicit Euler method, The global error of this first-order scheme depends on the quality of the approximation of the exact solution A(t) of (6) by a low-rank matrix (for t ∈ [0, T ]) and on properties of the right-hand side F, cf. [9]. Alternatively, any other explicit scheme, e.g., a fourth order Runge-Kutta method can be used to replace the increment A. We denote the latter method by PSI-RK.

Remark 1
The computational complexity of Algorithm 1 is dominated by the two products A V and A H U in lines 7 and 11, respectively, the two Q R-decompositions of matrices of dimension m × r and n × r , respectively, and the three matrix-matrix products in lines 8, 10, and 11. For an efficient implementation, it is important to compute the products A V and A H U without computing A explicitly.

Dynamical low-rank approximation of second-order matrix ODEs
Next, we devise a low-rank integrator for second-order matrix differential equations of the form (1). A straightforward practice would be to rewrite (1) as a first-order system , and to apply Algorithm 1. However, numerical tests showed that the quality of the numerical solutions deteriorates over time, possibly caused by the neglection of the structure of the the right-hand side of (7). We thus use a different ansatz and first split (7) into and then apply a standard Strang splitting. Solving the subproblems exactly leads to the well-known leapfrog or Störmer-Verlet scheme, cf. [6, Section 1.5]. If approximations to A = B are not required at full time steps, then the most economic implementation of the leapfrog scheme is to combine (9a) and (9c) via is computed from (9a). A low-rank integrator for the second-order equation (1) is now designed by approximating the exact flows of the subproblems in (8) by their respective low-rank flows using the projector-splitting method [12]. First, we determine initial values A 0 and B 0 as rank-r A and rank-r B best-approximations to A 0 and B 0 , respectively. After k time steps, the low-rank approximations The low-rank matrices A k+1 ≈ A(t k+1 ) and B k+ 1 ) are obtained by approximating the solutions of (8) by applying Algorithm 1 to Since the exact solutions of (10) read the increments B k− 1 2 and A k are given explicitly as The resulting dynamical low-rank integrator for second-order matrix ODEs will be named LRLF method, for low-rank leapfrog integrator. It is presented in Algorithm 2.
There is a close relationship between the LRLF and the leapfrog schemes: (8) with rank A 0 = r A , rank B 0 = r B stay of rank r A and r B , respectively, then the solutions of the LRLF and leapfrog schemes coincide.

Theorem 2 If for all t ∈ [0, T ] the exact solutions A(t) and B(t) of the subproblems
Proof The LRLF scheme is derived by exchanging the exact flows of the leapfrog scheme by their corresponding low-rank flows. Since the exact flows preserve the rank by assumption, the application of the exactness property of the projector-splitting integrator [12, Theorem 4.1] yields the desired result.
Clearly, in general the LRLF scheme does not inherit the exactness property of the projector-splitting integrator, because of the splitting error. A detailed error analysis of the LRLF iteration will be presented in Sect. 4.
In the same way, a variant of the LRLF method for the second-order matrix differential equation or equivalently (1), LRLF scheme, single time step

Algorithm 2 DLR integrator for second-order ODEs
Algorithm 3 DLR integrator for second-order ODEs (12), LRLF(ω) scheme, single time step is based on the non-staggered version of the leapfrog scheme (9), procuring approximations to A on the same time-grid as approximations to A. Here, the function F has to be evaluated twice in each time step so that the computational effort is larger. The resulting method is called the LRLF(ω) scheme. It is described in Algorithm 3.

Remark 3
A dynamical low-rank integrator for second-order matrix differential equations (1) can also be derived based on the two-step formulation of the leapfrog scheme, cf. [16,Section 4.2]. In numerical experiments it was observed that the resulting scheme does not yield better approximations than the LRLF method, and even produces larger errors in some examples.

Error analysis of LRLF
In the following, we analyze the error of the LRLF scheme given in Algorithm 2 when applied to (1) with a right-hand side F which is Lipschitz-continuous with a moderate Lipschitz constant L, i.e., F satisfies Our analysis relies on the error analysis in [9] for the projector-splitting integrator.
Recall that the LRLF scheme is derived from the leapfrog scheme (9), which is stable under the CFL condition, cf. [8], We therefore assume that the step size τ always satisfies (14).

Assumption 1 The exact solution
Additionally, there exist sufficiently large constants γ A and γ B such that (15) is also satisfied for all For the approximations A k ≈ A(t k ) and B k+ given in (11), we define By the triangle inequality, we have The analysis of the LRLF scheme is organized in two lemmas and a theorem. Our The constants C LF A and C LF B are given explicitly as Proof By Taylor series expansion, we have Hence for k = 0, we have by (20c), (11c), and (13). Employing (17) shows (18a). Using (18a) and (20a) yields Together with (17) this proves (18b).
In [9, Section 2.6.1] it was shown that the error between a time-dependent matrix A(t) satisfying (15a) and the rank-r A approximation Y 1 ≈ A(τ ) computed by Algorithm 1 started from X A (0) ∈ M r A is bounded by In the next lemma we eliminate E cf. (25). The estimate on the global error then follows from (17).
We are now able to prove a global error bound.
respectively, as long as t k+4 ≤ T .

Proof The bound for E k+1
A is a direct consequence of (31) with j = k + 1. The is obtained as for E k+1 A in the proof of Lemma 5, but starting from substituting (19b) into (19a).
The error of the LRLF scheme is hence a combination of two error contributions: an error caused by the low-rank approximations, and a time discretization error stemming from the leapfrog scheme. If the low-rank errors ρ A , ρ B , ρ A , and ρ B are small, i.e., the solutions A, B of (1) are well-approximated by low-rank matrices, the time discretization error dominates.
In the proofs of Lemma 5 and Theorem 6 we used the following Gronwall-type lemma.

Lemma 7
Let τ, L ≥ 0 and {M k } k≥0 a nonnegative, monotonically increasing sequence. If the nonnegative sequence {E k } k≥0 satisfies Proof Define ε k :=E k /M k for all k ≥ 0. The sequence {ε k } k≥0 is nonnegative and satisfies due to the monotonicity of {M k } k≥0 . The statement then follows from [1, Lemma 3.8].

Dynamical low-rank integrator for semilinear second-order matrix differential equations
We now consider semilinear second-order matrix differential equations of the form with given Hermitian, positive semidefinite matrices Ω 1 ∈ C m×m and Ω 2 ∈ C n×n of large norm and a Lipschitz continuous function f with a moderate Lipschitz constant. Our aim is to construct a dynamical low-rank integrator which exploits the semilinear structure of the right-hand side of (35) and yields smaller errors than the LRLF scheme. A dynamical low-rank integrator for stiff semilinear first-order equations was proposed in [15]. It is based on the important property of linear first-order equations, that the exact solution of an initial value problem stays in M r for all times if the initial value is in M r . This is in general not the case for linear second-order equations, so that additional considerations are required here.
We transform (35) into an equivalent first-order problem and split the right-hand side into two linear and one nonlinear part by introducing weights ω i ≥ 0, i = 1, 2, 3, A natural choice is ω 2 i = 1/3, i = 1, 2, 3, but different weightings are feasible. The split equations can be written as The solution of the linear problems (36a) and (36b) can be expressed in terms of the matrix exponential The full splitting scheme reads where φ Ω 1 τ 2 and φ Ω 2 τ 2 denote the numerical flows given by Algorithm 1 with step size τ 2 to the exact solutions of (36a) and (36b), respectively. φ S τ denotes the numerical flow of the LRLF(ω) scheme, described in Algorithm 3, with right-hand side f . The overall method (37) is called LRLF-semi.
When approximations at full time steps are dispensable, the last half step φ Ω 1 τ 2 in (37) can be combined with the first one of the next time step.

Homogeneous wave equation
In our first example we set Ω = [−π, π] 2 and γ = 0. With the initial values where k x , k y ∈ R, (38) models a planar wave traveling into the direction of k = (k x , k y ). For the discretization in space we use second-order finite differences with n and m discretization points in x-and y-direction, respectively. This yields the second-order matrix differential equation (35)  For m = n = 512, T = 10 and 2k x = k y = 2, we compute low-rank approximations to the exact solution of (35) with the LRLF and the LRLF-semi schemes and ranks r A = r B = 2. For the latter method, we use the weight configurations and ω 2 1 = sin 2 ϕ, ω 2 2 = cos 2 ϕ, ϕ = arctan k y k x .
While the first choice treats the two subproblems equally, the second takes the direction of motion into account. Moreover, we compute low-rank approximations with the projector-splitting integrator (Algorithm 1), applied to the equivalent first-order system (7) of (35) and rank r = 4. For the inner integration we use the explicit Euler method (PSI) and the classical Runge-Kutta method of order 4 (PSI-RK). In Fig. 1, the relative global error between the exact solution A and its low-rank approximations A at time t = T is displayed. Temporal convergence of order 2 is observed for both the LRLF and the LRLF-semi scheme. Compared to the LRLF method, the LRLF-semi scheme is more accurate and allows slightly larger step sizes. Choosing the weights ω i in consideration of the direction of motion improves the approximation substantially. In contrast, the PSI and PSI-RK methods require significantly smaller step sizes τ and come with larger errors than the schemes for second-order matrix differential equations. with l 0 = π/30 and w 0 = π/3. Due to the strong localization of the initial value, we discretize in space by a pseudospectral method in y-direction and fourth-order finite differences in x-direction. This yields the second-order matrix differential equation (35) with the cubic nonlinearity

Semilinear wave equation
where the symbol • denotes the entrywise matrix product. The matrix Ω 2 1 ∈ R m×m is given by where F m denotes the discrete Fourier transformation operator with m modes, and Ω 2 2 ∈ R n×n denotes the positive semidefinite, symmetric Toeplitz matrix with first row For m = 4096, n = 512, T = 0.5π , and r A = r B = 10, we compute low-rank approximations to the solution of (35) with the LRLF and the LRLF-semi scheme. A reference solution has been computed with the leapfrog scheme and step size τ 0 = 10 −4 T . For the LRLF-semi scheme, we again use two weight configurations, where the first treats all subproblems equally and the second takes the direction of propagation into account, and ω 2 1 = 2 3 , ω 2 = 0, ω 2 3 = 1 3 .
In Fig. 2, the relative global error of the low-rank approximations compared to the reference solution at time t = T is displayed. Both LRLF and LRLF-semi schemes show convergence of order 2. For both weight configurations the LRLF-semi scheme yields smaller errors, and the approximation quality improves if the weights are chosen according to the direction of motion. For all integrators, the error reaches a plateau for small τ . This is caused by the low rank best approximation error to the exact solution.
In fact, the rank-r A best-approximation A best r A to the exact solution A ∈ C m×n satisfies which is a lower bound on the relative error of an arbitrary low-rank approximation. The error plateau observed in Fig. 2 suggests to choose the rank adaptively such that the temporal convergence order is preserved for smaller step sizes. This issue has been addressed in [7].

Conclusion and outlook
In the present paper, we developed and analyzed dynamical low-rank integrators for second-order matrix differential equations of the forms (1) or (35). For LRLF we proved second-order convergence under reasonable assumptions. The theoretical results have been confirmed by numerical experiments.
Both the projector-splitting integrator and the unconventional robust integrator have been successfully adapted to first-order tensor differential equations, cf. [3,13,14] and references therein. We are confident that the constructed dynamical low-rank integrators for second-order matrix differential equations can be adapted to the tensor case as well. This is part of ongoing research and will be reported in the future.