Order reduction methods for solving large-scale differential matrix Riccati equations

TL;DR

This paper introduces a new order reduction method using rational Krylov subspaces for efficiently solving large-scale symmetric differential matrix Riccati equations, significantly reducing computational and memory costs.

Contribution

The paper presents a novel reduction process onto rational Krylov subspaces combined with a two-phase strategy for improved accuracy and efficiency in solving large-scale differential Riccati equations.

Findings

Significant computational savings over existing methods

Ability to solve larger problems with lower memory requirements

Effective numerical results on benchmark problems

Abstract

We consider the numerical solution of large-scale symmetric differential matrix Riccati equations. Under certain hypotheses on the data, reduced order methods have recently arisen as a promising class of solution strategies, by forming low-rank approximations to the sought after solution at selected timesteps. We show that great computational and memory savings are obtained by a reduction process onto rational Krylov subspaces, as opposed to current approaches. By specifically addressing the solution of the reduced differential equation and reliable stopping criteria, we are able to obtain accurate final approximations at low memory and computational requirements. This is obtained by employing a two-phase strategy that separately enhances the accuracy of the algebraic approximation and the time integration. The new method allows us to numerically solve much larger problems than in the current literature. Numerical experiments on benchmark problems illustrate the effectiveness of the procedure with respect to existing solvers.

Tables

Table 1. Table 1: Relevant information concerning the experimental data

Name	$n$	$p$/$s$/$q$	${\Vert A\Vert }_F$	${\Vert B\Vert }_F$	${\Vert C\Vert }_F$	${\Vert Z\Vert }_F$	${\Vert E\Vert }_F$
sym2d	$640000$	$5/1/1$	$3.6\cdot {10}^3$	$8.0\cdot {10}^2$	$1.8\cdot {10}^3$	$8.0\cdot {10}^2$	$8\cdot {10}^2$
nsym3d	$64000$	6/1/3	$2.0\cdot {10}^3$	$2.5\cdot {10}^2$	$6.2\cdot {10}^2$	$2.8\cdot {10}^2$	$2.5\cdot {10}^2$
Name	$n$	$p$/$s$/$q$	${\Vert \hat{A}\Vert }_F$	${\Vert \hat{B}\Vert }_F$	${\Vert \hat{C}\Vert }_F$	${\Vert Z\Vert }_F$	${\Vert \hat{E}\Vert }_F$
chip	$20082$	5/1/1	$2.2\cdot {10}^6$	$1.7\cdot {10}^2$	$3.3\cdot {10}^4$	$1.0\cdot {10}^2$	$2\cdot {10}^{-4}$
flow	$9669$	5/1/1	$4.5\cdot {10}^6$	$2.0\cdot {10}^4$	$1.2\cdot {10}^3$	$-$	$6.8\cdot {10}^0$
rail	$79841$	7/6/1	$7\cdot {10}^{-3}$	$1\cdot {10}^{-7}$	$6.2\cdot {10}^0$	$1.9\cdot {10}^2$	$8\cdot {10}^{-4}$

Table 2. Table 2: A breakdown of the computational time for the considered methods for the first two datasets.

		System	Orthogonalisation	Integration
Data	Method	solves (s)	steps (s)	steps (s)
sym2d	rksm-dre	6.1	6.9	0.4
	eksm-dre	8.6 (2.7)	12.1	1.3
nsym3d	rksm-dre	38.3	0.9	0.8
	eksm-dre	48.6 (43.5)	1.6	4.0

Table 3. Table 3: Data information for comparisons between projection-based methods and m.e.s.s.

Name	$n$	$p$/$s$/$q$	${\Vert A\Vert }_F$	${\Vert B\Vert }_F$	${\Vert C\Vert }_F$	${\Vert Z\Vert }_F$	${\Vert E\Vert }_F$
sym2d	$40000$	5/1/1	$1.3\cdot {10}^3$	$3.0\cdot {10}^2$	$6.7\cdot {10}^2$	$3.0\cdot {10}^2$	$2\cdot {10}^2$
nsym3d	$8000$	6/1/3	$6.1\cdot {10}^2$	$7.7\cdot {10}^1$	$1.9\cdot {10}^2$	$8.3\cdot {10}^1$	$2.8\cdot {10}^2$
Name	$n$	$p$/$s$/$q$	${\Vert \hat{A}\Vert }_F$	${\Vert \hat{B}\Vert }_F$	${\Vert \hat{C}\Vert }_F$	${\Vert Z\Vert }_F$	${\Vert \hat{E}\Vert }_F$
flow	$9669$	5/1/1	$4.5\cdot {10}^6$	$2.0\cdot {10}^4$	$1.2\cdot {10}^3$	$0$	$6.8\cdot {10}^0$
rail	$20209$	7/6/1	$4\cdot {10}^{-3}$	$2.1\cdot {10}^{-7}$	$6.2\cdot {10}^0$	$1.9\cdot {10}^2$	$2\cdot {10}^{-4}$

Table 4. Table 4: sym2d : Storage and computational time comparison of rksm-dre , eksm-dre and m.e.s.s. . Reduction phase performed with bdf (1,10), refinement phase with bdf (2,100). In m.e.s.s. only the approximate solution at the final time is stored, with no solutions at intermediate time instances returned.

	# $n$-long	Min/Max	Reduction	Refine	Tot CPU
Method	Vecs	rank	phase(s)	phase(s)	time(s)
rksm-dre	66	23/43	1.5	0.19	1.7
eksm-dre	144	23/43	1.9	2.8	4.7
m.e.s.s.-bdf(1,10)	988	58/75			319.9
m.e.s.s.-bdf(2,100)	1032	58/86			4005.4

Table 5. Table 5: nsym3d : Storage and computational time comparison of rksm-dre , eksm-dre and m.e.s.s. . Reduction phase performed with bdf (1,10), refinement phase with bdf (2,100). In m.e.s.s. only the approximate solution at the final time is stored, with no solutions at intermediate time instances returned.

	# $n$-long	Min/Max	Reduction	Refine	Tot CPU
Method	Vecs	rank	phase(s)	phase(s)	time(s)
rksm-dre	108	36/66	4.5	4.6	9.1
eksm-dre	196	36/66	2.8	5.7	8.5
m.e.s.s.-bdf(1,10)	1116	71/90			431.0
m.e.s.s.-bdf(2,100)	1152	67/94			4965.0

Table 6. Table 6: flow : Storage and computational time comparison of rksm-dre , eksm-dre and m.e.s.s. . Reduction phase performed with bdf (1,10), refinement phase with bdf (2,100). In m.e.s.s. only the approximate solution at the final time is stored, with no solutions at intermediate time instances returned.

	# $n$-long	Min/Max	Reduction	Refine	Tot CPU
Method	Vecs	rank	phase(s)	phase(s)	time(s)
rksm-dre	186	95/100	12	4.6	16.6
eksm-dre	372	95/100	30.8	23.7	54.5
m.e.s.s.-bdf(1,10)	1280	87/106			431.7

Table 7. Table 7: rail : Storage and computational time comparison of rksm-dre , eksm-dre and m.e.s.s. . Reduction phase performed with bdf (1,10), refinement phase with bdf (2,100). In m.e.s.s. only the approximate solution at the final time is stored, with no solutions at intermediate time instances returned.

	# $n$-long	Min/Max	Reduction	Refine	Tot CPU
Method	Vecs	rank	phase(s)	phase(s)	time(s)
rksm-dre	168	153/160	6.4	3.3	9.7
eksm-dre	462	153/160	39.2	5.7	44.9
m.e.s.s.-bdf(1,10)	6345	151/158			705.3
m.e.s.s.-bdf(2,100)	4023	124/158			3396.5

Table 8. Table 8: Storage and computational time comparison of rksm-dre and split-add4(500) . Reduction phase performed with bdf (1,10), refinement phase with bdf (3,1000). In split-add4 only the approximate solution at the final time is stored.

		# $n$-long	Min/Max	Tot CPU
Data ($n$)	Method	Vecs	rank	time (s)
sym2d $\left({10}^4\right)$	rksm-dre	8	3/6	0.6
	split-add4(500)	28	3/7	34.9
nsym3d $\left(8\cdot {10}^3\right)$	rksm-dre	10	4/7	2.2
	split-add4(500)	36	3/9	37.9
sym2d $\left(9\cdot {10}^4\right)$	rksm-dre	6	3/4	1.2
	split-add4(500)	28	3/7	330.0
nsym3d $\left(2.7\cdot {10}^4\right)$	rksm-dre	10	4/8	10.1
	split-add4(500)	36	3/9	127.8

