Rank-adaptive dynamical low-rank integrators for first-order and second-order matrix differential equations

Dynamical low-rank integrators for matrix differential equations recently attracted a lot of attention and have proven to be very efficient in various applications. In this paper, we propose a novel strategy for choosing the rank of the projector-splitting integrator of Lubich and Oseledets adaptively. It is based on a combination of error estimators for the local time-discretization error and for the low-rank error with the aim to balance both. This ensures that the convergence of the underlying time integrator is preserved. The adaptive algorithm works for projector-splitting integrator methods for first-order matrix differential equations and also for dynamical low-rank integrators for second-order equations, which use the projector-splitting integrator method in its substeps. Numerical experiments illustrate the performance of the new integrators.


Introduction
Dynamical low-rank integrators [17,19] have been designed for the approximation of large, time-dependent matrices which are solutions to first-order matrix differential equations whose solutions can be well approximated by low-rank matrices. The projector-splitting integrator introduced in [19] has particularly favorable properties. It is robust in the presence of small singular values, which appear in the case of over-approximation, i.e., when the approximation rank is chosen larger than the rank of the solution A(t) of (1). A variant of this method adapted to strongly dissipative problems was presented in [5]. Another variant for stiff first-order matrix differential equations was introduced in [24].
Recently, a novel dynamical low-rank integrator for second-order matrix differential equations has been constructed in [12]. Is is based on a Strang splitting of the equivalent firstorder formulation and the projector-splitting integrator [19]. This integrator is also robust in the case of over-approximation and shows second-order convergence in time. In particular, with a few amendments it is an effective method for semilinear second-order matrix differential equations, see Sect. 2.3. In applications, rank-adaptivity turns out to be essential for the efficiency of the algorithms. For first-order equations, in [4] a rank-adaptive variant of the unconventional integrator [5] was proposed. However, the approach from [4] is not applicable to the projector-splitting integrator [19]. In [8], rank-adaptivity for tensor methods for high-dimensional PDEs was based on a functional tensor train series expansion. For the special case of finite-dimensional parametrized Hamiltonian systems modeling non-dissipative phenomena, a rank-adaptive structure-preserving reduced basis method was introduced in [11]. Very recently, in [10] a Predictor-Corrector strategy for adaptivity was proposed, and in [31] the authors developed a rank-adaptive splitting method for the extended Fisher-Kolmogorov equation.
In the present paper, we discuss a general strategy for selecting the rank adaptively in the projector-splitting integrator. Increasing or decreasing the rank from one time step to the next was already proposed in [19] and quite recently in [13]. Our main contribution is a strategy for which the time step of the underlying splitting method, i.e., the Lie-Trotter splitting for first-order and the Strang splitting for second-order problems, is the only input parameter. We determine the rank such that the error of the dynamical low-rank approximation does not spoil the order of the underlying splitting method applied to the full matrix differential equation. This is achieved by propagating one additional singular value which is used for accepting or rejecting the current time step and for selecting the rank in the next step. The decision is based on an estimator of the global time-discretization error. This adaptivity control is also applicable to the dynamical low-rank integrators for stiff problems. The new dynamical low-rank integrator for (2) uses the projector-splitting integrator in the substeps of the Strang splitting, which allows to control the rank adaptively also for second-order equations. Moreover, it can be readily combined with the integrator from [5].
The paper is organized as follows: In Sect. 2, we briefly recall the projector-splitting integrator introduced in [19] and the LRLF scheme from [12]. Additionally, we sketch variants of both methods for (stiff) semilinear first-order and second-order differential equations, respectively. Section 3 is devoted to rank-adaptivity. In Sect. 4, numerical experiments illustrate the performance of the new schemes.
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 where = diag(σ 1 , . . . , σ r , 0, . . . , 0) and The Frobenius norm is denoted by · , and the Frobenius inner product by ·, · . The symbol • denotes the entrywise or Hadmard product of matrices. For a given step size τ we use the notation t k = kτ for any k with 2k ∈ N 0 .

