Solution formulas for differential Sylvester and Lyapunov equations

The differential Sylvester equation and its symmetric version, the differential Lyapunov equation, appear in different fields of applied mathematics like control theory, system theory, and model order reduction. The few available straight-forward numerical approaches when applied to large-scale systems come with prohibitively large storage requirements. This shortage motivates us to summarize and explore existing solution formulas for these equations. We develop a unifying approach based on the spectral theorem for normal operators like the Sylvester operator S(X)=AX+XB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {S}}(X)=AX+XB$$\end{document} and derive a formula for its norm using an induced operator norm based on the spectrum of A and B. In view of numerical approximations, we propose an algorithm that identifies a suitable Krylov subspace using Taylor series and use a projection to approximate the solution. Numerical results for large-scale differential Lyapunov equations are presented in the last sections.


Introduction
For coefficient matrices A ∈ C n×n and B ∈ C m×m , an inhomogeneity C ∈ C n×m , and an initial value D ∈ C n×m , we consider the differential matrix equatioṅ and provide formulas for the solution X with X (t) ∈ C n×m . Equation (1) is commonly known as differential Sylvester equation (as opposed to the algebraic Sylvester equation AX + B X +C = 0). In the symmetric case when B = A T , Eq. (1) and its algebraic counterpart are called differential and algebraic Lyapunov equations, respectively. In what follows, we will occasionally abbreviate Sylvester or Lyapunov equations by SLE. In particular the differential Lyapunov equation is a useful tool for stability analysis and controller design for linear time-varying systems [2]. Equilibrium points of the differential Lyapunov equation, namely solutions of the algebraic Lyapunov equation, are used to construct quadratic Lyapunov functions for asymptotically stable linear time-invariant systems [33,Thm. 7.4.7]. The controllability and observability problem for linear time-varying systems is strongly connected to the solution of the differential Lyapunov equation [9,, [19,. Other important applications lie in model order reduction [3] or in optimal control of linear time-invariant systems on finite time horizons [30]. Despite its importance, there have been but a few efforts to solve the differential Sylvester / Lyapunov or Riccati equation numerically, see [6,7,16,21,25,26,31,36]. These algorithms are usually based on applying a time discretization and solving the resulting algebraic equations. Thus, even if the algebraic SLE are solved efficiently, the storage needed for the discrete solution at the time instances makes these direct approaches infeasible for large scale settings. Recently, Krylov subspace based methods were proposed in [12][13][14]24].
In an attempt to overcome this shortage, we revisit known solution formulas, develop alternative solution representations, and discuss their suitability for numerical approximations. We start with deriving a spectral decomposition for the Sylvester operator S which allows functional calculus. We obtain formulas for the operator norm S as well as for S −1 and e S . This recovers previously known solution formulas. It will turn out that, in terms of efficiency, this solution representation is not well suited for approximation in general but, in special cases, allows for the construction or computation of exact solutions.
As a step towards efficient solution approximation, we use Taylor series expansions to identify suitable Krylov subspaces. For the differential Lyapunov equation with stable coefficient matrices and symmetric low-rank factored inhomogeneity, it is wellknown that the solution of the algebraic Lyapunov equation spans a Krylov subspace under these assumptions. We split the solution of the differential Lyapunov equation in a constant and time dependent part, where the constant part is the solution of an algebraic Lyapunov equation. We approximate the time dependent part using the subspace spanned by the solution of the algebraic Lyapunov equation. The resulting algorithm overcomes the essential problems with storage consumption. Numerical results are presented in Sect. 7 and "Appendices B and C".

Preliminaries
In this section, we introduce the considered equations and the Sylvester operator, set the notation, and recall basic results.
The spectrum of a matrix A ∈ C n×n is denoted by Λ(A) ⊆ C. A matrix is called stable if its spectrum is contained in the left open half plane C − . The Frobenius inner product on C n×m is given by A,  an open interval with t 0 ∈ I , A ∈ C(I , C n×n ), B ∈ C(I , C m×m ), C ∈ C(I , C n×m ) and D ∈ C n×m . The differential Sylvester equatioṅ Here Φ A (t, t 0 ) and Φ B H (t, t 0 ) are the unique state-transition matrices with respect to t 0 ∈ I defined byΦ An application of Theorem 1 specifically to the autonomous case with constant coefficients by simply replacing the state transition matrices with the matrix exponentials yields the next result.
Theorem 2 Let I ⊆ R be an open interval with t 0 ∈ I , A ∈ C n×n , B ∈ C m×m , C ∈ C(I , C n×m ) and D ∈ C m×n . The differential Sylvester equatioṅ has the unique solution Next, we restate basic properties of the Sylvester operator S : C n×m → C n×m , which, for given A ∈ C n×n and B ∈ C m×m , is defined by -e tS = e tH e tV for all t ∈ R, for any A ∈ C n×n and B ∈ C m×m .
Proof The first claim can be confirmed by direct computations. The second claim is a standard result for commuting linear operators.
By Lemma 1, formula (2) can be rewritten as

Spectral decomposition of the Sylvester operator
In this section we show that the Sylvester operator S, as defined in (3), is a normal operator if A and B are diagonalizable and if a suitably chosen inner product on a Hilbert space is considered. The inner product depends on the decomposition of A and B. Nevertheless, this approach will enable us to apply the spectral theorem and to derive a solution formula for the differential and algebraic SLE. This resembles the formulas of [11,Ch. 4.1.1], [18,34]. Those results were obtained by inserting the spectral decomposition into the SLE and by applying suitable algebraic manipulations and using the unrolled Kronecker representation of the SLE. Our strategy is to decompose the operator S first and then using functional calculus to obtain formulas for e S and S −1 . The eigenspaces of S can be constructed from the eigenspaces of A and B. The choice of the inner product ensures that the eigenvectors are orthonormal and S becomes a normal operator.
is an orthonormal basis of C n×m with respect to ·, · U ,V .
(iii) The adjoint operator S * : Proof It is straightforward to see that ·, · U ,V is an inner product on C n×m . Because of the identity the matrices u i v H j ∈ C n×m are orthogonal with respect to ·, · U ,V and therefore linearly independent. As dim(C n×m ) = n ·m, the tuple ( is an orthonormal basis of C n×m . Finally we show that, in this inner product, the adjoint of S is defined as This means S and S * commute and, therefore, by definition, S is normal. Now that we have an inner product on a Hilbert space for which S is normal, the second step is to compute the spectral decomposition of S. The spectral decomposition allows functional calculus and we get a formula for S −1 and e tS . Since for normal operators, the operator norm is its spectral radius, we directly get a formula for the induced operator norm of S. We remark that in the case that A and B are unitarily diagonalizable Lemma 2 also holds when ·, · U ,V is replaced by the standard Frobenius inner product ·, · F . Lemma 3 (Spectral Decomposition of S) Let the assumptions of Lemma 2 hold. Moreover let α 1 , . . . , α n ∈ C and β 1 , . . . , β m ∈ C be the diagonal entries of D A and D B H , respectively. Then we have Representing S(X ) ∈ C n×m as well as X ∈ C n×m with respect to the orthonormal basis u i v H j and exploiting linearity of S and ·, · U ,V yield (ii) The claim about the norm follows from a direct application of the fundamental functional analytical result on compact normal operators, see, e.g., [ , and obtain the formula for S −1 under the additional assumption that α i + β j = 0.
Using the spectral decomposition and functional calculus we find that, under the assumptions of Lemma 2, the solution of the differential Sylvester equatioṅ has the form with the involved scalar integrals given explicitly as:

Variation of constants
The application of the variation of constants formula leads to yet another solution formula for the SLE (1).

Lemma 4 (Variations of Constants
The differential Sylvester equatioṅ has the solution Proof Because of Λ(A) ∩ Λ(−B) = ∅, the inverse S −1 exists and we can rewrite the solution formula (4) as and confirm that From formula (6), we find that the solution can be written as the solution of the algebraic Sylvester equation and a time dependent part. We will make use of this fact in the numerical scheme that we propose in Sect. 6.

Solution as Taylor series
In this section we use Taylor series to derive a solution formula. From this we can read off suitable Krylov subspaces for our projection approach in the next section.
Lemma 5 (Taylor series solution formula) Let A ∈ C n×n , B ∈ C m×m , C ∈ C m×n , D ∈ C m×n . The differential Sylvester equatioṅ has the unique solution Proof First observe that The series converges absolutely and since (C n×m , || · ||) is a Banach space, the series converges for every t ∈ R. Therefore its radius of convergence is infinity, X is continuously differentiable and can be differentiated term-wise. Since, furthermore, If we assume that D and C are given in factored form, then we can exploit this to rewrite the truncated series in a closed form of a matrix product.
Then, having truncated the two parts of the series (7) after m 1 and m 2 summands, respectively, we can rewrite the solution approximation as with the explicit representation of the sums

Feasible numerical solution approaches
In this section, we briefly note that, for various reasons, all presented solution representations are not feasible for a straight-forward numerical approximation, in particular in a large-scale sparse setting. A common reason is that none of the formulas supports a sparse representation of the solutions such that exorbitant amounts of memory will be required.
Limitations in memory will doubly affect the solution representation through the spectral decomposition (5) since also the basis matrices U and V are generically dense matrices. Apart from that, the computation of a spectral decomposition is typically computationally expensive and can be ill conditioned. Nonetheless, the spectral decomposition formula is useful to construct exact solutions for given coefficients with known spectral decompositions.
Another issue is the unfeasible computation of the full matrix exponential in all variants (2), (4), and (6) of the variation of constants formula. A possible remedy is the approximation of the action of the matrix exponential on a low-rank matrix in a Krylov subspace.
The approach to the solution via a Taylor series (see Sect. 5) seems best suited for the large-scale case since, at least in the symmetric case, the formulas provided in Remark 1 allow for a solution representation in factored form with the original coefficients. One problem here is that the truncated Taylor series only leads to good approximations locally around the point of expansion.
We will, however, exploit and combine certain parts of the solution representations to propose an algorithm for fast and memory efficient solution approximations.
We consider the stable, linear time-invariant case, i.e. we assume that A ∈ R n×n , B ∈ R n× p , and Λ(A) ⊆ C − . We consider the differential Lyapunov equatioṅ By Lemma 4, we have that the solution splits into a constant part and a time dependent part. From Remark 1, we infer that the solution is contained in a Krylov subspace. We combine both observations in the following numerical solution approach: 1. Factors of the time constant part as Krylov space basis With A stable, the associated algebraic Lyapunov A T X + X A + B B T = 0 has a unique symmetric positive semi-definite solution X ∞ that can be written in factored form X ∞ = Z ∞ Z T ∞ , with Z ∞ ∈ R n×q and rank(Z ∞ ) = q ≤ n. Moreover, since A is stable, it holds (see [9,Ch. 13

Factors of the time dependent part evolve in the same Krylov space
And since by (8) also X (t) evolves in the same Krylov space 1 , as does the differencẽ X (t) and, thus,Z (t).
3. Orthogonalize the basis and compute the time dependent factors By means of a (reduced) Singular Value Decomposition of Z ∞ , we obtain orthogonal matrices Q ∞ ∈ R n×q and V ∞ ∈ R q×q with range(Q ∞ ) = range(Z ∞ ) and from which we deduce that We define z(t) = e t Q T ∞ A T Q ∞ S ∞ ∈ R q×q and find that z can be obtained by solvinġ which is a matrix valued ODE that can be solved column-wise or by computing the matrix exponential e t Q T ∞ A T Q ∞ . The solution of the differential Lyapunov equation is, thus, given by

Remark 2
The differential equation for z is of size q × q, which can be much smaller than n, if the solution of the algebraic Lyapunov equation X ∞ = Z ∞ Z T ∞ has low-rank. Moreover, the orthogonalization of the basis allows for the detection of the numerical rank and for a compression of Z ∞ through truncating singular values that are smaller than a certain threshold.
We further note that, with minor adjustments, all arguments also hold for the generalized differential Lyapunov equation with M ∈ R n×n nonsingular and AM −1 stable that can accommodate, e.g., a mass matrix from a finite element discretization. We have included the derivation of the relevant formulas, and an example that illustrates a fundamental difference when M is symmetric positive definite in "Appendix A". In summary, the proposed approach reads as written down in Algorithm 1.

Setup
To quantify and illustrate the performance of Algorithm 1, we consider differential Lyapunov equations that are used to define optimal controls for a finite element discretization of a heat equation; see [8] for the model description. Namely, we solve the differential Lyapunov equations: that are defined through matrices M, A ∈ R n×n that are symmetric, M is positive definite and A stable, B ∈ R n×7 and C ∈ R 6×n . For computing the error, we precomputed the spectral decomposition of (A, M) and constructed the exact solution Algorithm 1 Projection approach for generalized differential Lyapunov equations
7: % Truncate all singular values smaller than tolerance and get truncated low-rank factor: 12: % Compute projected system and solve: 13: if M is symmetric positive definite then 14: We solve the resulting projected ODE system column-wise using MATLAB ODE solvers ode45, ode23, ode113, ode15s, ode23s, ode23t and ode23tb. The solvers ode45, ode23, ode113 are nonstiff solvers, whereas the rest are stiff solvers. We have used the following parameters for the odeset function:

Projection approach
The initial step of Algorithm 1 requires the solutions to the associated algebraic Lyapunov equations. For this task we call MEX-M.E.S.S. that iteratively computes the solutions up to the following absolute and relative residuals and The achieved values for the different test setups as well as the number of columns of the corresponding Z ∞ after truncation (see Step 7 of Algorithm 1), that define the dimension of the reduced model for the time dependent part, are listed in Tables 1  and 2. We report the absolute and relative errors where X is the numerical solution obtained from Algorithm 1 with various ODE solvers and where the reference solution X re f was obtained from the spectral decomposition As the matrices A and M are real symmetric and M is positive definite, the eigenproblem for (A, M) is a generalized symmetric definite one. Therefore all eigenvalues are real and the corresponding system of eigenvectors can be chosen real and orthogonal with respect to M, that is V T M V = E n,n . The representation of the reference solution simplifies as follows We plot the numerical errors and ||X re f (t)|| 2 on the initial short time interval [0, 10], where most of the evolution is happening, and on the full time interval [0, 4500] in "Appendix B", Figs. 4a-d, 5a-d, 6a-d, 7a-d, 8a-d and 9a-d.
In view of the performance of the different ODE solvers, we can interpret the presented numbers and plots as follows: the solver ode15s, which is a stiff solver of variable order, performs best in time and accuracy. Due to the stiffness, the error is oscillating for the solvers ode45, ode23, ode113. Note that the discrete Laplacian that is encoded in the coefficient matrix A becomes stiffer with a finer space discretization, i.e. for larger n. Accordingly, the computational times for the non-stiff solvers grow with n at a higher rate than the stiff solvers; see Fig. 2.
The solution of (DLE-2) itself is large in norm which makes the relative error stagnate around the prescribed tolerance and the absolute error comparatively large; see the plots in "Appendix B.2". Particularly, the non-stiff solvers (and ode15s) achieve this error level and the oscillations as stiffness is dominated by the approximation error. This might be the reason, why for (DLE-2), despite the fact that the coefficient matrices are still stiff, the non-stiff solvers perform better than the stiff solvers (except ode15s). Nonetheless, again the computation times for the non-stiff solver grow at a higher rate with the increasing stiffness that comes with increasing n; see Fig. 2.
Because X (t) → 0 for t 0, the plots for the relative error spread out for small times.
Except ode15s, there is no general rule which solver performs better in terms of computational time; see Fig. 2. The timings may change, when different relative and absolute error tolerances are used in the MATLAB odeset function.
Finally, we want to make the following remark: instead of integrating the projected ODE (10) which is linear with constant coefficients of moderate dimension, one may consider using the Schur-Parlett [17,Ch. 10] algorithm to compute the matrix expo- Fig. 1 Link to code and data nential. The initialization efforts of the Schur-Parlett algorithm will pay off, if the matrix exponential has to be evaluated for many different values of t. Also, asymptotically, the storage requirements for z(t) will be lifted, since the matrix exponential for a given t can be computed on demand. Nonetheless, we used the ODE approach, which we think is more efficient because of the sophisticated MATLAB ODE solver implementations that come, e.g, with step size selection methods integrated.
The code of the implementation with the precomputed spectral decompositions needed to construct the exact solution is available as mentioned in Fig. 1.

Comparison with backward differentiation formulas
To benchmark our method, we have run comparison tests with the MATLAB implementation of the Backward Differentiation Formulas / Alternating Direction Implicit (BDF/ADI) scheme as developed in [29] (see also the "Appendix C" for a short summary of the algorithm). The numerical experiments were conducted on the same machine, with the same MATLAB version and the same model as described in Sect. 7.1.
In contrast to the Algorithm 1 that needs to solve only a single algebraic Lyapunov equation for the initialization, the BDF/ADI approach solves a Lyapunov equation in every time step. Moreover, the numerical solution is stored in L DL T -format, i.e., in terms of the factors of X (t k ) ≈ L k D k L T k and L k ∈ R n×l k , which grow at least linearly with the size of n and the number of time steps. For this reason, we had to restrict our numerical experiments with BDF/ADI to the interval [0, 100] and consider only the model of the smallest size n = 1357. As for the test of our Algorithm 1 in "Appendix B", for computing the error, we used an exact solution X re f based on the spectral decomposition of (A, M).
We have compared various BDF methods that we abbreviated as BDF1, BDF2, BDF3, BDF4, BDF5 and BDF6, where the number denotes the order s of the method. We used the constant time step sizes τ k := h ∈ {2 −4 , 2 −6 , 2 −8 } for our computations.
In "Appendices C.1 and C.2", we plot the relative and absolute errors compared to the solution obtained by the spectral decomposition, c.f. Figs. 10a-f and 11a-f. For comparison, we also plot the error of the numerical solution obtained by Algorithm 1 and the MATLAB ODE solver ode15s. We list the computational times for performing the BDF/ADI method based computations in Sect. 7.3.
As the actual solution X converges towards the solution of the algebraic Lyapunov equation, also the numerical errors show a decay towards a certain error level. Decreasing the step size lowers this level accordingly. This effect is visible only for methods of lower order (s ≤ 2). For higher orders, the error level stagnates and the error rather shows oscillations, which are likely due to errors in the solution of the Lyapunov equations. Since all BDF methods are stiff solvers, there are no oscillations of higher frequency to observe. The error levels are comparable to the level reached by our approach with ode15s.
The computational timings for the BDF/ADI solvers are nearly the same for same step sizes; cf. Fig. 3. Accordingly, the higher order methods clearly outperform the low-order methods.
As for the comparison to our approach, we note that the reported timings were for the longer time-interval [0, 4500]. However, since the transient behavior of the solution is confined to a short initial phase and since the ODE solvers have a step size control, a restriction to a shorter interval would not change the timings by much.
For the test case with (DLE-1), the BDF/ADI schemes achieved the same accuracy as our approach already with the coarsest step size; see Fig. 10a-f. Thus, we compare the timings in the left most columns in Fig. 2 (n = 1357) and Fig. 3 (h = 2 −4 ) to conclude that our approach with the best performing ODE solver is about 25 times faster.
As for the test case with (DLE-2), the BDF/ADI schemes reached the same accuracy only for the finest chosen step size; see Fig. 11a-f. Thus, comparing the right most column in Fig. 3 (h = 2 −8 ) to the left column in the right plot of Fig. 2 (n = 1357), we find that our Algorithm 1, again, is faster by a factor of 30.
To conclude the comparison, we add the following remarks. With the costs of solving Lyapunov equations, apart from the increased memory requirements, the BDF/ADI scheme will also become significantly more costly for larger system sizes. As opposed to our Algorithm 1 however, the BDF/ADI scheme also applies for the time varying case as well as for nonzero initial conditions. Moreover, as ode15s in fact uses the Fig. 3 Timings of the different BDF/ADI solvers BDF formulas but with variable order s and variable time step, one may consider integrating similar error control mechanisms in the BDF/ADI approach to improve performance.

Conclusions
We presented several solution formulas for the differential Sylvester and Lyapunov equations. For the autonomous stable differential Lyapunov equation, we proposed a numerical algorithm that combined certain aspects derived from the solution representations. The main feature of the algorithm is the projection of the time dependent part onto a suitable subspace of lower dimension. Only this makes the numerical solution feasible in terms of memory requirements. As for the computational time, the projected system can be directly solved with optimized ODE solvers like the ode-suite in MATLAB. Moreover, its structure allows for column-wise computation of the factors and, thus, for straight-forward parallelization.
We illustrated the performance of the algorithm in an example that was derived from a finite-element discretization of a heat equation. The achieved accuracy is fully satisfactory. The greatest benefit, also in contrast to existing numerical schemes, is the low memory requirement.
The possible extension to the unstable or unsymmetric as well as to the nonautonomous case is not straightforward since the space in which the solution evolves has not been found to span a suitable invariant Krylov subspace and is possibly of high dimension. Moreover, the solution of AX = C (which is the special case of AX + X B = C with B = 0) is unlikely to span an A-invariant subspace at all, meaning that the span of the solution of the algebraic Sylvester equation is not suited for the differential equation. Thus, for the general case, a different ansatz is needed. The same is true for time-dependent coefficients or non-zero initial conditions. Here, a remedy might be structured approaches for the case that time dependence comes, e.g., from a low-rank update.
In the unstable case, apart from the fact that the algebraic Lyapunov equation may not have a unique solution, there can be modes that grow exponentially in time. For this case it may be worth investigating whether a projector that identifies the stable part of the underlying mathematical model can be efficiently incorporated in the proposed solution approach.  (13) is given by

Example 1 Let
A reduced SVD of Z ∞ is given by System (17) reads now as which is no ordinary differential equation and has multiple solutions. In contrast system (18) is given byż The unique solution of the generalized differential Lyapunov equation (12) is given by
The s-step BDF method applied to the DLE 9 is given by where α j and β are coefficients of the BDF method and can be found in [15]. The parameter s is the order of the BDF method. We recall that for s > 6, the method is not numerical stable, and for s = 1, the BDF method coincides with the implicit Euler method. A minor rearrangement shows that the current iterate X k can be obtained as the solution of the algebraic Lyapunov equation Since for s ≥ 2, certain coefficients α j , j ≥ 1 are positive, the algebraic Lyapunov equation (19) has a symmetric but possibly indefinite right-hand side, which makes the standard ADI method infeasible. For this reason a L DL T -decomposition for the numerical solution is proposed and suitable modifications of the ADI method have been developed; [5,26]: Assume that X i ≈ L i D i L T i for i = 0, . . . , k − 1, L i ∈ R n×l i , D i ∈ R l i ×l i and l i n, then the right hand side can be factored as Now the L DL T -type ADI method can be used to determine X k ≈ L k D k L T k . The BDF/ADI methods can also be extended to generalized differential Lyapunov equations in a similar way [28].