Table 9. Table 9: sym2d (of size 10 6 superscript 10 6 10^{6} ): Results with rksm-dre , using different refinement strategies. Reduction phase performed with bdf (1,10) and tolerance 10 − 7 superscript 10 7 10^{-7} .

Refinement	# $n$-long	Soln.	Reduction	Refinement	Tot CPU
Method	Vecs	rank	phase(s)	phase(s)	time(s)
bdf(2,100)	72	55	37.7	1.1	38.8
bdf(3,1000)	72	55	37.7	9.6	47.3
bdf(4,10000)	72	55	37.7	95.5	133.2
split-add4(500)	72	55	37.7	25.9	63.6
split-add8(500)	72	55	37.7	59.6	97.3
split-adapt8	72	55	37.7	160.7	198.4

Equations

$$\dot{X}(t)=A^{T}X(t)+X(t)A-X(t)B{B}^{T}X(t)+C^{T}C,X(0)=X_0,$$

$$\dot{x}(t)=Ax(t)+Bu(t),y(t)=Cx(t),x(0)=x_0$$

$$J(u)=x(t_f{)}^T{P}_{f}x(t_f)+{\int }_0^{t_f}\left(x(t{)}^T{C}^{T}Cx(t)+u(t{)}^{T}u(t)\right)dt.$$

$$\dot{P}(t)=A^{T}P(t)+P(t)A-P(t)B{B}^{T}P(t)+C^{T}C,P(t_f)=P_f.$$

$$X(t)=e^{t{A}^T}Z{Z}^T{e}^{tA}+{\int }_0^t{e}^{(t-s)A^T}\left(C^{T}C-X(s)B{B}^{T}X(s)\right)e^{(t-s)A}ds,$$

$$0=A^T{X}_{\infty }+X_{\infty }A-X_{\infty }B{B}^T{X}_{\infty }+C^{T}C.$$

$$\dot{\hat{x}}(t)=A_m\hat{x}(t)+B_{m}u(t),y(t)=C_m\hat{x}(t),\hat{x}(0)={\hat{x}}_0$$

$${\mathcal{K}}_m(A,N)$$

$${\mathcal{E}\mathcal{K}}_m(A,N)$$

$${\mathcal{R}\mathcal{K}}_m(A,N,\mathbf{s})$$

$$A^T{\mathcal{V}}_m={\mathcal{V}}_m{\mathcal{T}}_m^T+{\nu }_{m+1}{\tau }_m^T,$$

$${\mathcal{T}}_m={\mathcal{V}}_m^{T}A{\mathcal{V}}_m,B_m={\mathcal{V}}_m^{T}B,Z_m={\mathcal{V}}_m^{T}Z\text{mbox}and{C}_m=C{\mathcal{V}}_m.$$

$$\begin{array}{rl}{\dot{Y}}_m(t)& ={\mathcal{T}}_m^T{Y}_m(t)+Y_m(t){\mathcal{T}}_m-Y_m(t)B_m{B}_m^T{Y}_m(t)+C_m^T{C}_m\\ Y_m(0)& =Z_m{Z}_m^T,\end{array}$$

$$X_m(t)={\mathcal{V}}_m{Y}_m(t){\mathcal{V}}_m^T\approx X(t),t\in [0,t_f].$$

$${\mathcal{R}}_m\left(t\right):={\dot{X}}_m(t)-A^T{X}_m(t)-X_m(t)A+X_m(t)B{B}^T{X}_m(t)-C^{T}C$$

$${\mathcal{V}}_m^T{\mathcal{R}}_m\left(t\right){\mathcal{V}}_m=0,t\in [0,t_f],$$

$$\mathcal{V}_{m}^{T}\big{(}\mathcal{V}_{m}\dot{Y}_{m}(t)\mathcal{V}_{m}^{T}-A^{T}\mathcal{V}_{m}Y_{m}(t)\mathcal{V}_{m}^{T}-\mathcal{V}_{m}Y_{m}(t)\mathcal{V}_{m}^{T}A+\mathcal{V}_{m}Y_{m}(t)\mathcal{V}_{m}^{T}BB^{T}\mathcal{V}_{m}Y_{m}(t)\mathcal{V}_{m}^{T}-C^{T}C)\mathcal{V}_{m}=0,$$

$$\mathbf{F}\left(Y^{(k+1)}\right)={\mathcal{T}}_m^T{Y}^{(k+1)}+Y^{(k+1)}{\mathcal{T}}_m-Y^{(k+1)}{B}_m{B}_m^T{Y}^{(k+1)}+C_m^T{C}_m,$$

$$Y^{(k+1)}=\sum _{i=0}^{b-1}{\alpha }_i{Y}^{(k-i)}+h\beta \mathbf{F}(Y^{(k+1)}),$$

$$-Y^{(k+1)}+h\beta \left({\mathcal{T}}_m^T{Y}^{(k+1)}+Y^{(k+1)}{\mathcal{T}}_m-Y^{(k+1)}{B}_m{B}_m^T{Y}^{(k+1)}+C_m^T{C}_m\right)+\sum _{i=0}^{b-1}{\alpha }_i{Y}^{(k-i)}=0,$$

$${\hat{\mathcal{T}}}_m^T{Y}^{(k+1)}+Y^{(k+1)}{\hat{\mathcal{T}}}_m-Y^{(k+1)}{\hat{B}}_m{\hat{B}}_m^T{Y}^{(k+1)}+{\hat{Q}}_m=0.$$

$${\hat{\mathcal{T}}}_m=h\beta {\mathcal{T}}_m-\frac{1}{2}{I}_m,{\hat{B}}_m=\sqrt{h\beta }{B}_m,{\hat{Q}}_m=h\beta C_m^T{C}_m+\sum _{i=0}^{b-1}{\alpha }_i{Y}^{(k-i)}.$$

$${\partial }_{t}w=\Delta w-10x{w}_x-100y{w}_y,w{|}_{\partial \Omega }=0$$

$$\Vert {\mathcal{R}}_m(t){\Vert }_F^2=\Vert R_m^{(I)}(t){\Vert }_F^2+2\Vert R_m^{(O)}(t){\Vert }_F^2,$$

$$\begin{array}{rl}{\mathcal{R}}_m(t)={\mathcal{V}}_m{\dot{Y}}_m(t){\mathcal{V}}_m^T& -A^T{\mathcal{V}}_m{Y}_m(t){\mathcal{V}}_m^T-{\mathcal{V}}_m{Y}_m(t){\mathcal{V}}_m^{T}A\\ & +{\mathcal{V}}_m{Y}_m(t){\mathcal{V}}_m^{T}B{B}^T{\mathcal{V}}_m{Y}_m(t){\mathcal{V}}_m^T-C^{T}C.\end{array}$$

$$\begin{array}{rl}{\mathcal{R}}_m(t)={\mathcal{V}}_m{\dot{Y}}_m(t){\mathcal{V}}_m^T& -({\mathcal{V}}_m{\mathcal{T}}_m^T+{\nu }_{m+1}{\tau }_m^T)Y_m(t){\mathcal{V}}_m^T-{\mathcal{V}}_m{Y}_m(t)({\mathcal{T}}_m{\mathcal{V}}_m^T+{\tau }_m{\nu }_{m+1}^T)\\ & +{\mathcal{V}}_m{Y}_m(t){\mathcal{V}}_m^{T}B{B}^T{\mathcal{V}}_m{Y}_m(t){\mathcal{V}}_m^T-{\mathcal{V}}_m{C}_m^T{C}_m{\mathcal{V}}_m^T.\end{array}$$

$$\mathcal{J}_{m}(t)=\left[\begin{array}[]{c|c}\dot{Y}_{m}(t)-\mathcal{T}_{m}^{T}Y_{m}(t)-Y_{m}(t)\mathcal{T}_{m}+Y_{m}(t)B_{m}B_{m}^{T}Y_{m}(t)-C_{m}^{T}C_{m}&Y_{m}(t)\tau_{m}\\ \hline\cr\tau_{m}^{T}Y_{m}(t)&\mathbf{0}\end{array}\right].$$

$$\mathcal{J}_{m}(t)=\left[\begin{array}[]{c|c}R_{m}^{(I)}(t)&Y_{m}(t)\tau_{m}\\ \hline\cr\tau_{m}^{T}Y_{m}(t)&\mathbf{0}\end{array}\right].$$