First-order differential equations
In the dynamical low-rank approximation of the solution to first-order matrix differential equations (1), the approximation A ≈ A is determined as solution of the projected differential equation Here, P A(t) denotes the orthogonal projector onto the tangent space of the low-rank manifold M r at A(t) ∈ M r , i.e., cf. [17,Lemma 4.1], where A(t) ∈ M r is decomposed in a non-unique fashion resembling the singular value decomposition into For the initial value A 0 , typically the rank-r best-approximation to A(0) computed by a truncated SVD is used. The dynamical low-rank integrator developed in [19], also called the projector-splitting integrator, is constructed by performing a Lie-Trotter splitting on the right-hand side of (3) and solving the three subproblems on the low-rank manifold. This approach yields an efficient time-stepping algorithm for computing the desired low-rank approximations. One time-step of the projector-splitting integrator is given in Algorithm 1.

Algorithm 1
Projector-splitting integrator for low-rank approximations to the solution A(t) of (1), single time step, cf. [19, Section 3.2] S ∈ C r ×r , functions for matrix-vector multiplication with A and A H , where 4: compute Q R-decomposition V S H = L 12: 13: return U , S, V , L 14: {output: factors U , S, V of rank-r approximation A = U S V H ≈ A(t + τ ) and L = V S H (optional), 15: with U ∈ V m,r , V ∈ V n,r , S ∈ C r ×r } 16: end function

Second-order differential equations
For the second-order problem (2), a novel dynamical low-rank integrator named LRLF (low-rank leapfrog) scheme was presented in [12,Section 3]. Given r A , r B ∈ N and approximations A k ≈ A(t k ) of rank r A and B k− 1 , respectively. This integrator is based on the first-order formulation of (2), combined with a Strang splitting. The subproblems are first-order matrix differential equations, (5a, 5b) which can be solved exactly. The low-rank matrices A k+1 and B k+ 1 2 are obtained by approximating the solutions of (5a) by application of Algorithm 1 to This leads to the dynamical low-rank integrator LRLF shown in Algorithm 2.

Semilinear problems
A fixed-rank dynamical low-rank integrator for the stiff semilinear first-order problem where the norms of L 1 ∈ C m×m and L 2 ∈ C n×n are large and f is a Lipschitz continuous function with moderate Lipschitz constant, was introduced in [24]. It is based on the subproblems (7a, 7b) The solution to (7a) is given by Note that the rank of the initial value is preserved for all times [24,Section 3.2]. In contrast, the rank of the solution to the nonlinear subproblem (7b) may vary in time.
In [24], a low-rank approximation has been computed by applying a Lie-Trotter splitting to (6) and solving the subproblems (7) by the projector-splitting integrator in Algorithm 1. This method is called PSI-stiff in the following.
For semilinear second-order equations of the form with Hermitian, positive semidefinite matrices 1 ∈ C m×m , 2 ∈ C n×n and f again Lipschitz continuous, a dynamical low-rank integrator named LRLF-semi was proposed in [12,Section 5]. It is based on the equivalent first-order formulation of the second-order problem (9), where the right-hand side is split into . (10) The weights ω i ≥ 0, i = 1, 2, 3, can be chosen arbitrarily such that ω 2 1 +ω 2 2 +ω 2 3 = 1. A natural choice is ω 2 i = 1/3. The linear subproblems can be solved exactly. Low-rank approximations to these solutions are obtained by application of the projector-splitting integrator. The nonlinear subproblem is solved approximately with a variant of the LRLF scheme, cf. [12,Algorithm 3]. Denoting the numerical flows of the linear subproblems by φ 1 τ and φ 2 τ , and the numerical flow of the nonlinear subproblem as φ S τ , respectively, one step of the LRLF-semi scheme reads

Rank adaptivity
In many applications, an appropriate rank for computing a low-rank approximation to the exact solution of (2) or (1) is not known a priori and it may also vary with time. If the rank is chosen too small, the low-rank approximation lacks accuracy. Conversely, if the rank is chosen too large, the algorithm becomes inefficient.
In the following, we develop rank-adaptive variants of the projector-splitting integrators PSI and PSI-stiff for first-order problems and for the LRLF and the LRLF-semi schemes for second-order problems.

