Numerical solutions to large-scale differential Lyapunov matrix equations

TL;DR

This paper introduces two novel numerical methods for solving large-scale differential Lyapunov matrix equations with low rank terms, combining integral approximations and Krylov subspace projections to improve efficiency and accuracy.

Contribution

The paper proposes two new approaches for solving large-scale differential Lyapunov equations, integrating integral approximations and Krylov subspace techniques.

Findings

New theoretical results on differential Lyapunov equations

Numerical experiments demonstrating method effectiveness

Efficient algorithms for large-scale problems

Abstract

In the present paper, we consider large-scale differential Lyapunov matrix equations having a low rank constant term. We present two new approaches for the numerical resolution of such differential matrix equations. The first approach is based on the integral expression of the exact solution and an approximation method for the computation of the exponential of a matrix times a block of vectors. In the second approach, we first project the initial problem onto a block (or extended block) Krylov subspace and get a low-dimensional differential Lyapunov matrix equation. The latter differential matrix problem is then solved by the Backward Differentiation Formula method (BDF) and the obtained solution is used to build the low rank approximate solution of the original problem. The process being repeated until some prescribed accuracy is achieved. We give some new theoretical results and present some numerical experiments.

Tables

Table 1. Table 1: Coefficients of the p 𝑝 p -step BDF method with p ≤ 3 𝑝 3 p\leq 3 .

$p$	$\beta$	${\alpha }_0$	${\alpha }_1$	${\alpha }_2$
1	1	1
2	2/3	4/3	-1/3
3	6/11	18/11	-9/11	2/11

Table 2. Table 2: runtimes and residual norms for EBA-exp and EBA+BDF(2)

size($A$)	EBA-exp	EBA-BDF(2)	Residual norm
$2500\times 2500$	$3.23$ s	$31.9$ s	$𝒪\left({10}^{-9}\right)$ $\left(m=16\right)$
$6400\times 6400$	$5.2$ s	$81.6$ s	$𝒪\left({10}^{-9}\right)$ $\left(m=19\right)$
$10000\times 10000$	$5.6$ s	$168$ s	$𝒪\left({10}^{-8}\right)$ $\left(m=19\right)$
$22500\times 22500$	$11.8$ s	$1546$ s	$𝒪\left({10}^{-8}\right)$ $\left(m=23\right)$

Table 3. Table 3: runtimes and residual norms for EBA-exp and EBA-BDF(2)

size($A$)	EBA-exp	EBA-BDF(2)	Residual norms
$2500\times 2500$	$1.0$ s	$8.0$ s	$𝒪\left({10}^{-11}\right)$ $\left(m=11\right)$
$6400\times 6400$	$4.9$ s	$14.4$ s	$𝒪\left({10}^{-14}\right)$ $\left(m=11\right)$
$10000\times 10000$	$11.5$ s	$29.7$ s	$𝒪\left({10}^{-13}\right)$ $\left(m=11\right)$
$20000\times 20000$	$11.8$ s	$173.4$ s	$𝒪\left({10}^{-13}\right)$ $\left(m=11\right)$

Table 4. Table 4: Optimal Cooling of Steel Profiles: runtimes and residual norms for EBA-exp and EBA-BDF(2)

size($A$)	EBA-exp	EBA-BDF(2)	Residual norms
$1357\times 1357$	$11.1$ s	$515.5$ s	$𝒪\left({10}^{-8}\right)$ $\left(m=6\right)$
$5177\times 5177$	$148.7$ s	$1721$ s	$𝒪\left({10}^{-7}\right)$ $\left(m=39\right)$

Equations