$$\begin{array}{rl}\Vert {\mathcal{R}}_m(t){\Vert }_F^2& =\Vert {\mathcal{V}}_{m+1}{\mathcal{J}}_m(t){\mathcal{V}}_{m+1}^T{\Vert }_F^2=\Vert {\mathcal{J}}_m(t){\Vert }_F^2\\ & =\text{mbox}Tr\left(R_m^{(I)}(t{)}^T{R}_m^{(I)}(t)+2(Y_m(t){\tau }_m)({\tau }_m^T{Y}_m(t))\right),\end{array}$$

$${\dot{Y}}_m(t)-{\mathcal{T}}_m^T{Y}_m(t)-Y_m(t){\mathcal{T}}_m+Y_m(t)B_m{B}_m^T{Y}_m(t)-C_m{C}_m^T-R_m^{(I)}(t)=0.$$

$${\int }_0^{t_f}\Vert R_m^{(O)}(\gamma ){\Vert }_{F}d\gamma \approx \sum _{j=1}^{\ell }\frac{t_f}{\ell }\Vert R_m^{(O)}(t_j){\Vert }_F=:{\rho }_m$$

Full text

1 Introduction

We consider the solution of the continuous-time differential matrix Riccati equation (DRE in short) of the form

[TABLE]

in the unknown matrix $X(t)\in\mathbb{R}^{n\times n}$, where $X_{0}=ZZ^{T}$ and $t\in[0,t_{f}]$. Here, $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times s}$, $C\in\mathbb{R}^{p\times n}$ and $Z\in\mathbb{R}^{n\times q}$ are time invariant, and $s,p,q\ll n$. The matrix $A$ is assumed to be large, sparse and nonsingular, whereas $B$, $C$ and $Z$ have full rank. In particular, we consider low-rank DREs, where both matrices $C^{T}C$ and $X_{0}$ have very low rank compared to $n$. Even though the matrix $A$ is sparse, the solution $X(t)$ is typically dense and impossible to store when $n$ is large. Under the considered hypotheses, numerical evidence seems to indicate that $X(t)$ usually has rapidly decaying singular values, hence a low-rank approximation to $X(t)$ may be considered, see e.g., [48]. For completeness, we also refer the reader to [21, 20] for results on the existence of low-rank solutions for the algebraic Sylvester and Lyapunov (linear) equations.

The DRE plays a fundamental role in optimal control theory, filter design theory, model reduction problems, as well as in differential games [2, 7, 11, 13, 38]. Equations of the form Eq. 1 are crucial in the numerical treatment of the linear quadratic regulator (LQR) problem [2, 13, 29]: given the state equation

[TABLE]

consider the finite horizon case, where the finite time cost integral has the form

[TABLE]

The matrix $P_{f}$ is assumed to be symmetric and nonnegative definite. Assuming that the pair $(A,B)$ is stabilizable and the pair $(C,A)$ is detectable, the optimal input $\tilde{u}(t)$, minimizing Eq. 3, can be determined as $\tilde{u}(t)=-B^{T}P(t)\tilde{x}(t)$, and the optimal trajectory is subject to $\dot{\tilde{x}}=(A-BB^{T}P(t))x(t)$. The matrix $P(t)$ is the solution to the DRE

[TABLE]

Using a common practice, we can transform Eq. 4 into the initial value problem Eq. 1 via the change of variables $X(t_{f}-t)=P(t)$.

Under certain assumptions, the exact solution of Eq. 1 can be expressed in integral form as (see e.g., [28, Theorem 8])

[TABLE]

so that when $t\rightarrow\infty$ the DRE reaches a steady state solution satisfying the algebraic Riccati equation (ARE)

[TABLE]

In the framework of differential equations, the DRE is characterized by both fast and slow varying modes, hence it is classified as a stiff ordinary differential equation (ODE). The stiffness and the nonlinearity of the DRE are responsible for the difficulties in its numerical solution even on a small scale ($n<10^{3}$). Several stiff integrators have been investigated, including the matrix generalizations of implicit ODE solvers [15, 12], linearization methods [14] and more recently matrix versions of splitting methods [35, 46, 47]. These methods are feasible on a small scale but fail to be efficient when $n$ is large. In [41], iterative methods are implemented within the matrix generalization of standard implicit methods allowing for the computation of an approximate solution to the DRE when $n\gg 10^{3}$. These algorithms require the solution of a large algebraic Riccati equation at each timestep, which again raises big concerns as of storage and computational efforts.

A promising idea is to rely on a model order reduction strategy typically used in linear and nonlinear dynamical systems. In this setting, the original system is replaced with

[TABLE]

where $A_{m},B_{m}$ and $C_{m}$ are projections and restrictions of the original matrices onto a subspace of small dimension. The differential Riccati equation associated with this reduced order problem is solved, yielding an optimal corresponding cost function. This strategy allows for a natural low-rank approximation to the sought after DRE solution $X(t)$, obtained by interpolating the reduced order solution at selected time instances. One main feature is that a single space is used for all time snapshots, so that the approximate solutions can be kept in factored form with few memory allocations. We refer the reader to, e.g., [4] for a general presentation of algebraic reduction methods for linear dynamical systems, and to [43] for a detailed discussion motivating the reduction approach in the context of the algebraic Riccati equation.

A key ingredient for the success of the reduction methodology is the choice of the approximation space onto which the algebraic reduction is performed; [4] presents a comprehensive description of various space selections in the dynamical system setting. Following strategies already successfully adopted for the algebraic Riccati equations, the authors of [27] and [23] have independently used polynomial and extended Krylov subspaces as approximation space, respectively, in the differential setting. A major characteristic of these spaces is that their dimension can be expanded iteratively, so that if the determined approximate solution is not sufficiently accurate, the Krylov space can be enlarged and the process continued. Several questions remain open in the methods proposed in [27],[23]. On the one hand, it is well known that polynomial Krylov subspaces require a very large dimension to satisfactorily solve real application problems, thus destroying the reduction advantages. On the other hand, the multiple timestepping proposed in the method in [23] only provides an accurate approximation at $t=t_{f}$, except when $X_{0}=0$. For $X_{0}=ZZ^{T}\neq 0$ of low rank, memory requirements of the extended method grow significantly. These problems can be satisfactorily solved by using a general rational Krylov subspace, which is shown in various applications to be able to supply good spectral information on the involved matrices with much smaller dimensions than the polynomial and extended versions. Such gain has been experimentally reported in the literature in the solution of the algebraic Riccati equation. We show that great computational and memory savings can be obtained when projecting onto the fully rational Krylov subspace, and that with an appropriate implementation the extended Krylov subspace may also be competitive with certain data.

A related issue that has somehow been overlooked in the available literature is the expected final accuracy and thus the stopping criterion. Time dependence of the DRE makes the reduced problem trickier to handle than in the purely algebraic case; in particular, two intertwined issues arise: i) The accuracy of the approximate solution may vary considerably within the time interval $[0,t_{f}]$; ii) Throughout the reduction process the reduced ODE cannot be solved with high accuracy and, quite the opposite, low-order methods should be used to make the overall cost feasible. We analyze these difficulties in detail, and by exploiting the inherent structure of the reduced order model, we derive a two-phase strategy that first focuses on the reduction and then on the integration, in a way that is efficient for memory and CPU time usage, but also in terms of final expected accuracy.

We also discuss several algebraic properties of the approximate solution and its relation both with the solution $X(t)$ for $t\in[0,t_{f}]$, and with the steady state solution $X_{\infty}$. These results continue a matrix analysis started in [27], where positivity and monotonicity properties of the approximate solution obtained by certain reduction methods are explored.

The paper is organized as follows. In Section 2 we introduce reduction methods and discuss the use of Krylov subspace based strategies. Matrix-oriented BDF methods are recalled for the solution of the projected problem in Section 3. In Section 4 we devise a stopping criterion for the order reduction methods and illustrate its key role in the implementation. Section 5 is devoted to the analysis of matrix properties of the solution, as well as the reduced model, from a control theory perspective. Several numerical experiments are reported in Section 6, where the new methods are also compared with state-of-the-art procedures. Our conclusions are discussed in Section 7. Finally, in Appendix A and Appendix B we review some properties of the extended and rational Krylov subspaces.

Notation and definitions. Throughout the paper, the matrix $I_{n}$ will denote the $n\times n$ identity matrix. In terms of norms, $\|\cdot\|$ refers to any induced matrix norm, where in particular the Frobenius norm is denoted by $\|\cdot\|_{F}$. A matrix $A$ is stable (sometimes also called Hurwitz) if all its eigenvalues are contained in the left half open complex plane. A linear dynamical system, $\dot{x}=Ax,$ is called dissipative if the real matrix $A$ has its field of values contained in the left half open complex plane.

All reported experiments were performed using MATLAB 9.4 (R2018b) ([33]) on a MacBook Pro with 8-GB memory and a 2.3-GHz Intel core i5 processor.