Selecting the rank
We first discuss the projector-splitting integrator for the first-order problem (1). The general idea of our rank-adaptive strategy is to approximate the exact solution of (1) by a low-rank solution of rank r k in the kth time step, but to propagate a solution of rank r k + 1. The additional information is used as an indicator whether to accept or reject the current time step, and for selecting the rank r k+1 for the next time step. Given We then compute the singular value decomposition of S k+1 , with P k+1 , Q k+1 ∈ V r k +1,r k +1 , and σ 1 ≥ · · · ≥ σ r k +1 ≥ 0 so that is the singular value decomposition of A k+1 . Given a tolerance tol, we determine r k such that by distinguishing three cases: 1. Augmentation case If σ r k +1 ≥ tol, the step is rejected and recalculated with rank r k + 2. The ranks of the initial values U k , S k , V k of the current integration step are increased by adding a zero entry to S k , This choice has been motivated by [19,Section 5.2]. The matrices U k and V k are augmented by unit vectors u ∈ C m and v ∈ C n such that Numerical tests indicate that choosing u and v as random vectors and orthonormalizing them against U k and V k is reliable and robust, but other choices are also possible. Clearly, thus the initial value of the current integration step has not changed. However, the numerical approximation has been enabled to evolve to rank r k + 2. The step is recomputed with the new initial values U * , S * , V * , and it is again checked if the new smallest singular value is sufficiently small for accepting the step. This procedure is repeated until (11) is satisfied and the step is finally accepted, see Algorithm 3 for details. 2. Reduction case If σ r k < tol, this indicates that a sufficiently accurate approximation is available with a smaller rank. The step is accepted, but the rank for the next step is set to i.e., the rank is reduced by either 1 or 2. Thus, the rank may only decay slowly. Sudden rank-drops are prohibited. For the initial values in the next time step we use To prevent rank-oscillations, rank reduction is prohibited within the first 10 steps after an augmentation step. 3. Persistent case If σ r k ≥ tol > σ r k +1 , the time step is accepted and the same rank r k+1 = r k is used for the next one.

Choice of tolerance
It remains to find a suitable tolerance parameter tol. The global error of our integrator is a combination of a time-discretization error and a low-rank approximation error. We suggest to choose the rank such that the low-rank error is of about the same size as the time discretization error, but does not exceed the latter. If the time discretization error is large, the low-rank error is allowed to be large as well, and hence the approximation rank might be chosen small. If the time discretization error is small, then the approximation rank needs to be sufficiently large for an equally small low-rank error. By this procedure, we hope to ensure that the rank-adaptivity does not impair the convergence order of the integrator. To balance the low-rank error with the time discretization error, we need to approximate the global time discretization error. This comprises an estimate of the local error S ∈ C (r +1)×(r +1) , functions for products with A, tolerance tol} 4: 5: ready = False 6: while not ready do 7: r = r + 1 8: choose u ∈ C m orthonormal to U (e.g., random) 9: choose v ∈ C n orthonormal to V (e.g., random) 10: compute U = U * , S = S * , V = V * as in (12)  11: compute singular values σ 1 , . . . , σ r +1 of S 13: ready = ( σ r +1 < tol) 14: end while 15: S ∈ C (r +1)×(r +1) } 18: end function and a simple model for the global error. While the error analysis of the projector-splitting integrator [16] and of the LRLF scheme [12,Theorem 6] show exponential error growth w.r.t. time, numerical experiments indicate that this is a far too pessimistic bound, and that piecewise linear accretion in time is a more realistic scenario. Since the estimation of the local error induces some computational overhead, we keep an estimate for M steps before we recompute it. To be more precise, let e be an approximation to the local error at time t M+1 . Since we assume that the local error is constant for the next M time steps, the global error at time t M+ j is approximated by E + je , j = 1, . . . , M, where E is defined recursively by This simple technique worked very well in numerous numerical simulations. Note that the linear model is a conservative choice in the sense that if the error growth is faster than linear (e.g., quadratic or even exponential), then we underestimate the global error which enforces the integrator to use a larger rank. We thus still compute a numerical solution where the low-rank approximation does not impair the time integration error. The tolerance threshold tol is then determined heuristically by the following steps: 1. Estimation of the local error (every M steps): Starting from an approximation A M ≈ A(t M ), we compute an approximation A M+1 to A(t M+1 ) by performing one integration step with step size τ and rank r . Additionally, we perform two time steps with step size τ 2 and the same rank r , starting again from the initial value A M . By this, we obtain an alternative approximationȂ M+1 ≈ A(t M+1 ). Assuming that the method converges with order p ∈ N in time, we apply Richardson extrapolation [9, Section II.4] to estimate the local error e at t M+1 as cf. [7,Section 5]. The global error is modeled as 2. Estimation of the low-rank error: If σ 1 ≥ · · · ≥ σ n ≥ 0 are the singular values of the exact solution A(t k+1 ), then the rank-r k+1 best-approximation A best k+1 to A(t k+1 ) fulfills so that 3. Tolerance threshold: The parameter tol k is set by equating the right-hand sides of (13) and (14) for k = M + j, j = 0, . . . , M − 1. This yields the condition 4. Initial rank: An obvious choice for the initial rank for the integration is r 0 = rank(A 0 ). However, if the rank of A 0 is very small, this may not necessarily hold for the rank of the exact solution A(t), even for small t. On the other hand, if the rank of A 0 is large, this choice is also questionable. In our implementation, we first perform ν integration steps (with ν small, e.g., ν = 5) with an initial rank r 1 given by the user (say r 1 = 5). Rank reduction is disabled in this phase. Then let r * denote the number of singular values of A ν greater than or equal to tol ν defined in (15). If r * < r 1 , we continue the integration with r ν+1 = r * . Otherwise, we multiply r 1 by 2 and rerun the initializing process, until r * < r 1 holds.