$$\left\{\begin{array}[]{l}\dot{X}(t)=A(t)\,X(t)+X(t)\,A^{T}(t)+B(t)B(t)^{T};\;(DLE)\\ \;X(t_{0})=X_{0},\;\;t\in[t_{0},\,T_{f}],\end{array}\right.$$

$$\dot{X}(t)=A(t)X+XA(t{)}^T+M(t);X(t_0)=X_0$$

$$X(t)={\Phi }_A(t,t_0)X_0{\Phi }_A^T(t,t_0)+{\int }_{t_0}^t{\Phi }_A(t,\tau )M(\tau ){\Phi }_A^T(t,\tau )d\tau .$$

$${\dot{\Phi }}_A(t,t_0)=A(t){\Phi }_A(t,t_0),{\Phi }_A(t_0,t_0)=I.$$

$$X(t)=e^{(t-t_0)A}{X}_0{e}^{(t-t_0)A^T}+{\int }_{t_0}^t{e}^{(t-\tau )A}M(\tau )e^{(t-\tau )A^T}d\tau .$$

$$\left\{\begin{array}[]{c c l}\dot{x}(t)&=&\mathcal{A}(t)x(t)+b(t)\\ x_{0}&=&vec(X_{0})\end{array}\right.$$

$${\mathbb{K}}_m(A,B)=\mathrm{range}(B,AB,\dots ,A^{m-1}B)$$

$${\mathcal{K}}_m(A,B)=\mathrm{range}(A^{-m},\dots ,A^{-1}B,B,AB,\dots ,A^{m-1}B)$$

$${\mathcal{K}}_m(A,B)={\mathbb{K}}_m(A,B)+{\mathbb{K}}_m(A^{-1},A^{-1}B).$$

$$A{\mathcal{V}}_m$$

$$=$$

$$X(t)=e^{(t-t_0)A}{X}_0{e}^{(t-t_0)A^T}+{\int }_{t_0}^t{e}^{(t-\tau )A}B{B}^T{e}^{(t-\tau )A^T}d\tau ,$$

$$Z_m(t)={\mathcal{V}}_m{e}^{(t-\tau ){\mathcal{T}}_m}{\mathcal{V}}_m^{T}B$$

$$e^{(t-\tau )A}B{B}^T{e}^{(t-\tau )A^T}\approx Z_m(t)Z_m(t{)}^T.$$

$$X_m(t)={\mathcal{V}}_m{G}_m(t){{\mathcal{V}}_m}^T,$$

$$G_m(t)={\int }_{t_0}^t{\tilde{G}}_m(\tau ){\tilde{G}}_m^T(\tau )d\tau ,$$

$${\dot{G}}_m(t)={\mathcal{T}}_m{G}_m(t)+G_m(t){{\mathcal{T}}_m}^T+B_m{B}_m^T,t\in [t_0,T_f]$$

$${\mathcal{V}}_m^T{R}_m(t){\mathcal{V}}_m$$

$$F(z)=a_0+\sum _{i=1}^p\frac{a_i}{z-{\theta }_i},$$

$${\tilde{G}}_m(\tau )\approx a_0{B}_m+\sum _{i=1}^p{a}_i[(t-\tau ){\mathcal{T}}_m-{\theta }_{i}I{]}^{-1}{B}_m.$$

$$\parallel R_m(t)\parallel =\parallel T_{m+1,m}{\bar{G}}_m(t)\parallel ,$$

$${\dot{X}}_m(t)=(A-F_m)X_m+X_m(A-F_m{)}^T+B{B}^T.$$

$$X_m(t)\approx {\tilde{Z}}_m(t){\tilde{Z}}_m(t{)}^T.$$

$${\dot{E}}_m(t)=A{E}_m(t)+E_m(t)A^T-R_m(t),$$

$$E_m(t)=e^{(t-t_0)A}{E}_{m,0}{e}^{(t-t_0)A^T}+{\int }_{t_0}^t{e}^{(t-\tau )A}{R}_m(\tau )e^{(t-\tau )A^T}d\tau ,t\in [t_0,T_f].$$

$$\parallel E_m(t)\parallel \le \parallel T_{m+1,m}\parallel \parallel {\bar{G}}_m{\parallel }_{\infty }\frac{e^{2(t-t_0){\mu }_2(A)}-1}{2{\mu }_2(A)},$$

$$\parallel E_m(t)\parallel \le {\int }_{t_0}^t\parallel e^{(t-\tau )A}{\parallel }^2\parallel R_m(\tau )\parallel d\tau .$$

$$\parallel E_m(t)\parallel$$

$$\Vert X(t)-X_m(t)\Vert$$

Full text

1 Introduction

In the present paper, we consider the differential Lyapunov matrix equation (DLE in short) of the form

[TABLE]

where the matrix $A(t)\in\mathbb{R}^{n\times n}$ is assumed to be nonsingular and $B(t)\in\mathbb{R}^{n\times s}$ is a full rank matrix, with $s\ll n$. The initial condition $X_{0}$ is assumed to be a symmetric and positive low-rank given matrix.

Differential Lyapunov equations play a fundamental role in many areas such as control, filter design theory, model reduction problems, differential equations and robust control problems [1, 5]. For those applications, the matrix $A$ is generally sparse and very large. For such problems, only a few attempts have been made to solve (1).

Let us first recall the following theoretical result which gives an expression of the exact solution of (1).

Theorem 1**.**

[1]** The unique solution of the general Lyapunov differential equation

[TABLE]

is defined by

[TABLE]

where the transition matrix $\Phi_{A}(t,t_{0})$ is the unique solution to the problem

[TABLE]

Futhermore, if $A$ is assumed to be a constant matrix, then we have

[TABLE]

We notice that the problem (1) is equivalent to the linear ordinary differential equation

[TABLE]

where $\mathcal{A}=I\otimes A(t)+A(t)\otimes I$, $x(t)=vec(X(t))$ and $b(t)=vec(B(t)B(t)^{T})$, where $vec(Z)$ is the long vector obtained by stacking the columns of the matrix $Z$. For moderate size problems, it is then possible to use an integration method to solve (5). However, this approach is not adapted to large problems. In the present paper, we will consider projection methods onto extended block Krylov (or block Krylov if $A$ is not invertible) subspaces associated to the pair $(A,B)$. These subspaces are defined as follows

[TABLE]

for block Krylov subspaces, or

[TABLE]

for extended block Krylov subspaces. Notice that the extended Krylov subspace ${\cal K}_{k}(A,B)$ is a sum of two block Krylov subspaces

[TABLE]

To compute an orthonormal basis $\{V_{1},\ldots,V_{m}\}$, where $V_{i}$ is of dimension $n\times s$ for the block Krylov and $n\times 2s$ in the extended block Krylov case, two algorithms have been defined: the first one is the well known block Arnoldi algorithm and the second one is the extended block Arnoldi algorithm [7, 25]. These algorithms also generate block Hessenberg matrices ${\bar{\cal T}_{m}}={\cal V}_{m+1}^{T}\,A\,{\cal V}_{m}$ satisfying the following algebraic relations

[TABLE]

where ${\cal T}_{m}={\bar{\cal T}_{m}}(1:d,:)={\cal V}_{m}^{T}A{\cal V}_{m}$ and where $T_{i,j}$ is the $(i,j)$ block of ${\bar{\cal T}_{m}}$ of size $d\times d$, and $E_{m}=[O_{d\times(m-1)d},I_{d}]^{T}$ is the matrix of the last $d$ columns of the $md\times md$ identity matrix $I_{md}$ with $d=s$ for the block Arnoldi and $d=2s$ for the extended block Arnoldi.

When the matrix $A$ is nonsingular and when the computation of $W=A^{-1}V$ is not difficult (which is the case for sparse and structured matrices), the use of the extended block Arnoldi is to be preferred.

The paper is organized as follows: In Section 2, we present a first approach based on the approximation of the exponential of a matrix times a block using a Krylov projection method. We give some theoretical results such as an upper bound for the norm of the error and an expression of the exact residual. A second approach,presented in Section 3, for which the initial differential Lyapunov matrix equation is projected onto a block (or extended block) Krylov subspace. Then, the obtained low dimensional differential Lyapunov equation is solved by using the well known Backward Differentiation Formula (BDF). In Section 4, an application to balanced truncation method for large scale linear-time varying dynamical systems is presented. The last section is devoted to some numerical experiments.

2 The first approach: using an approximation of the matrix exponential

In this section, we give a new approach for computing approximate solutions to large differential equations (1). The expression of the exact solution as

[TABLE]

suggests the idea of computing $X(t)$ by approximating the factor $e^{(t-\tau)A}B$ and then using a quadrature method to compute the desired approximate solution.

As computing the exponential of a small matrix is straightforward , this is not the case for large scale problems, as $e^{(t-\tau)A}$ could be dense even though $A$ is sparse. However, in our problem, the computation of $e^{(t-\tau)A}$ is not needed as we will rather consider the product $e^{(t-\tau)A}\,B$, for which approximations via projection methods onto block or extended block Krylov subspaces are well suited.

Krylov subspace projection methods generate a sequence of nested subspaces (Krylov or extended Krylov subspaces). Let $\mathcal{V}_{m}=[V_{1},\ldots,V_{m}]$ be the orthogonal matrix whose columns form an orthonormal basis of the subspace $\mathit{K}_{m}$, Following [21, 22, 27], an approximation to $Z=e^{(t-\tau)A}\,B$ can be obtained as

[TABLE]

where $\mathcal{T}_{m}=\mathcal{V}^{T}_{m}A\mathcal{V}_{m}$. Therefore, the term appearing in the integral expression (8) can be approximated as

[TABLE]

If for simplicity, we assume $X_{t_{0}}=0$, an approximation to the solution of the differential Lyapunov equation (8) can be expressed as

[TABLE]

where

[TABLE]

and ${\widetilde{G}}_{m}(\tau)=e^{(t-\tau)\mathcal{T}_{m}}B_{m}$.

The next result shows that the matrix function $G_{m}$ is the solution of a low-order differential Lyapunov matrix equation.

Theorem 2**.**

Let $G_{m}(t)$ be the matrix function defined by (12), then it satisfies the following low-order differential Lyapunov matrix equation

[TABLE]

Proof.

The proof can be easily derived from the expression (12) and the result of Theorem 1. ∎

As a consequence, intruducing the residual $R_{m}(t)=\displaystyle{\dot{X}}_{m}(t)-A\,X_{m}-X_{m}\,A^{T}-BB^{T}$ associated to the approximation $X_{m}$, we have the following relation

[TABLE]

which shows that the residual satisfies a Petrov-Galerkin condition.

As mentioned earlier, once ${\widetilde{G}}_{m}(\tau)$ is computed, we use a quadrature method to approximate the integral (12) in order to approximate $G_{m}(t)$.

We now briefly discuss some practical aspects of the computation of $e^{(t-\tau)\mathcal{T}_{m}}B_{m}$ where $B_{m}={\cal V}_{m}^{T}B$, when $m$ is small and $\mathcal{T}_{m}$ is a an upper block Hessenberg matrix.

In the last decade, many approximation techniques such as the use of partial fraction expansions or Padé approximation have been proposed, see for example [9, 22]. However, it was remarked that a good way for evaluating the exponential of matrix times by a vector by using rational approximation to the exponential function. One of the main advantages of rational approximations as compared to polynomial approximations is the better stability of their integration schemes. Let us consider the rational function

[TABLE]

where the $\theta_{i}$’s are the poles of the rational function $F$. Then, the approximation to ${\widetilde{G}}_{m}(\tau)=e^{(t-\tau)\mathcal{T}_{m}}$ is given by

[TABLE]

One of the possible choices for the rational function $F$ is based on Chebychev approximation of the function $e^{x}$ on $[0,\,\infty[$, see [22]. We notice that for small values of $m$, one can also directly compute the matrix exponential $e^{(t-\tau){\mathcal{T}}_{m}}$ by using the well-known ’scaling and squaring method for the matrix exponential’ method, [13]. This method was associated to a Padé approximation and is implemented in the expm Matlab routine.

From now on, we assume that the basis formed by the orthonormal columns of ${\cal V}_{m}$ is obtained by applying the block Arnoldi or the extended block Arnoldi algorithm to the pair $(A,B)$.

The computation of $X_{m}(t)$ (and of $R_{m}(t)$) becomes expensive as $m$ increases. So, in order to stop the iterations, one has to test if $\parallel R_{m}\parallel<\epsilon$ without having to compute extra products involving the matrix $A$. The next result shows how to compute the residual norm of $R_{m}(t)$ without forming the approximation $X_{m}(t)$ which is computed in a factored form only when convergence is achieved.

Theorem 3**.**

Let $X_{m}(t)={\cal V}_{m}G_{m}(t){\cal V}_{m}^{T}$ be the approximation obtained at step $m$ by the block (or extended block) Arnoldi method. Then the residual $R_{m}(t)$ satisfies

[TABLE]

where $\bar{G}_{m}$ is the $d\times md$ matrix corresponding to the last $d$ rows of $G_{m}$ where $d=s$ when using the block Arnoldi and $d=2s$ for the extended block Arnoldi.

Proof.

The proof of this theorem comes directly from (11) and the fact that $G_{m}$ solves the low dimensional problem (13). ∎

The result of Theorem 3 is very important in practice, as it allows us to stop the iterations when convergence is achieved without computing the approximate solution $X_{m}(t)$.

The following result shows that the approximation $X_{m}$ is an exact solution of a perturbed differential Lyapunov equation.

Theorem 4**.**

Let $X_{m}(t)$ be the approximate solution given by (11). Then we have

[TABLE]

where $F_{m}=V_{m}\,T_{m+1,m}^{T}\,V_{m+1}^{T}$.

Proof.

The proof is easily obtained from (13) and the expression (11) of the approximate solution $X_{m}(t)$.

∎

Remark 1**.**

The solution $X_{m}(t)$ can be given as a product of two low rank matrices. Consider the eigen-decomposition of the symmetric and positive matrix $md\times md$ $G_{m}(t)=U\,D\,U^{T}$ where $D$ is the diagonal matrix of the eigenvalues of $G_{m}(t)$ sorted in decreasing order and $d=s$ for the block Arnoldi or $d=2s$ for the extended block Arnoldi. Let $U_{l}$ be the $md\times l$ matrix of the first $l$ columns of $U$ corresponding to the $l$ eigenvalues of magnitude greater than some tolerance $dtol$. We obtain the truncated eigen-decomposition $G_{m}(t)\approx U_{l}\,D_{l}\,U_{l}^{T}$ where $D_{l}={\rm diag}[\lambda_{1},\ldots,\lambda_{l}]$. Setting ${\widetilde{Z}}_{m}(t)={\cal V}_{m}\,U_{l}\,D_{l}^{1/2}$, it follows that

[TABLE]

Therefore, one has to compute and to store only the matrix ${\widetilde{Z}}_{m}(t)$ which is usually the required factor in some control problems such as in the balanced truncation method for model reduction in large scale dynamical systems. This possibility is very important for storage limitations in the large scale problems.

The next result states that the error matrix $X(t)-X_{m}(t)$ satisfies a differential Lyapunov matrix equation.

Theorem 5**.**

Let $X(t)$ be the exact solution of (1) and let $X_{m}(t)$ be the approximate solution obtained at step $m$. The error $E_{m}(t)=X(t)-X_{m}(t)$ satisfies the following equation

[TABLE]

and

[TABLE]

where $E_{m,0}=E_{m}(0)$.

Proof.

The result is easily obtained by subtracting the residual equation from the initial differential Lyapunov equation (1). ∎

Next, we give an upper bound for the norm of the error in the case where $A$ is a stable matrix.

Theorem 6**.**

Assume that $A$ is a stable matrix and $X(t_{0})=X_{m}(t_{0})$. Then we have the following upper bound

[TABLE]

where $\mu_{2}(A)=\displaystyle\frac{1}{2}\lambda_{max}(A+A^{T})<0$ is the 2-logarithmic norm and $\parallel\bar{G}_{m}\parallel_{\infty}=\displaystyle\max_{\tau\in[t_{0},\,t]}\parallel\bar{G}_{m}(\tau)\parallel$. The matrix $\bar{G}_{m}$ is the $d\times md$ matrix corresponding to the last $d$ rows of $G_{m}$ where $d=s$ when using the block Arnoldi and $d=2s$ for the extended block Arnoldi.

Proof.

We first remind that if $A$ is a stable matrix, then the logarithmic norm provides the following bound $\parallel e^{tA}\parallel\leq e^{\mu_{2}(A)t}$. Therefore, using the expression (19), we obtain the following relation

[TABLE]

Therefore, using (15) and the fact that $\parallel e^{(t-\tau)A}\parallel\leq e^{(t-\tau)\mu_{2}(A)}$, we get

[TABLE]

which gives the desired result.

∎

Notice that if $\parallel T_{m+1,m}\parallel$ is close to zero, which is the case when $m$ is close to the degree of the minimal polynomial of $A$ for $B$, then Theorem 6 shows that the error $E_{m}(t)$ tends to zero.

Next, we give another error bound for the norm of the error for every matrix $A$.

Theorem 7**.**

Let $X(t)$ be the exact solution to (1) and let $X_{m}(t)$ be the approximate solution obtained at step $m$. Then we have

[TABLE]

where $\mu_{2}(A)=\lambda_{max}((A+A^{T})/2)$, $Z(\tau)=e^{(t-\tau)A}B$ and $Z_{m}(\tau)={\cal V}_{m}e^{(t-\tau)\mathcal{T}_{m}}B_{m}$ with $B_{m}={\cal V}_{m}^{T}B$.

Proof.

From the expressions of $X(t)$ and $X_{m}(t)$, we have

[TABLE]

Therefore, using the fact that $\mu_{2}(\mathcal{T}_{m})=\lambda_{max}((\mathcal{T}_{m}+\mathcal{T}_{m}^{T})/2)\leq\lambda_{max}((A+A^{T})/2)=\mu_{2}(A)$, where $\mathcal{T}_{m}=\mathcal{V}_{m}^{T}A\mathcal{V}_{m}$, it follows that

[TABLE]

∎

When using a block Krylov subspace method such as the block Arnoldi method, then one can generalize to the block case the results already stated in many papers; see [7, 9, 12, 22]. In particular, we can easily generalize the result given in [22] for the case $s=1$ to the case $s>1$. In this case, we have the following upper bound.

[TABLE]

where $\rho=\|A\|$

The rupper bound (21) could be used in Theorem 7 to obtain a new upper bound for the norm of the error. In that case, we obtain the following upper bound

[TABLE]

We summarize the steps of our proposed first approach (using the extended block Arnoldi) in the following algorithm

3 A second approach: Projecting and solving with BDF

3.1 Low-rank approximate solutions via BDF

In this section, we show how to obtain low rank approximate solutions to the differential Lyapunov equation (1) by projecting directly the initial problem onto small block Krylov or extended block Krylov subspaces.

We first apply the block Arnoldi algorithm (or the extended block Arnoldi) to the pair $(A,B)$ to get the matrices ${\cal V}_{m}$ and ${\cal T}_{m}={\cal V}^{T}_{m}A{\cal V}_{m}$. Let $X_{m}(t)$ be the desired low rank approximate solution given as

[TABLE]

satisfying the Petrov-Galerkin orthogonality condition

[TABLE]

where $R_{m}(t)$ is the residual $R_{m}(t)=\displaystyle{\dot{X}}_{m}(t)-A\,X_{m}(t)-X_{m}(t)\,A^{T}-BB^{T}$. Then, from (23) and (24), we obtain the low dimensional differential Lyapunov equation

[TABLE]

with ${\cal T}_{m}={\cal V}^{T}_{m}A{\cal V}_{m}$ and ${B}_{m}={\cal V}_{m}^{T}\,B$. The obtained low dimensional differential Lyapunov equation (25) is the same as the one given by (13). For this second approach, we have to solve the latter low dimensional differential Lyapunov equation by some integration method such as the well known Backward Differentiation Formula (BDF).

Notice that we can also compute the norm of the residual without computing the approximation $X_{m}(t)$ which is also given, when convergence is achieved, in a factored form as in (17). The norm of the residual is given as

[TABLE]

where $\bar{Y}_{m}$ is the $d\times md$ matrix corresponding to the last $d$ rows of $Y_{m}$ where $d=s$ when using the block Arnoldi and $d=2s$ for the extended block Arnoldi.

3.2 BDF for solving the low order differential Lyapunov equation (25)

In this subsection, we will apply the Backward Differentiation Formula (BDF) method for solving, at each step $m$ of the block (or extended) block Arnoldi process, the low dimensional differential Lyapunov matrix equation (25). We notice that BDF is especially used for the solution of stiff differential equations.

At each time $t_{k}$, let $Y_{m,k}$ of the approximation of $Y_{m}(t_{k})$, where $Y_{m}$ is a solution of (25). Then, the new approximation $Y_{m,k+1}$ of $Y_{m}(t_{k+1})$ obtained at step $k+1$ by BDF is defined by the implicit relation

[TABLE]

where $h_{k}=t_{k+1}-t_{k}$ is the step size, $\alpha_{i}$ and $\beta_{i}$ are the coefficients of the BDF method as listed in Table 1 and ${\mathcal{F}}(X)$ is given by

[TABLE]

The approximate $Y_{m,k+1}$ solves the following matrix equation

[TABLE]

which can be written as the following Lyapunov matrix equation

[TABLE]

We assume that at each time $t_{k}$, the approximation $Y_{m,k}$ is factorized as a low rank product $Y_{m,k}\approx Z_{m,k}{Z_{m,k}}^{T}$, where $Z_{m,k}\in\mathbb{R}^{n\times m_{k}}$, with $m_{k}\ll n$. In that case, the coefficient matrices appearing in (28) are given by

[TABLE]

The Lyapunov matrix equation (28) can be solved by applying direct methods based on Schur decomposition such as the Bartels-Stewart algorithm [3, 11]. We notice that for large problems, many Krylov subspace type methods have been proposed to solve (28); [8, 14, 15, 16, 17, 25, 21].

Remark 2**.**

The main difference between Approach 1 and Approach 2 is the fact that in the first case, we compute an approximation of an integral using a quadrature formulae while in the second case, we have to solve a low dimensional differential Lyapunov equation using the BDF method. Mathematically, the two approaches are equivalent and they differ only in the way of computing numerically the low-order approximations: $G_{m}$ in the first approach and $Y_{m}$ in the second one.

‘

We summarize the steps of our proposed first approach (using the extended block Arnoldi) in the following algorithm

4 Application: Balanced truncation for linear time-varying dynamical systems

In this section, we assume that the coefficient matrices $A$ and $B$ are time-dependent. It is the case for example when we are dealing with Multi-Input Multi-Output (MIMO) linear-time varying (LTV) dynamical systems

[TABLE]

where $x(t)\in\mathbb{R}^{n}$ is the state vector, $u(t)\in\mathbb{R}^{p}$ is the control and $y(t)\in\mathbb{R}^{p}$ is the output. The matrices $A(t)\in\mathbb{R}^{n\times n}$, $B(t)\in\mathbb{R}^{n\times p}$ and $C(t)\in\mathbb{R}^{p\times n}$ are assumed to be continuous and bounded for all $t\in[t_{0},\,T_{f}]$.

The LTV dynamical system (29) can also be denoted as

[TABLE]

In many applications, such as circuit simulation, or time dependent PDE control problems, the dimension $n$ of ${\Sigma}$ is quite large, while the number of inputs and outputs is small $p\ll n$. In these large-scale settings, the system dimension makes the computation infeasible due to memory, time limitations and ill-conditioning. To overcome these drawbacks, one approach consists in reducing the model. The goal is to produce a low order system that has similar response characteristics as the original system with lower storage requirements and evaluation time.

The reduced order dynamical system can be expressed as follows

[TABLE]

where $x_{m}\in\mathbb{R}^{m}$, $y_{m}\in\mathbb{R}^{p}$, $A_{m}\in\mathbb{R}^{m\times m}$, $B\in\mathbb{R}^{m\times p}$ and $C_{m}\in\mathbb{R}^{p\times m}$ with $m\ll n$. The reduced dynamical system (31) is also represented as

[TABLE]

The reduced order dynamical system should be constructed in order that

• •

The output $y_{m}(t)$ of the reduced system approaches the output $y(t)$ of the original system.

• •

Some properties of the original system such as passivity and stability are preserved.

• •

The computation methods are steady and efficient.

One of the well known methods for constructing such reduced-order dynamical systems is the balanced truncation method for LTV systems [23, 24, 26]; see also [10, 19, 20] for the linear time-independent case. This method requires the LTV controllability and observability Gramians $P(t)$ and $Q(t)$ defined as the solutions of the differential Lyapunov matrix equations

[TABLE]

and

[TABLE]

Using the formulae (3), the differential Lyapunov equation (33) has the unique symmetric and positive solution $P(t)$ given by

[TABLE]

where the transition matrix $\Phi_{A}(t,\tau)$ is the unique solution of the problem

[TABLE]

The observability Gramian is given by

[TABLE]

The two LTV controllability and observability Gramians $P(t)$ and $Q(t)$ are then used to construct a new balanced system such that $\tilde{P}(t)=\tilde{Q}(t)=diag(\sigma_{1}(t),\ldots,\sigma_{n}(t))$ where the Hankel singular values are given by $\sigma_{i}(t)=\sqrt{\lambda_{i}(P(t)Q(t)}$, $i=1,\ldots,n$ and order in decreasing order.

The concept of balancing is to a transform the original LTV system to an equivalent one in which the states that are difficult to reach are also difficult to observe, which is finding an equivalent new LTV system such that the new Gramians ${\widetilde{P}}$ and ${\widetilde{Q}}$ are such that

[TABLE]

where $\sigma_{i}$ is the $i$-th Hankel singular value of the LTV system; i.e.

[TABLE]

Consider the Cholesky decompositions of the Gramians $P$ and $Q$:

[TABLE]

and consider also the singular value decomposition of $L_{c}(t)^{T}L_{o}(t)$ as

[TABLE]

where $Z(t)$ and $Y(t)$ are unitary $n\times n$ matrices and $\Sigma$ is a diagonal matrix containing the singular values. The balanced truncation consists in determining a reduced order model by truncating the states corresponding to the small Hankel singular values. Under certain conditions stated in [24], one can construct the low order model ${\Sigma}_{m}(t)$ as follows: We set

[TABLE]

where $\Sigma_{m}(t)=diag(\sigma_{1}(t),\ldots,\sigma_{m}(t))$; $Z_{m}(t)$ and $Y_{m}(t)$ correspond to the leading $m$ columns of the matrices $Z(t)$ and $Y(t)$ given by the singular value decomposition (36). The matrices of the reduced LTV system

[TABLE]

The use of Cholesky factors in the Gramians $P(t)$ and $Q(t)$ is not applicable for large-scale problems. Instead, one can compute low rank approximations of $P(t)$ an $Q(t)$ as given by (17) and use them to construct an approximate balanced truncation model.

As $A$, $B$ and $C$ are time-dependent, the direct application of the two approaches we developed is too expensive. Instead, we can apply directly an integration method such as BDF to the differential Lyapunov matrix equations (33) and (34). Then, at each iteration of the BDF method, we obtain a large Lyapunov matrix equation that can be numerically solved by using the extended block Arnoldi algorithm.

Consider the differential matrix equation (33), then, at each iteration of the BDF method, the approximation $P_{k+1}$ of $P(t_{k+1})$ where $P$ is the exact solution of (33), is given by the implicit relation

[TABLE]

where $h_{k}=t_{k+1}-t_{k}$ is the step size, $\alpha_{i}$ and $\beta_{i}$ are the coefficients of the BDF method as listed in Table 1 and ${\mathcal{G}}(X)$ is given by

[TABLE]

The approximate solution $P_{k+1}$ solves the following matrix equation

[TABLE]

which can be written as the following continuous-time algebraic Riccati equation

[TABLE]

Assuming that at each timestep, $P_{k}$ can be approximated as a product of low rank factors $P_{k}\approx\tilde{Z}_{k}\tilde{Z}_{k}^{T}$, $\tilde{Z}_{k}\in\mathbb{R}^{n\times m_{k}}$, with $m_{k}\ll n$, the coefficient matrices are given by

[TABLE]

A good way for solving the Lyapunov matrix equation (40) is by using the block or extended block Arnoldi algorithm applied to the pair $(\mathcal{A}_{k},\mathcal{B}_{k})$. This allows us to obtain low rank approximate solutions in factored forms. The procedure is as follows: applying for example the block Arnoldi to the pair $(\mathcal{A}_{k},\mathcal{B}_{k})$ we get, at step $m$ of the Arnoldi process, an orthonormal basis of the extended block Krylov subspace formed by the columns of the matrices: $\{V_{1,k},\ldots,V_{m,k}\}$ and also a block upper Hessenberg matrix $\mathbb{H}_{m,k}$. Let $\mathbb{V}_{m,k}=[V_{1,k},\ldots,V_{m,k}]$ and $\mathbb{H}_{m,k}=\mathbb{V}_{m,k}^{T}\mathcal{A}_{k}\mathbb{V}_{m,k}$. Then the obtained low rank approximate solution to the solution $P_{k+1}$ of (40) is given as $P_{m,k}=\mathbb{V}_{m,k}\mathbb{Y}_{m,k}\mathbb{V}_{m,k}^{T}$ where $\mathbb{Y}_{m,k}$ is solution of the following low order Lyapunov equation

[TABLE]

where $\mathcal{\tilde{B}}_{k}=\mathbb{V}_{m,k}^{T}\,\mathcal{B}_{k}$. As stated in Remark 1, the approximate solution can be given in a factored form.

5 Numerical examples

In this section, we compare the two approaches presented in this paper. The exponential approach (EBA-exp) summarized in Algorithm 1, which is based on the approximation of the solution to (1) applying a quadrature method to compute the projected exponential form solution (12). We used a scaling and squaring strategy, implemented in the MATLAB expm function; see [12, 18] for more details. The second method (Algorithm 2) is based on the BDF integration method applied to the projected Lyapunov equation (25). The basis of the projection subspaces were generated by the extended block Arnoldi algorithm for both methods. All the experiments were performed on a laptop with an Intel Core i7 processor and 8GB of RAM. The algorithms were coded in Matlab R2014b.

Example 1. The matrix $A$ was obtained from the 5-point discretization of the operators

[TABLE]

on the unit square $[0,1]\times[0,1]$ with homogeneous Dirichlet boundary conditions. The number of inner grid points in each direction is $n_{0}=$ and the dimension of the matrix $A$ was $n=n_{0}^{2}=$. Here we set $f_{1}(x,y)=10xy$, $f_{2}(x,y)=e^{x^{2}y}$, $f_{3}(x,y)=100y$, $f_{4}(x,y)={x^{2}y}$ , $g_{1}(x,y)=20y$ and $g_{2}(x,y)=x\,y$. The time interval considered was $[0,\,2]$ and the initial condition $X_{0}=X(0)$ was choosen as the low rank product $X_{0}=Z_{0}Z_{0}^{T}$, where $Z_{0}=0_{n\times 2}$. For both methods, we used projections onto the Extended Block Krylov subspaces

[TABLE]

and the tolerance was set to $10^{-10}$ for the stop test on the residual. For the EBA-BDF method, we used a 2-step BDF scheme with a constant timestep $h$. The entries of the matrix $B$ were random values uniformly distributed on the interval $[0,\,1]$ and the number of the columns in $B$ was $s=2$.

literature. In order to check if our approaches produce reliable results, we began comparing our results to the one given by Matlab’s ode23s solver which is designed for stiff differential equations. This was done by vectorizing our DLE, stacking the columns of $X$ one on top of each other. This method, based on Rosenbrock integration scheme, is not suited to large-scale problems. Due to the memory limitation of our computer when running the ode23s routine, we chose a size of $100\times 100$ for the matrix $A$.

In Figure 1, we compared the component $X_{11}$ of the solution obtained by the methods tested in this section, to the solution provided by the ode23s method from Matlab, on the time interval $[0,\,2]$, for $size(A)=100\times 100$ and a constant timestep $h=10^{-3}$.

We observe that all the considered methods give similar results in terms of accuracy. The relative error norms $\displaystyle{\frac{\|X_{EBA-exp}(t_{f})-X_{ode23s}(t_{f})\|}{\|X_{ode23s(t_{f})}\|}}$ and $\displaystyle{\frac{\|X_{EBA-BDF(2)}(t_{f})-X_{ode23s}(t_{f})\|}{\|X_{ode23s(t_{f})}\|}}$ at final time $t_{f}=2$ were equal to $1.8\times 10^{-10}$ and $9.1\times 10^{-11}$ respectively. The runtimes were respectively 0.59s, 5.1s for the EBA-exp and EBA-BDF(2) methods and 1001s for the ode23s routine.

In Table 2, we give the obtained runtimes in seconds, for the resolution of Equation (1) for $t\in[0,\,2]$, with a timestep $h=0.001$ and the Frobenius norm of the residual at the final time.

The results in Table 2 illustrate that the EBA-exp method clearly outperforms the EBA-BDF(2) method in terms of computation time even though both methods are equally accurate. In Figure 2, we featured the norm of the residual at final time $t=2$ for both EBA-exp and EBA-BDF(2) methods for size($A$)$=6400\times 6400$ in function of the number $m$ of extended Arnoldi iterations. We observe that the plots coincide for both methods.

Example 2. This example comes from the autonomous linear-quadratic optimal control problem of one dimensional heat flow

[TABLE]

Using a standard finite element approach based on the first order B-splines, we obtain the following ordinary differential equation

[TABLE]

where the matrices $M$ and $K$ are given by:

[TABLE]

Using the semi-implicit Euler method, we get the following discrete dynamical system

[TABLE]

We set $A=(M-\Delta tK)^{-1}\,M$ and $B=\Delta t\,(M-\Delta tK)^{-1}\,F$. The entries of the $n\times s$ matrix $F$ and the $s\times n$ matrix $C$ were random values uniformly distributed on $[0,\,1]$. In our experiments we used $n=$, $s=2$, $\Delta t=0.01$ and $\alpha=0.05$.

In Table 3, we give the obtained runtimes in seconds, for the resolution of Equation (1) for $t\in[0,\,2]$, with a timestep $h=0.001$ and the Frobenius norm of the residual at the final time.

The figures in Table 3 illustrate the gain of speed provided by the EBA-exp method. Again, both methods performed similarly in terms of accuracy. In figure 3, we considered the case size($A$)$=100\times 100$ and plotted the upper bound of the error norms as stated in Formula (20) at the final time $T_{f}$ against the computed norm of the errors, taking the solution given by the integral formula (8) as a reference, in function of the number $m$ of Arnoldi iterations for the EBA-exp method.

Example 3 In this last example, we applied the EBA-BDF(1) method to the well-known problem Optimal Cooling of Steel Profiles. The matrices were extracted from the IMTEK collection 111https://portal.uni-freiburg.de/imteksimulation/downloads/benchmark. We compared the EBA-BDF(2) method to the EBA-exp method for problem sizes $n=1357$ and $n=5177$, on the time interval $[0\,,1000]$. The initial value $X_{0}$ was chosen as $X_{0}=0$ and the timestep was set to $h=0.01$. The tolerance for the Arnoldi stop test was set to $10^{-7}$ for both methods and the projected low dimensional Lyapunov equations were numerically solved by the solver (lyap from Matlab) at each iteration of the extended block Arnoldi algorithm for the EBA-BDF(2) method.

In Table 4, we listed the obtained runtimes which again showed the advantage of the EBA-exp method in terms of execution time and similar accuracy for both methods.

6 Conclusion

We presented in the present paper two new approaches for computing approximate solutions to large scale differential Lyapunov matrix equations. The first one comes naturally from the exponential expression of the exact solution and the use of approximation techniques of the exponential of a matrix times a block of vectors. The second approach is obtained by first projecting the initial problem onto a block Krylov (or extended Krylov) subspace, obtain a low dimensional differential Lyapunov equation which is solved by using the well known BDF integration method. We gave some theoretical results such as the exact expression of the residual norm and also upper bounds for the norm of the errors. An application in model reduction for linear time-varying dynamical systems is also given. Numerical experiments show that both methods are promising for large-scale problems, with a clear advantage for the EBA-exp method in terms of computation time.