2 Order reduction with Krylov-based subspaces

In this section, we review Krylov-based order reduction methods and show how they are applied to the DRE. Krylov subspaces that have been explored in the past years have the form

[TABLE]

where $N$ is a tall matrix associated with the given problem. In the rational subspace, ${\boldsymbol{s}}=\{s_{2},\ldots,s_{m}\}$ is a set of properly chosen real or complex shifts, whose computation can be performed a priori or dynamically during the generation of the subspace; we refer the reader to [44, 18] for more complete descriptions.

Krylov-based projection methods (in short generically denoted as ${\cal K}_{m}$) were first applied to ARE’s in [25] (polynomial spaces) and later improved in [24] (extended space) and [45] (rational spaces). The two rational spaces prove to be far superior to the polynomial Krylov space in most reduction strategies where they are applied in the literature, as long as solving linear systems at each iteration is feasible. The differential Riccati equation has been attacked in [23] with the extended space, and in [27] with the polynomial space; here we close the gap, as far as Krylov subspaces are concerned. In addition, we address several implementation issues to make the final method computationally reliable and, to the best of our knowledge, a great competitor among the available methods for large-scale DRE problems.

While for the algebraic Riccati equation $N=C^{T}$, in the differential context the starting matrix for generating these spaces is given by $N=[C^{T},Z]$, where $X_{0}=ZZ^{T}$. Both matrices $C$ and $Z$ play a crucial role in the closed-form DRE solution matrix and are thus included to generate the projection space. The idea of reduction methods is to first project the large DRE onto the smaller subspace $\mathcal{K}_{m}$, then solve the projected equation, and finally expand the solution back to the original space.

Let the columns of $\mathcal{V}_{m}\in\mathbb{R}^{n\times d}$ span the considered Krylov subspace. Then the following Arnoldi-type relation holds,

[TABLE]

where the actual values of $\nu_{m+1}\in\mathbb{R}^{n}$ and $\tau_{m}^{T}$ depend on the chosen subspace. Moreover, setting ${\cal V}_{m+1}=[\mathcal{V}_{m},\nu_{m+1}]$ we have that $\mathcal{K}_{m+1}={\rm range}({\cal V}_{m+1})$, which shows that Krylov subspaces are nested, that is $\mathcal{K}_{m}\subseteq\mathcal{K}_{m+1}$, resulting in a dimension increase after each iteration. Matrix relations leading to Eq. 8 for the extended and rational Krylov subspaces are recalled in Appendix A.

Assume that ${\cal V}_{m}$ has orthonormal columns. Following similar reduction methods in the dynamical system contexts, see, e.g., [4], the reduction process consists of first projecting and restricting the original data onto the approximation space as

[TABLE]

Then the following low order differential Riccati equation needs to be solved,

[TABLE]

for $t\in[0,t_{f}]$. This low-dimensional DRE admits a unique solution for $t_{f}<\infty$ , see e.g., [28]. Restrictions on the data to allow for positive, stabilizing solutions are discussed in more detail in section 5.1. An approximation to the sought after solution is then written as

[TABLE]

We stress that $X_{m}(t)$ is never explicitly computed, but always referred to via the matrix $\mathcal{V}_{m}$ and the set of matrices $Y_{m}(t)$ at given time instances. In fact, the matrices $Y_{m}(t)$ may also be numerically low rank, so that at the end of the whole process a further reduction can be performed by truncating the eigendecomposition of $Y_{m}(t)$ for each $t$.

Remark \thetheorem.

The approach we have derived is solely based on the order reduction of the dynamical system (2). Nonetheless, and with some abuse of notation, the reduced DRE could have been formally obtained by means of a Galerkin condition on the differential equation. For $t\in[0,t_{f}]$ let

[TABLE]

be the residual matrix for $X_{m}(t)={\cal V}_{m}Y_{m}{\cal V}_{m}^{T}$. The matrix $Y_{m}(t)$ is thus determined by imposing that the residual satisfies the following Galerkin condition

[TABLE]

that is, $\mathcal{R}_{m}(t)\perp\mathcal{K}_{m}$ in a matrix sense, so that the residual is forced to belong to a smaller and smaller subspace as $\mathcal{K}_{m}$ grows. Substituting $X_{m}(t)=\mathcal{V}_{m}Y_{m}(t)\mathcal{V}_{m}^{T}$ into the residual matrix, the application of the Galerkin condition results in the projected system

[TABLE]

*which corresponds to (9). This is rigorous as long as ${\dot{X}}_{m}={\cal V}_{m}{\dot{Y}}_{m}{\cal V}_{m}^{T}$ holds. $\,\,\square$ *

It is crucial to realize that, as opposed to some available methods in the literature (such as, for instance, [41],[47] and the time-invariant algorithms in [30]), the approximation space is independent of the time stepping, that is a single space range$(\mathcal{V}_{m})$ is used for all time steps. This provides enormous memory savings whenever the approximate solution is required at different time instances in $[0,t_{f}]$. Theoretical motivation for keeping the approximation space independent of the time-stepping is contained in [5], where it is shown that the solution of the DRE lives in an invariant Krylov-subspace111 We also refer the reader to the recent manuscript [6], which appeared on-line briefly before the first round of our revision..

The class of numerical methods we used for solving the reduced DRE is described in the next section. In the rest of this paper, we specialize the generic derivation above to the extended and rational Krylov subspaces, which greatly outperformed polynomial spaces both in terms of CPU time and memory requirements. More information on these spaces and their properties are given in Appendix A; in particular, we discuss the generation of a real rational Krylov basis in the presence of non-real shifts.

3 BDF methods for the DRE

The numerical solution of the small-scale DRE is a well-studied topic, see, e.g., [34, 14, 35, 46]. Among the explored methods are matrix generalizations of the BDF methods [34, 15], which are computationally appealing only for small problems. Due to the reduction strength of rational Krylov subspaces, we expect the reduced DRE in Eq. 9 to be small enough to allow for efficient use of matrix-based BDF methods, which we are going to summarize next. For simplicity of exposition, in the rest of this section we omit the subscript in $Y_{m}$, and denote $Y^{(k+1)}=Y(t_{k+1})$. If we define

[TABLE]

the approximation of $Y^{(k+1)}$ is given by the implicit relation

[TABLE]

where $h=t_{k+1}-t_{k}$ is the stepsize and the respective $\alpha_{i}$’s and $\beta$ are the coefficients of the $b$-step BDF method for $b\leq 3$ and are given below.

[TABLE]

Substituting Eq. 12 into Eq. 13 results in the following nonlinear matrix equation

[TABLE]

which can be reformulated as the following continuous-time ARE

[TABLE]

The coefficient matrices are given by

[TABLE]

The Riccati equation Eq. 14 can be solved using “direct” methods; see, e.g., [9]. In our experiments we used the MATLAB solver care from the control systems toolbox. A brief sketch of the $b$-step BDF method is reported in Algorithm 1; other approaches are discussed, e.g., in [34, 30].

We conclude the section by depicting the typical convergence behavior of the BDF methods in our context. We consider an example from [35], where the $n\times n$ matrix $A$ stems from the spatial finite difference discretization of the following advection-diffusion equation

[TABLE]

on $\Omega=(0,1)^{2}$ with homogeneous Dirichlet boundary conditions. The choices of $B\in\mathbb{R}^{n\times 1}$ and $C\in\mathbb{R}^{1\times n}$ are given binomially as described in [35]. The initial condition is taken to be the zero matrix, that is $Z=\mathbf{0}_{n\times 1}$. We compare the obtained solution with a “reference” numerical solution $Y_{ref}(t)$ computed by an accurate but expensive method (the MATLAB function ode23s in our experiments), so that $n$ is kept small, $n=49$. The convergence behavior for $b=1,2,3$ and $\ell$ timesteps, with $\ell=10,100,1000$ is displayed in Fig. 1. The left plot shows the error $\|Y(t)-Y_{ref}(t)\|$ as a function of $t$, for different values of $\ell$. The right plot shows the evolution of the (1,1) component of the solution throughout the time span for the most accurate choice of BDF method, compared with that of the reference solution.

These plots illustrate that we cannot expect an overall high accuracy of the projection method as long as the reduced differential equation is not solved with sufficiently good accuracy. The importance of this is discussed in more detail in the following section.

4 Stopping criterion and the complete algorithm

To complete the reduction algorithm of section 2, we need to introduce a stopping criterion. We found that it is crucial to take into account the accuracy of the numerical method employed to solve the reduced DRE, as discussed in section 3.