Rank-adaptive algorithms
The rank-adaptive version of the projector-splitting integrator Algorithm 1 is called RAPSI for rank-adaptive projector-splitting integrator in the following. A single step of the RAPSI scheme is given in Algorithm 4. The rank-adaptive version of the LRLF scheme is derived by replacing the PSI routines by the RAPSI routines. We name this new integrator rank-adaptive LRLF (RALRLF) method. It is presented in Algorithm 5.
The rank-adaptive counterpart of the PSI-stiff scheme is named RAPSI-stiff. Since the linear subproblem preserves the rank, rank-adaptivity is only applied in the integration of the nonlinear subproblem (7b).

Numerical experiments
We now report on numerical experiments for matrix differential equations resulting from space discretizations of PDEs on a rectangular domain For simplicity, we use a uniform mesh with n grid points in x-and m grid points in y-direction, i.e., Errors of the low-rank solutions are measured w.r.t. numerically computed reference solutions. Details are given in the respective subsections. Since we are only interested in the time discretization error, reference solutions and low-rank solutions are computed compute rank-(r A + 1) best-approximation A = U S V H , S = diag( σ 1 , . . . , σ r A , σ r A +1 ) to A 0 6: compute rank-(r B + 1) best-approximation B = T R W H , R = diag(ρ 1 , . . . , ρ r B , ρ r B +1 ) to B 0 7: t 0 = 0 8: for k = 1, . . . , n do 9: t k = t k−1 + τ 10: 11: B-step: A-step: return U , S, V 21: {output: factors U , S, V of rank-r approximation to exact solution A(t n ) of (2) with U ∈ V m,r , 22: V ∈ V n,r , S ∈ C r ×r } 23: end function on the same spatial grid. The computation of the tolerance threshold as explained in Sect. 3.2 is performed with M = 100, i.e., every 100 steps we perform four additional steps with step size τ 2 , so that we increase the computational effort by 4%. Coosing a smaller value of M may sometimes be advantageous to reduce the global error.
All algorithms have been implemented in Python and were performed on a computer with an Intel(R) Core(TM) i7-7820X @ 3.60GHz CPU and 128 GB RAM storage. The codes are available from [25].
In [29], (18) was solved with the linearized second-order backward differential scheme (LBDF2). A fixed-rank dynamical low-rank integrator for (18) was proposed in [32], based on considerations from [24]. We compute low-rank solutions of (18) with PSI-stiff and RAPSI-stiff. The solution to the linear subproblem of (18), which is of form (8), is computed with the Krylov subspace method proposed in [18].
An efficient implementation of products F( A)E of the right-hand side F(A) in (18) with a skinny matrix E is crucial. For the linear part, this is achieved by computing the matrix products in successively from the right to the left. The implementation of the cubic nonlinear part of F is more involved, see "Appendix A.2".
In our first experiment, we use the same parameter values as in [32], namely L x = L y = 10, m = n = 512, ν = η = κ = ξ = γ = 1, T = 1, α = 1.2, β = 1.9, and the initial value The (full-rank) reference solution is computed with LBDF2 on the same spatial grid, with time step size τ = 10 −4 . Figure 1 shows the relative global errors between the reference solution A and the respective low-rank solutions A at t = T for different step sizes τ . Convergence order one is observed for the PSI-stiff scheme. For large step sizes, the approximations computed with the RAPSI-stiff method exhibit large errors. This is explained by the behavior of the singular values, cf. Fig. 1 (right picture). The time-discretization error is overestimated in this experiment, so that the tolerance threshold becomes so large that the second largest singular value is discarded. The induced low-rank error is then of magnitude 10 −1 . If the parameter M is reduced to 10, this unfortunate rank reduction vanishes. However, reducing M increases the workload for the updates of tol, while M = 100 worked well in all other experiments and also in this experiment for smaller step sizes. For our second example, we choose the second parameter set from [32], L x = L y = 8, n = m = 512, ν = κ = 1, η = 0.5, ξ = −5, γ = 3, α = 1.2, β = 1.9, and the initial values . . . , m − 1, j = 1, . . . , n − 1. The relative global errors at T = 1 are displayed in Fig. 2. In contrast to the previous experiment, order one is observed for both integrators. Now the error curves for the fixed-rank integrator and its rank-adaptive variant align almost perfectly, and the singular values follow nicely the trajectories of the singular values of the reference solution.
In our experiment, the reference solution to the problem was again computed with the LBDF2 method, using the step size τ = 2 · 10 −5 . Figure 3 shows the results for the parameter values from [30], L x = L y = 10, n = m = 512, η = 1, ξ = −2, T = 0.2, α = 1.2, β = 1.9, and the initial value Again, the relative global error curves match almost perfectly for both low-rank methods, and are also clearly indicating convergence of order one. The results for other choices of α and β are available in [25].