To derive our stopping criterion we were inspired by those in [27, 23], however, we made some important modifications. In both cited references, the authors assume that the inner problem Eq. 9 is solved exactly, which is not true in general. We thus consider that the numerical method solves the reduced problem with residual matrix $R_{m}^{(I)}(t):=\dot{Y}_{m}(t)-\mathbf{F}(Y_{m}(t))$, so that the final DRE residual can be split into two components.

Proposition 4.1**.**

Let $X_{m}(t)=V_{m}Y_{m}(t)V_{m}^{T}$ be the Krylov-based approximate solution after $m$ iterations, where $Y_{m}(t)$ approximately solves the reduced problem Eq. 9. With the previous notation, the residual matrix $\mathcal{R}_{m}(t)=\dot{X}_{m}(t)-\mathbf{F}(X_{m}(t))$ satisfies

[TABLE]

*where $R_{m}^{(I)}(t)=\dot{Y}_{m}(t)-\mathbf{F}(Y_{m}(t))$ and $R_{m}^{(O)}(t)=\tau_{m}^{T}Y_{m}(t)$ with $\tau_{m}$ as in (8). *

Proof 4.2**.**

Substituting Eq. 10 into the residual $\mathcal{R}_{m}(t)$ we obtain

[TABLE]

Since $C^{T}$ belongs to range$({\cal V}_{m})$, we can write $C^{T}={\cal V}_{m}C_{m}^{T}$. Using Eq. 8, we get

[TABLE]

Since $\mathcal{V}_{m+1}=[\mathcal{V}_{m},\nu_{m+1}]$, we can write $\mathcal{R}_{m}(t)=\mathcal{V}_{m+1}\mathcal{J}_{m}(t)\mathcal{V}_{m+1}^{T}$, where

[TABLE]

Let $R_{m}^{(I)}(t)$ be the residual of the numerical ODE inner solver. Then

[TABLE]

Since the columns of $\mathcal{V}_{m+1}$ are orthonormal,

[TABLE]

*that is, $\|\mathcal{R}_{m}(t)\|_{F}^{2}=\|R_{m}^{(I)}(t)\|_{F}^{2}+2\|\tau_{m}^{T}Y_{m}(t)\|_{F}^{2}$, and the result follows. *

The expression for ${\cal J}_{m}(t)$ emphasizes that at each iteration $m$ the matrix $Y_{m}(t)$ is the exact solution of

[TABLE]

Hence, as long as $\|R_{m}^{(I)}(t)\|_{F}$ is not very small, the increase of $m$ aims at more and more accurately approximating a “nearby” differential problem to the truly projected one, with a term $R_{m}^{(I)}(t)$ that varies with $m$. Hence, $X_{m}(t)={\cal V}_{m}Y_{m}{\cal V}_{m}^{T}$ is an approximation not to $X(t)$, but to the solution of a differential problem with an additional term whose projection onto the space is $R_{m}^{(I)}(t)$.

Proposition 4.1 also implies that we cannot expect an overall small residual norm if either of the two partial residual norms $\|R_{m}^{(I)}(t)\|_{F}$, $\|R_{m}^{(O)}(t))\|_{F}$ is not small. In particular, we observe that the two residuals can be made small independently. Therefore we propose the following practical strategy:

1. (i)

Run the algorithm as presented, with a low-order cheap ODE inner solver (i.e., BDF(1, $\ell$) with $\ell$ relatively small) and use ${R}_{m}^{(O)}(t)$ in the stopping criterion; 2. (ii)

Once completed step (i) after $\widehat{m}$ iterations, use the matrices $T_{\widehat{m}},C_{\widehat{m}},B_{\widehat{m}}$ and $Z_{\widehat{m}}$ to refine the ODE inner solution by using a higher-order ODE solver for the projected system.

The final matrix $Y_{\widehat{m}}(t)$ obtained in step (ii) will provide a more accurate solution matrix than what would have been obtained at the end of step (i). We emphasize that any ODE method for small and medium scale DREs could be used at steps (i) and (ii). Our choice of BDF(1, $\ell)$ is due to its good trade-off between accuracy and computational effort; other approaches could be considered.

To complete the description of the stopping criterion, we recall that $R_{m}^{(O)}(t)$ depends on $t$, so that we need to estimate the integral over the whole time interval by means of a quadrature formula, that is

[TABLE]

where the interval $[0,t_{f}]$ has been divided into $\ell$ intervals with nodes $t_{j}$.

The overall algorithm222A Matlab implementation of both algorithms will be made available upon publication of this work. based on the rational Krylov subspace method is reported in Algorithm 2, while the algorithm based on the extended method is postponed to Appendix B. Several implementation issues of Algorithm 2 are also described in Appendix A, such as the use of a real basis in case of complex shifts $s_{j}$ in the basis construction.

5 Stability analysis and error bounds

In this section, we provide a few results on the spectral and convergence properties of the obtained approximate solution We first inspect some properties of the asymptotic matrix solution, which solves the algebraic Riccati equation. Then we propose a bound for the error matrix, in an appropriate functional norm.

5.1 Properties of the (steady state) algebraic Riccati equation

Properties associated with the algebraic Riccati equation – as asymptotic solution to the DRE – are well known in linear quadratic optimal control, see, e.g., [29, 13]. In particular, classical uniqueness and stabilization properties of the solution (see, e.g., [13, Lemma 12.7.2]), can directly be extended to the reduced DRE Eq. 9.

Corollary 5.1**.**

Let $({\cal T}_{m},B_{m},C_{m})$ be stabilizable and detectable system. Let $Y_{m}(t)$ be the solution of Eq. 9 at time $t$ and let $Y_{m}^{\infty}=\lim_{t\rightarrow\infty}Y_{m}(t)$. Then $Y_{m}^{\infty}$ is the unique symmetric nonnegative definite solution and the only stabilizing solution to the (reduced) algebraic Riccati equation

[TABLE]

*Moreover, if the pair $(C_{m},{\cal T}_{m})$ is observable, $Y_{m}^{\infty}$ is strictly positive definite. *

We notice that the stabilizability and detectability properties of $({\cal T}_{m},B_{m},C_{m})$ are not necessarily implied by those on $(A,B,C)$. Nevertheless, it is shown in [43] that if there exists a feedback matrix $K$, such that the linear dynamical system $\dot{x}=(A-BK)x$ is dissipative, then the pair $({\cal T}_{m},B_{m})$ is stabilizable. A similar result can be formulated for the detectability of $(C_{m},{\cal T}_{m})$, since by duality reasoning, $(C_{m},{\cal T}_{m})$ is detectable if $({\cal T}_{m}^{T},C_{m}^{T})$ is stabilizable. The question regarding the existence of such a feedback matrix, with respect to $A$ and $B$ ($A^{T}$ and $C^{T}$), is addressed in [22].

With these results, we can relate the asymptotic solution of the original and projected problems. Let $X_{m}(t)={\cal V}_{m}Y_{m}(t){\cal V}_{m}^{T}$ and $X_{m}^{a}={\cal V}_{m}Y_{m}^{\infty}{\cal V}_{m}^{T}$ respectively be approximate solutions to Eqs. 1 and 6 by a projection onto $\mbox{range}({\cal V}_{m})$. If there exist matrices $K$ and $L$ such that the systems $\dot{x}=(A-BK)x$ and $\dot{x}=(A^{T}-C^{T}L)x$ are dissipative, then

[TABLE]

that is, $X_{m}^{a}$ is the steady state solution of $X_{m}(t)$ when projected onto the same basis.

Under the hypotheses that $(A,B,C)$ is a stabilizable and detectable system, there exists a unique non-negative and stabilizing solution $X_{\infty}$ to Eq. 6 (see, e.g., [28, Theorem 5]). In [43] a bound was derived for the error $X_{\infty}-X_{m}^{a}$ in terms of the matrix residual norm. Here we complete the argument by stating that in exact arithmetic and if the whole space can be spanned, the obtained approximate solution equals $X_{\infty}$.

Proposition 5.2**.**

*Suppose $(A,B,C)$ is stabilizable and detectable. Assume it is possible to determine $m_{*}$ such that $\mbox{dim}(range({\cal V}_{m_{*}}))=n$, and let $X_{m_{*}}^{a}={\cal V}_{m_{*}}Y_{m_{*}}^{\infty}{\cal V}_{m_{*}}^{T}$ be the obtained approximate solution of Eq. 6 after $m_{*}$ iterations. Then, $X_{m_{*}}^{a}=X_{\infty}$. *

Proof 5.3**.**

Since ${\cal V}_{m_{*}}$ is square and orthogonal the projected ARE is given by

[TABLE]

*From $(A,B,C)$ stabilizable and detectable it follows that $({\cal V}_{m_{*}}^{T}A{\cal V}_{m_{*}},{\cal V}_{m_{*}}^{T}B,C{\cal V}_{m_{*}})$ is also stabilizable and detectable, so that $Y_{m_{*}}^{\infty}\geq 0$ and stabilizing. Multiplying by ${\cal V}_{m*}$ (by ${\cal V}_{m*}^{T}$) from the left (right), we obtain $0=A^{T}X_{m_{*}}^{a}+X_{m_{*}}^{a}A-X_{m_{*}}^{a}BB^{T}X_{m_{*}}^{a}+C^{T}C,$ that is, $X_{m_{*}}^{a}\geq 0$ is a solution to the original ARE. Since $X_{\infty}$ is the unique nonnegative definite solution, it must be $X_{m_{*}}^{a}=X_{\infty}$. *

5.2 Error bound for the differential Riccati equation

In this section we derive a bound for the maximum error obtained by the reduction process, in terms of the residual

[TABLE]

Note that $R_{m}(t)$ is the residual matrix with respect to the exact solution of the reduced differential problem, that is, it also includes the discretization error. A similar bound on the error has been derived for the nonsymmetric DRE in [3], which used matrix perturbation techniques from [26].

Proposition 5.4**.**

For $t\in[0,t_{f}]$ let ${\cal E}_{m}(t)=X(t)-X_{m}(t)$ and assume that ${\cal A}(t):=A-BB^{T}X(t)$ is stable for all $t\in[0,t_{f}]$. Denote

[TABLE]

where $\!\Phi_{{\cal A}}\!$ is the state-transition matrix satisfying $\frac{\partial{\Phi}_{{\cal A}}(t,s)}{\partial t}\!=\!{\cal A}(t)\Phi_{{\cal A}}(t,s),\Phi_{\cal A}(s,s)=I$.

If $4\nu^{2}\|B\|^{2}\|R_{m}\|_{\infty_{t}}<1$, then

[TABLE]

*where $\|L\|_{\infty_{t}}=\max_{t\in[0,t_{f}]}\|L(t)\|$ for any continuous matrix function $L(t)$. *

Proof 5.5**.**

By subtracting Eq. 20 from Eq. 1 and manipulating terms we obtain

[TABLE]

with ${\cal E}_{m}(0)=0$. Therefore, by the variation of constants formula (see, e.g., [28])

[TABLE]

Taking norms yields

[TABLE]

so that $\|{\cal E}_{m}(t)\|_{\infty_{t}}\leq\nu\left(\|R_{m}(t)\|_{\infty_{t}}+\|{\cal E}_{m}(t)\|_{\infty_{t}}^{2}\|B\|^{2}\right).$ Solving this quadratic inequality yields

[TABLE]

*The result follows from multiplying and dividing by $(1+\sqrt{1-4\nu^{2}\|B\|^{2}\|R_{m}\|_{\infty_{t}}})$ and noticing that at the denominator this quantity can be bounded from below by 1. *

We conclude with a remark on the intuitive fact that if the approximation space spans the whole space, the obtained solution by projection necessarily coincides with the sought after solution of the DRE.

Remark 5.6**.**

If it is possible to determine $m_{*}$ such that $\mbox{dim}({\cal V}_{m_{*}})=n$, then the approximate solution $X_{m_{*}}(t)$ coincides with $X(t)$ for all $t\geq 0$. Indeed, let us write $X_{m_{*}}(t)={\cal V}_{m_{*}}Y_{m_{*}}(t){\cal V}_{m_{*}}^{T}$, where ${\cal V}_{m_{*}}$ is square and orthogonal. The reduced DRE is given by

[TABLE]

with $Y_{m_{*}}=Y_{m_{*}}(t)$. Multiplying by ${\cal V}_{m*}$ (by ${\cal V}_{m*}^{T}$) from the left (right), we obtain

[TABLE]

*hence, $X_{m_{*}}(t)\geq 0$ is a solution of Eq. 1. Since $X(t)$ is the unique nonnegative definite solution of Eq. 1 for any $X_{0}\geq 0$ (see, e.g., [28]), then $X_{m_{*}}(t)=X(t)$ for $t\geq 0$. $\,\,\square$ *

6 Numerical experiments

In this section we report on our numerical experience with the developed techniques. We consider two artificial symmetric and nonsymmetric model problems, as well as three (of which two are nonsymmetric) standard benchmark problems. Information about the considered data is contained in Table 1. For the first two datasets displayed in Table 1, the matrix $A$ stems from the finite difference discretization with homogenous Dirichlet boundary conditions on the unit square and unit cube, respectively. The first matrix (sym2d) comes from the finite difference discretization of the two-dimensional Laplace operator in the unit square with homogeneous boundary conditions, while the second matrix (nsym3d) stems from the finite difference discretization of the three-dimensional differential operator

[TABLE]

in the unit cube, with homogeneous boundary conditions. For both datasets, the matrices $B,C$ and $Z$ are selected randomly with normally distributed entries. The realizations of the random matrices are fixed for both examples using the MATLAB command rng: for $B,C$ and $Z$ we use rng(7), rng(2) and rng(3), respectively. The following two datasets ( chip and flow) are taken from [1], and all coefficient matrices ($\widehat{A}$, $\widehat{B}$, $\widehat{C}$ and $\widehat{E}$) are contained in the datasets, which stem from the dynamical system

[TABLE]

Since $\widehat{E}$ is diagonal and nonsingular, it is incorporated as $A=\widehat{E}^{-\frac{1}{2}}\widehat{A}\widehat{E}^{-\frac{1}{2}}$, while $\widehat{B}$ and $\widehat{C}$ are updated accordingly to form $B$ and $C$.

The final considered dataset (rail) stems from a semi-discretized heat transfer problem for optimal cooling of steel profiles333Data available at http://modelreduction.org/index.php/Steel_Profile [8]. We consider the largest of the four available discretizations (file rail_79841_c60 containing $\widehat{A},\widehat{B},\widehat{C}$ and $\widehat{E}$) with $n=79841$. The symmetric and positive definite mass matrix $\widehat{E}$ has a sparsity pattern very similar to $\widehat{A}$. Both matrices are therefore reordered by the same approximate minimum degree (rksm-dre) or reverse Cuthill-McKee (eksm-dre) permutation to limit fill-in. The state-space transformation is done using the Cholesky factorization of $\widehat{E}$. More precisely, let $\widehat{E}=\widehat{E}_{L}\widehat{E}_{L}^{T}$ with $\widehat{E}_{L}$ lower triangular, and consider the transformed state $x=\widehat{E}_{L}^{T}\widehat{x}$. Then

[TABLE]

with $A=\widehat{E}_{L}^{-1}\widehat{A}\widehat{E}_{L}^{-T}$, $B=\widehat{E}_{L}^{-1}\widehat{B}$ and $C=\widehat{C}\widehat{E}_{L}^{-T}$. These matrices are never explicitly formed, rather they are commonly applied implicitly by solves with the factor $\widehat{E}_{L}$ at each iteration; see, e.g., [18, 42].

The initial low-rank factors are selected as the zero vector for flow, $Z=\sin{g}$ for chip and $Z=\cos{g}$ for rail, where $g\in\mathbb{R}^{n\times 1}$ is a vector with entries in $[0,2\pi]$. Other sufficiently general choices were tried during our numerical investigation however results did not significantly differ from the ones we report.

Performance of the projection methods. We first investigate the convergence behavior of the outer solver. The quantity we monitor in our stopping criterion is the backward error in an integral norm given by

[TABLE]

with $\rho_{m}$ as in Eq. 17 and

[TABLE]

The integrals are approximated by a quadrature formula in a similar fashion to Eq. 17, and we note that $\xi_{m}$ can be cheaply computed by using the Arnoldi-type relation.

For all datasets, the stopping tolerance was chosen as $10^{-7}$. For the first four datasets, $t_{f}=1$ and bdf(1,10) is used as inner solver. For rail, $t_{f}=4500$ (see e.g., [8] for further details about the setting) and bdf(1,45) is used as inner solver. Figs. 2 to 6 display the convergence of the rational Krylov subspace method (Algorithm 2, rksm-dre) and of the extended Krylov subspace method (Algorithm 3, eksm-dre). The left plots report the history of the backward error as the approximation space dimension increases, while the right plots display the same history versus the total computational time (in seconds) as the iterations proceed. We notice that the cost of the refinement step is not taken into account in these first tests.