Laser-plasma interaction
As an example for second-order problems, we consider a reduced model of laserplasma interaction from [14,15,26]. It is given by a wave equation with spacedependent cubic nonlinearity on a bounded, rectangular domain given in (16) with periodic boundary conditions. After space discretization according to [26,Section 4.1.3], we obtain the matrix differential equation with initial values A(0) = A 0 and A (0) = B 0 given by where x j , y i are defined in (17) and i = 1, . . . , m, j = 1, . . . , n. The discrete Laplacian L acts on A(t) via where D x ∈ R n×n denotes the symmetric Toeplitz matrix with first row  Equation (19) describes the propagation of a laser pulse with wavelength λ 0 in the direction of the positive y-axis through vacuum and through a strongly localized plasma barrier. The plasma is located between y = 50π and y = 300π and has constant density 0.3. The localization is modeled by the matrix χ ∈ R m×n with entries The interaction between the pulse and the plasma is modeled by a cubic nonlinearity. As in [15] we use the parameters λ 0 = π , l 0 = 10π , w 0 = 100π , L x = 300π , and L y = 600π .
Our numerical experiments were carried out for t ∈ [0, 600π ] with n = 1024 and m = 8192 discretization points in transversal and longitudinal direction, respectively. The reference solution was computed with the Gautschi-type method from [26], with step size τ 0 = L y /(80m). For the low-rank solutions we used the step sizes τ = 2 k τ 0 , k = 2, . . . , 6. The algorithms LRLF and LRLF-semi were performed with fixed ranks r A = r B = 4.
In [12] it was obeserved that choosing the weights ω i in LRLF-semi according to the direction of motion can improve the approximation significantly. Since the laser pulse moves mainly in longitudinal direction, we therefore used the weights , ω 2 2 = 0, ω 2 3 = 1 3 The left picture in Fig. 4 shows the relative global error at T = 600π between the reference solution and the different low-rank integrators. Second-order convergence is observed in all cases, and the integrators designed for semilinear problems yield better approximations than the other methods. The accuracy of the rank-adaptive schemes is comparable to those of the fixed-rank integrators, showing nicely that the heuristics works well for this example, cf. Fig. 5.
In physics, the maximal intensity max i, j |A i j (t)| 2 of the propagating pulse over time is sometimes of higher interest than A itself. In the right picture of Fig. 4, the absolute error between the maximal intensity of the numerical pulse computed with the Gautschi-type integrator and the maximal intensity of the approximations obtained by the low-rank integrators is displayed.