For the dataset sym2d, the large algebraic linear system in rksm-dre was iteratively solved by implementing a block conjugate gradient algorithm, with an inner tolerance of $10^{-10}$, preconditioned with an incomplete Cholesky factorization with drop tolerance $10^{-4}$. For all other datasets, the MATLAB built-in backslash operator was used. For eksm-dre the coefficient matrix $A$ used to generate the Krylov space remains constant, hence a sparse reordered Cholesky (for sym2d and rail) or LU (for all other datasets) factorization was performed once and for all at the start of the algorithm. Therefore, only sparse triangular solves are required at each iteration. Clearly, the cost of the initial factorization depends on the size and density of the coefficient matrix. These two cost stages are particularly noticeable in the right plots of Fig. 3 and Fig. 4, where the eksm-dre curve starts towards the right of the plot, while the rest of the computation throughout the iterations is significantly faster.

In the implementation of rksm-dre it is possible to decide a priori whether to use only real or generically complex shifts. Our experiments showed that complex shifts were unnecessary for sym2d and nsym3d and, in fact, slowed down convergence when used. On the other hand, the use of general complex shifts proved to be crucial for the efficient convergence of rksm-dre for chip and flow. For the symmetric data in rail no complex shifts were used. We mention in passing that both algorithms are implemented so that the inner solves of Eq. 9 and the residual computations are performed at each iteration; for more demanding data we would advise a user to perform these computations only periodically to save on computational time.

Comparing performance, we observe that the two algorithms have alternating leadership in terms of computational time, but that rksm-dre almost consistently requires half the space dimension of eksm-dre. This is expected as the space dimension of eksm-dre increases with twice the number of columns per iteration, in comparison to rksm-dre. This observation is crucial at the refinement step, where it could be considerably more expensive to accurately integrate a DRE of dimension $2m(p+q)$ in comparison to a DRE with approximately half the dimension.

To have a clearer picture of how the various steps influence the performance of the methods, Table 2 depicts the overall computational time for the system solves, the orthogonalization steps and the integration of the reduced systems for each algorithm. For eksm-dre the CPU time required for the Cholesky and $LU$ factorizations are included in the solving time, but indicated in brackets as well. It is particularly interesting to notice the small percentage of time required by rksm-dre in comparison to eksm-dre for integrating the reduced system, confirming the comment made in the previous paragraph.

Comparisons with other BDF based methods. We compare the two projection methods rksm-dre and eksm-dre with low-rank methods that have been developed following different strategies. The package m.e.s.s. [41], for instance, can solve Lyapunov and Riccati equations, and perform model reduction of systems in state space and structured differential algebraic form, with time-variant and time-invariant data. For our purposes, the solvers in m.e.s.s. first discretize the time interval, and then solve the algebraic Riccati equation resulting from the ODE solver at each time step. Therefore, the approximation strategy employed at each time iteration to solve the algebraic problem is completely independent, and the obtained low-rank numerical solution needs to be stored separately. More precisely, if $\ell$ timesteps are performed, the procedure requires solving at least $\ell$ AREs of large dimensions, delivering the corresponding low-rank approximate solutions. Moreover, the rank of the constant term in the ARE increases with the time step, due to the way the ODE solver is structured, further increasing the complexity of the ARE numerical treatment. In our experiments with m.e.s.s. we only requested the approximate solution at the final stage. If the whole approximate solution matrix is requested at different time instances, the memory requirements will grow linearly with that. The overall strategy appears to be memory and computational time consuming, therefore we considered datasets of reduced size for our comparisons, as displayed in Table 3. The considered timespans were left unchanged.

Our experimental results are displayed in Tables 4 to 7; we remark that now also the refinement cost is taken into account in the projection methods. In all tables, the code bdf($b,\ell$) refers to the BDF method implemented in the refinement procedure of the reduction methods and in the time discretization procedure of m.e.s.s.

The tables show the storage requirements in terms of $n$-length vectors, the minimum and maximum approximate solution rank (with a truncation tolerance $10^{-8}$ for the projection methods) within the set of solutions, the CPU time break out of projection and refinement phases for the two projected methods, and finally the total CPU time. The stopping tolerance for all algebraic methods – that is the two projection methods and the Newton–Kleinmann-type method used in m.e.s.s. to solve each ARE – is set to $10^{-7}$.

In the m.e.s.s. software the user can either select a stopping tolerance (to be used for all solvers within the Newton–Kleinmann strategy) or a maximum number of iterations. We have experimented with both cases, where the maximum number of iterations was detected (a-posteriori) as the maximum number of iterations required within m.e.s.s to reach the tolerance of $10^{-7}$. It was observed that, in the majority of cases, avoiding the residual computation may, in fact, slow down the computational procedure. This is due to the possibility of performing several unnecessary iterations at some timesteps after the desired accuracy has in fact been reached. We therefore only report the results of the more realistic, reliable case where a stopping tolerance is selected beforehand. Galerkin acceleration is used to boost the performance of Newton–Kleinmann.

All numbers in the tables illustrate the large computational costs of m.e.s.s., as expected by the strategy “first time-discretize, then solve”, whereas both projection methods require just a few seconds of CPU in most cases.

The storage requirements of both reduction methods is independent of the number of timesteps where the solution is required. This is due to the fact that only a few $n$-long basis vectors need to be generated and stored, while only the reduced problem solution $Y_{m}(t)$ changes at the timesteps $t$. The memory requirements of m.e.s.s. are measured as the dimensions of the low-rank factor returned by the Newton-Kleinmann procedure, before column compression, at the final timestep. The dimension decreases significantly with the column compression. In our experiments we only stored the approximate solution at the last time step, however memory will be correspondingly higher if the whole approximation matrix is required at more instances (memory will thus grow linearly with the number of time instances to be monitored).

Between the two projection methods, we observe that the extended space yields a significantly larger basis than the actual approximate solution rank it produces. This means that the approximate solution belongs to a much smaller space than the one constructed by eksm-dre. This is far less so with rksm-dre. The different behavior confirms what has been already observed for the two methods in the ARE case [45].

Comparisons with splitting methods. We next compare rksm-dre with the fourth order additive splitting method (split-add4($\ell$)) developed in [47]. The method is based on splitting the DRE into the linear and non-linear subproblems, for which respective closed form solutions exist and are explicitly approximated. The numerical solutions to the subproblems are then recombined to approximate the solution to the full problem, by means of an additive splitting scheme. The main computational effort is due to the repeated evaluation of matrix exponentials, which has been resolved by using a Krylov-based matrix exponential approximation. Similar to the issue discussed with m.e.s.s. in the previous section, the $\ell$ (factored) solution matrices are independently calculated at each timestep, leading to significant memory requirements.

To ensure that we are comparing methods with similar approximation accuracies, we generate reference solutions $X_{ref}(t_{j})$ for the selected time instances $t_{j}$. This is done by using rksm-dre with a stopping tolerance of $10^{-10}$, plus a refinement process with bdf($4,10^{4}$) from [41]. To allow for such accurate approximations, we consider slightly smaller problem dimensions for the first two datasets, and we set $p=s=1$ and $X_{0}=0$.

The input parameters are tailored so that the approximate solutions from different methods have relatable accuracies. In particular, rksm-dre is solved with an outer stopping tolerance of $10^{-6}$ and with bdf(3,1000) in the refinement process. The number of timesteps utilized in split-add4 is selected as $\ell=500$. The expected approximation errors relative to the reference solution, measured as

[TABLE]

are illustrated in Fig. 7 (dataset sym2d in the left plot, dataset nsym3d in the right plot).

The figures indicate that we compare methods having approximation errors of similar order. The performance results are contained in Table 8 for two different discretizations of sym2d and nsym3d.

All numbers indicate the competitiveness of rksm-dre in terms of storage and computational time. The memory requirements for split-add4 is measured as the dimension of the solution factor at the final timestep, before column compression. If the solution is required at more time instances, then these memory requirements will increase accordingly.

We also mention that we have experimented with the dynamic splitting methods introduced in [35], however the algorithms proposed by the authors444We thank Chiara Piazzola for providing us with her Matlab implementation of the method. in [35] appeared to be better suited for small to medium size problems. Discussion on the refinement step. In previous sections, we have stressed that the two approximation stages of the projection method are independent, and we have focused on determining an effective approximation space. Here we linger over the accuracy of the second stage, the refinement step. Exploiting the far smaller problem size of the reduced problem, it is possible to allow for a much more accurate integration phase than what was done during the iteration of the reduction step. This crucial fact is already illustrated in the time break down of Tables 4 to 6, where especially for rksm-dre the refinement phase employs a fraction of the overall computational time, while still allowing for a rather accurate final solution.

We next explore in more detail these advantages with rksm-dre on sym2d, where the discretization was further refined to get a coefficient matrix of dimension $10^{6}$. The dimensions of the other corresponding matrices remain as presented in Table 1. We investigate the time taken by DRE solvers with different accuracies to emphasize the advantages and flexibility of the refinement procedure. Table 9 reports the timings for a refinement step performed by three different bdf methods and three splitting methods. The 8th order adaptive splitting method (split-adapt8) also comes from [47] and is performed with a tolerance of $10^{-7}$, i.e., the same as that of the reduction procedure. We emphasize that in the refinement phase we have utilized some of the most accurate integrators available, and nevertheless the high-dimensional ($n=10^{6}$) problem is approximated in less than 200 seconds for all integrators.

7 Conclusions and open problems

We have devised a rational Krylov subspace based order reduction method for solving the symmetric differential Riccati equation, providing a low-rank approximate solution matrix at selected time steps. A single projection space is generated for all time instances, and the space is expanded until the solution is sufficiently accurate. We stress that our approach is very general, and that it could be applied to subspaces other than Krylov-based ones, as long as the spaces are nested, so that they keep growing as the iterations proceed. This methodology could then be employed for more complex settings, such as parameter dependent problems, where the involved approximation space may require the inclusion of some parameter sampling.

Like in typical model order reduction strategies, in our methodology time stepping is only performed at the reduced level, so that the integration cost is drastically lower than what one would have by applying the time stepping on the original large dimensional problem. We have derived a new stopping criterion that takes into account the different approximation behavior of the algebraic and differential portions of the problem, together with a refinement procedure that is able to improve the final approximate solution by using a high-order integrator. These enhancement strategies have also been applied to the extended Krylov subspace approach. We have analyzed the asymptotic behavior of the reduced order solution, so as to ensure that the generated approximation behaves like the sought after time-dependent solution.

Although our numerical results are promising, there are still several open issues associated with the reduced order solution of the DRE. In particular, while stability and other matrix properties associated with the solutions $X(t)$ have been thoroughly studied [10, 19, 37, 16], the analysis of corresponding properties for the approximate solution $X_{m}(t)={\cal V}_{m}Y_{m}(t){\cal V}_{m}^{T}$ for $t\in[0,t_{f}]$ is still a largely open problem. In [27] some interesting monotonicity properties have been shown when the polynomial Krylov subspace is used together with particular ODE solvers; a complete analysis for $X_{m}(t)$ in a more general setting would be desirable.

Acknowledgements

We thank Khalide Jbilou for helpful explanations on [3], and Jens Saak for assistance in using the software M.E.S.S. v.1.0.1 [41]. We are also grateful to Chiara Piazzola for making her code from [35] available and to Tony Stillfjord for pointing us towards his codes from [47]. We thank Maximilian Behr for some helpful criticism on a previous version of this manuscript. The first author acknowledges the insightful comments and kind hospitality of the group of Peter Benner at the Max Planck Institute in Magdeburg (D), during a short visit in February 2019. Furthermore, we thank two anonymous referees for their careful reading and insightful comments.

This research is supported in part by the Alma Idea grant of the Alma Mater Studiorum Università di Bologna, and by INdAM-GNCS under the 2018 Project “Metodi numerici per equazioni lineari, non lineari e matriciali con applicazioni”. The second author is a member of INdAM-GNCS.

Appendix A Krylov subspace properties

In this Appendix we review some properties of extended and rational Krylov subspaces. As in Section 2 we denote $N=[C^{T},Z]$. Extended Krylov subspace. The extended Krylov subspace $\mathcal{EK}_{m}(A^{T},N)$ takes the form discussed in Section 2. The orthonormal basis $\mathcal{V}_{m}\in\mathbb{R}^{n\times 2m(p+q)}$ spanning the subspace is formed using the extended Arnoldi algorithm [17]. Let

[TABLE]

where $\mathcal{V}_{m+1}=[\mathcal{V}_{m}\hskip 2.84544ptV_{m+1}]\in\mathbb{R}^{n\times 2(m+1)(p+q)}$ and $E_{2m}$ is the last $2(p+q)$ columns of $I_{2m(p+q)}$. The extended Arnoldi algorithm produces the Arnoldi-type relation

[TABLE]

Rational Krylov subspace. The rational Krylov subspace was originally proposed in the eigenvalue context in [39]. Its use in our context is motivated by [45] and later [43], where its effectiveness in the solution of the algebraic Riccati equation is amply discussed.

Assume that $A$ is Hurwitz. Given ${\boldsymbol{s}}=\{s_{1},s_{2},\dots\}$, with $s_{j}\in\mathbb{C}^{+}$, the rational Krylov subspace is given by ${\cal RK}_{m}(A,N,{\boldsymbol{s}})$ as defined in Section 2. The approximation effectiveness of this subspace depends on the choice of shifts ${\boldsymbol{s}}$, and this issue has been investigated in the literature; see, e.g., [36], [18]. The adaptive choice of shifts was tailored to the ARE in [32] by the inclusion of information of the term $BB^{T}$ during the shift selection; see also [43] for a more detailed discussion555The Matlab code of the rational Krylov subspace method for ARE is available at http://www.dm.unibo.it/$\,\tilde{\,}$simoncin/software. In our numerical experiments we used this last adaptive strategy, where the approximate solution at timestep $t_{f}$ is used.

The algorithm presented in [18] forms a complex basis, when the shifts are not all real. In short, when $s_{j}\in\mathbb{C}^{+}$, the original approach would be to use the shift $s_{j}$ to form the next block $V_{j}$ and to then let the following shift be given by $s_{j+1}=\overline{s}_{j}$, where $\overline{s}_{j}$ denotes the complex conjugate of $s_{j}$. This results in both $V_{j}$ and $V_{j+1}$ being complex. As a consequence, the reduced DRE has complex coefficient matrices, although the final resulting approximations $X_{m}(t)$ will be real. Standard ODE solvers do not handle complex arithmetic well, hence we implemented an all-real basis using the method introduced in [40], which works as follows. If the shift $s_{j}$ is complex then the block $W_{j}=(A-s_{j}I)^{-1}V_{j-1}$ is also complex, hence we split it into its real and complex parts, that is $W_{j}=W_{j}^{(r)}+W_{j}^{(c)}\imath$. The block $V_{j}$ is then formed by orthogonalizing $W_{j}^{(r)}$ with respect to all vectors in the already computed basis, after which $V_{j+1}$ is formed by orthogonalizing $W_{j}^{(c)}$ with respect to all previous vectors in in the computed basis, and in $V_{j}$. This determines the same space, since span$\{W_{j},\bar{W}_{j}\}=$span$\{V_{j},V_{j+1}\}$. The resulting real basis of the rational Krylov subspace is given by $\mathcal{V}_{m}=[V_{1},\dots,V_{m}]\in\mathbb{R}^{n\times m(p+q)}$. We also define the matrices $\mathcal{V}_{m+1}=[\mathcal{V}_{m},V_{m+1}]\in\mathbb{R}^{n\times(m+1)(p+q)}$ and the matrix

[TABLE]

where $r_{m+1,m}\in\mathbb{R}^{(p+q)\times(p+q)}$ and $E_{m}$ holds the last $(p+q)$ columns of $I_{m(p+q)}$. The matrix $\widetilde{\mathcal{H}}_{m}$ contains the orthogonalization coefficients obtained during the rational Arnoldi algorithm.

Let $\mathcal{T}_{m}^{T}=\mathcal{V}_{m}^{T}A^{T}\mathcal{V}_{m}\in\mathbb{R}^{m(p+q)\times m(p+q)}$. The rational Krylov basis satisfies the Arnoldi-type relation

[TABLE]

where $G_{m}^{T}=\mathbf{\gamma}r_{m+1,m}E_{m}^{T}\mathcal{H}_{m}^{-1}$ and the matrix $\widehat{V}_{m+1}$ is an orthonormal matrix such that

[TABLE]

is the $QR$ decomposition of the matrix on the right (see [18, 31]). The rational Krylov procedure requires as an extra input the (usually real) values $s_{0}^{(1)},s_{0}^{(2)}$, which form a rough approximation of a spectral region used to compute the next shift. The reader is referred to [18, 43] for implementation details. Further, for the computation of the term $G_{m}^{T}Y_{m}(t)$ contained in the residual computation of rksm-dre, we follow an accelerated computation technique presented in [18].

Appendix B Extended Krylov subspace based method

The extended Krylov (eksm-dre) subspace method for solving Eq. 1 is presented in Algorithm 3.