Sine-Gordon equation
In our last experiment, we consider the two-dimensional sine-Gordon equation on the domain from (16) with homogeneous Neumann boundary conditions [1]. Using finite differences of second order on the grid (17) with L x = L y = 7, n = m = 1001 in both x-and y-direction, we obtain the semi-discretized matrix differential equation Here, sin(A) denotes the entrywise evaluation of the sine function. The matrix D is given by The storage-economical evaluation of the products sin( A)E and sin( A) H E is presented in "Appendix A". There is no preferred direction of propagation, so that we used the weights in (10)  The reference solution is computed by the leapfrog scheme on the same spatial grid with time step size τ = 2.5 · 10 −5 . Figure 6 shows the relative global errors between the low-rank approximations and the reference solution. Convergence order two is observed for all methods. The fixed-rank integrators are slightly more accurate than their rank-adaptive pendants, probably because they use a higher rank. In a second setting, we consider the symmetric perturbation of a static line soliton [2, Section 3.1.2] with i j = 1, (B 0 ) i j = 0, and , i, j = 0, . . . , m. Figure 7 shows a similar behavior of the methods as in the first setting. For the RALRLF scheme however, the initial rank is rather large, and drops significantly after a few steps. This is caused in the routine for determining an appropriate initial rank. As explained in Sect. 3.2, the initial guess r 1 = 5 is doubled repeatedly until the criterion for continuing the integration beyond the first ν steps is satisfied. In this experiment, an initial rank of ∼ 23 is adequate. Therefore, the guesses 5, 10, and 20 are rejected, until r 1 = 40 is accepted and rank reduction applies in the subsequent integration steps.

Conclusion and outlook
In the present paper, we developed adaptive schemes for dynamical low-rank integrators for first and second-order matrix differential equations. The performance of these schemes have been illustrated by numerical experiments. Both the projector-splitting integrator and the unconventional robust integrator have been successfully adapted to first-order tensor differential equations, cf. [6,20,21] and references therein. We are confident that the strategy for the adaptive algorithm for matrix differential equations also works for the tensor case. This is part of ongoing research and will be reported in the future.
where U jk = U j • U k • U , The computational cost is further reduced by exploiting the symmetry in j and , The product ( A • A • A) H E with E ∈ C m×r can be computed analogously.
Additional simplifications apply to the product ( χ • A • A • A)E with χ ∈ C m×n given in (20). It satisfies where 0 n , 1 n are the vectors of length n filled with zeros and ones, respectively, and 0 n× p and 1 n× p the matrices of dimension n × p with all entries being zeros and ones, respectively. From (22) and (25), we obtain where ϑ( U j ) denotes the restriction of U j to its ηth to ξ th entries. Here, we made use of the following property of the Hadarmard product, cf. [27, Section 2]: If A ∈ C m×n , x ∈ R m , and y ∈ R n , then Hence, it suffices to compute and sum up the small matrices ϑ( U j )( V H j E). Likewise, using (24), we have The implementation of ( χ • A • A • A) H E, E ∈ C m×r , is realized in the same manner.