Large-scale discrete-time algebraic Riccati equations— Doubling algorithm and error analysis

Journal of Computational and Applied Mathematics 277 (2015) 115–126

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Large-scale discrete-time algebraic Riccati equations— Doubling algorithm and error analysis Eric King-wah Chu a , Peter Chang-Yi Weng b,∗ a

School of Mathematical Sciences, Building 28, Monash University, VIC 3800, Australia

b

Institute of Statistical Science, Academia Sinica, Taipei 115, Taiwan

article

info

Article history: Received 17 October 2013 Received in revised form 29 May 2014 MSC: 15A24 65F30 93C05

abstract We consider the numerical solution of large-scale discrete-time algebraic Riccati equations. The doubling algorithm is adapted, with the iterates for A not computed explicitly but recursively. The resulting algorithm is efficient, with computational complexity and memory requirement proportional to the size of the problem, and essentially converges quadratically. An error analysis, on the truncation of iterates, and some numerical results are presented. © 2014 Elsevier B.V. All rights reserved.

Keywords: Discrete-time algebraic Riccati equation Doubling algorithm Large-scale problem

1. Introduction The nineteenth and twentieth centuries have seen tremendous achievements in science and technology. One hallmark of these achievements is the ability to influence and control dynamical systems and machines in our daily lives, ranging from economies of countries to domestic electrical appliances. All these are results of advances and sophistication in science, mathematics and engineering in general, and dynamical systems and control theory in particular. Various control tools have been proved successful, such as PID control, H∞ control, pole assignment and linear quadratic Gaussian (LQG) optimal control [1–6]. Our paper concerns large-scale discrete-time algebraic Riccati equations associated with the LQG optimal control problem. Any control process starts from modelling the system under investigation. Based on some experience of the system, a model, be it discrete- or continuous-time, time-variant or invariant, linear or nonlinear, is determined in terms of a mathematical model. This is then analysed with various tools in mathematics, towards its understanding and ultimate control, according to some performance criteria and within associated physical constraints. This model and system identification process is usually and necessarily a trade off between the accurate and detailed description of the system on one hand, and the ease of analysis and the ability to derive applicable results on the other. Too simple a model will not be accurate in reflecting the actual behaviour of our system, and too much complexity will harm the chance of achieving its understanding and ultimate control. Many systems have been assumed to be linear and time-invariant [1,7,5], which can be considered to be local approximations to more general nonlinear and time-variant systems. The theory of linear time-invariant (LTI) systems has been vital in control system theory, and has contributed successfully to many real-life applications.

∗

Corresponding author. Tel.: +886 963625216. E-mail addresses: [email protected] (E.K.-w. Chu), [email protected], [email protected] (P.C.-Y. Weng).

http://dx.doi.org/10.1016/j.cam.2014.09.005 0377-0427/© 2014 Elsevier B.V. All rights reserved.

116

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

1.1. Optimal control of linear time-invariant systems For continuous- or discrete-time LTI systems, one well-developed and much used technique is the LQG optimal control [1,2,6]. An energy or cost functional is minimized, constrained by the linear differential equations modelling the system, so that an optimal control can be obtained. With dynamic programming, calculus of variation or the Pontryagin principle [1,8,2,9,6], we obtained the optimality condition in the form of the algebraic Riccati equations [10,11] (ARE; for optimal control within finite time horizon, the differential or difference Riccati equations are obtained.) The optimal control we seek can be realized in terms of the unique maximal stabilizing solutions of these AREs. For discrete-time systems, the associated discrete-time ARE also plays an important role in the Kalman filter [8], reputedly one of the brightest jewels in twentieth century mathematics and system theory. All these contributed to tremendous amount of research for AREs in the past fifty years, since Kalman’s pioneering work in 1960 [12,10]. The numerical solution of AREs has been an active area of research, due to the importance of the equations in optimal control and filtering. Many practitioners in control theory and applied mathematics considered the problem, contributed dozens of methods [13,14,7], if not more. Each has its own favourites. Classical methods made use of canonical forms, determinants and polynomial manipulation, and state-of-the-art methods solve AREs in a numerical stable manner; see [13,14] and the references therein for more details. One favourite approach reformulates the AREs as eigenvalue value problems [15] and has been implemented in software packages like MATLAB [16] (as the commands care and dare). The other favourite is the Newton–Kleinman method [17]. Some of these methods work well for AREs from systems of moderate dimensions. For large-scale AREs, from applications such as PDE boundary control, efficient numerical methods have only been proposed recently; please consult [18–27] and references therein. Lastly, the numerics for linear systems benefits tremendously from advances in matrix theory and numerical linear algebra, associated with masters like Gantmacher, Bellman and others [1,8,2,3,11,4,5,28,6]. In many applications, various matrix operators have many exploitable structures, leading to many specialized efficient algorithms; for example, see [29,7,30,31]. 1.2. Discrete-time algebraic Riccati equations Consider the discrete-time algebraic Riccati equation (DARE):

D (X ) ≡ −X + A⊤ X (I + GX )−1 A + H = 0,

(1)

with the large and sparse system matrix A, where G = BR−1 B⊤ and H = CT −1 C ⊤ have low ranks, R ∈ Rm×m and T ∈ Rl×l are positive definite, B ∈ Rn×m and C ∈ Rn×l are orthogonal, m, l ≪ n, and (A, B) and (A, C ) are respectively stabilizable and detectable. DAREs arise often in linear-quadratic optimal control problems [11]. The solution of DAREs and continuous-time algebraic Riccati equations (CAREs) has been an extremely active area of research; see, e.g., [11]. The usual solution methods such as the Schur vector method [15], symplectic SR methods [32,18], the matrix sign or disk function [33,34] or the doubling method [13,14] have not made full use of the sparsity and structure in A, G and H. Requiring O(n3 ) flops and workspace of size O(n2 ), they are obviously inappropriate for large-scale problems. For control problems for parabolic PDEs and the balancing based model order reduction of large linear systems, largescale CAREs, DAREs, Lyapunov and Stein equations have to be solved; see [35,2,33,18–24,36–40]. Importantly, without solving the corresponding algebraic Riccati equations, alternative solutions to the optimal control problem require the deflating subspace of the corresponding Hamiltonian matrices or symplectic pencils which are prohibitively expensive to compute. Benner and his collaborators have done much on large-scale algebraic Riccati equations; see [34,18,41,19,40] and the references therein. They built their methods on (inexact) Newton’s methods with inner iterations for the associated Lyapunov and Stein equations. (The initialization of Newton’s method and the choice of parameters for the ADI inner iterations remain difficult.) Consult also [42,21–24,36,37] on the related invariant and Krylov subspace methods. We shall adapt the structure-preserving doubling algorithm (SDA) [13,14], making use of the sparsity in A and the fact that G, H and X possess (numerically) low ranks. Results on the numerical rank of the iterates of the SDA, neglected here, can be found in the technical report in [26, Section 2], which is an extended version of this paper. This work has also been applied and extended successively to other large-scale problems in [43,25,38,44,27]. The main results of the paper, for the error analysis in Section 3, have not appeared elsewhere. We shall use ∥ · ∥ to denote the 2-norm, unless otherwise stated. 1.3. Main contributions The paper considers the efficient solution of large-scale discrete-time algebraic Riccati equations using the doubling method. Our algorithm avoids the difficulties for the alternative Galerkin–Newton-ADI method [45], which involves the difficult initial stabilization for Newton’s method, or the difficult choice of parameters for the ADI acceleration. A detailed error analysis is provided.

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

117

1.4. Organization of paper We have presented some basics of control systems and AREs in the Introduction. The doubling algorithm and its modification for large-scale DAREs are described in Section 2, with its feasibility proved and computational complexity discussed. Sections 3 and 4 contain respectively an error analysis and some illustrative numerical experiments. The paper is concluded in Section 5. 2. Doubling algorithm From [46], we have the theorem on the solvability of the DARE (1): Theorem 2.1. Suppose that B and C are of full column rank, and R and T are positive. Assume that the pairs (A, B) and (A, C ) are stabilizable and detectable, respectively. Then the DARE (1) has a unique symmetric positive semi-definite solution X with ρ((I + GX )−1 A) < 1. Furthermore, if (A, C ) is completely reconstructible then the solution X of DARE is positive. We shall adapt the SDA [14] for large-scale DAREs. The main idea behind the SDA is to apply the doubling transformation to the symplectic matrix-pencil associated with the DARE. The doubling transformation preserves the deflating subspaces, while squaring the spectrum. After convergence, the stable deflating subspaces, and subsequently the maximal stabilizing solution X to the DARE, can be easily retrieved. The SDA [14] has the following form: Gk+1 = Gk + Ak (In + Gk Hk )−1 Gk A⊤ k ,

−1 Hk+1 = Hk + A⊤ Ak , k Hk (In + Gk Hk )

Ak+1 = Ak (In + Gk Hk )−1 Ak ,

(2)

assuming (In + Gk Hk )−1 exist, with A0 = A,

G0 = BR−1 B⊤ ,

H0 = CT −1 C ⊤ .

(3)

We have the following convergence theorem on the SDA [14]: Theorem 2.2. Assume that X , Y > 0 satisfy (1) and its dual equation Y = AY (I + HY )−1 A⊤ + G,

(4)

where G and H are given in (1), and let S = (I + GX )−1 A and S˘ = (I + HY )−1 A⊤ . Then the sequences {Ak }, {Gk } and {Hk } generated by the SDA (2) satisfy k

(a) Ak = (I + Gk X )S 2 ; k k (b) H 6 Hk 6 Hk+1 6 X , X − Hk 6 (S ⊤ )2 (X + XYX )S 2 ; and k k ⊤ 2 2 (c) G 6 Gk 6 Gk+1 6 Y , Y − Gk 6 (S˘ ) (Y + YXY )S˘ . 2.1. Large-scale SDA For large-scale DAREs, the solution X has low numerical rank. This is naturally exploited in the SDA, leading to an efficient algorithm which computes and stores the factors of approximations to X using smaller matrices. Let B0 = B, C0 = C , R0 = R−1 > 0 and T0 = T −1 > 0, where B, C , R and T are given in (1). Suppose that {Gk } and {Hk } are generated by SDA (2), Gk and Hk have the forms: Gk = Bk Rk B⊤ k ,

Hk = Ck Tk Ck⊤ ,

(k ≥ 0)

(5)

where Bk ∈ R , Ck ∈ R , R k ∈ R and Tk ∈ R . Applying the Sherman–Morrison–Woodbury formula (SMWF), we have n×mk

n×lk

mk ×mk

lk × lk

 −1 ⊤ (In + Gk Hk )−1 = In − Gk Ck Tk I + Ck⊤ Gk Ck Tk Ck  −1 ⊤ ⊤ = In − Bk I + Rk Bk Hk Bk Rk Bk Hk .

(6)

From the first glance, the iteration for A in the SDA in (2) appears doomed, with O(n ) operations for the products of full matrices. However, with the low rank and sparse-plus-low-rank (splr) forms in (1) and (6), we can organize Ak = A2k−1 − 3

(1)



(2)

Dk Sk Dk

⊤

(k ≥ 1). Eq. (6) then implies

Ak+1 = Ak (In + Gk Hk )−1 Ak  −1 ⊤  = Ak In − Gk Ck Tk I + Ck⊤ Gk Ck Tk Ck A k

    −1 = Ak In − Bk I + Rk B⊤ Rk B⊤ k Hk Bk k Hk Ak ⊤  = A2k − D(k1+)1 Sk+1 D(k2+)1 ,

(7)

118

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

with (1)

(2) Dk+1 = A⊤ k Ck ,

(1)

(2) Dk+1 = A⊤ k Hk Bk ,

Dk+1 = Ak Gk Ck , Dk+1 = Ak Bk ,

Sk+1 = Tk I + Ck⊤ Gk Ck Tk



Sk+1 = I + Rk B⊤ k Hk Bk



−1

−1

;

or

(8a)

Rk ;

(8b)

still involving O(n3 ) operations for a dense A. Note that the operation counts can be reduced to O(n) with the assumption that the maximum number of nonzero components in any row or column of A (denoted by nA ) is much less than and independent of n (see the flop counts in Table 1 in Section 2.4). One possibility is not to form Ak explicitly. Note that we have to store all the Bi , Ci , Ri and Ti for i = 0, 1, . . . , k − 1 to facilitate the multiplication of matrices of low ranks by Ak or A⊤ k . The induction proof of Ak in (7)–(8) is completed by the trivial k = 0 case.    Applying (6) to the SDA, we can group and organize the factors for the iterates  Gk+1 ≡  Bk+1 Rk+1 B⊤ k+1 and Hk+1 ≡ Ck+1 Tk+1

 Ck⊤+1 in terms of those in the previous iteration, with   Bk+1 = [Bk , Ak Bk ], Ck+1 = [Ck , A⊤ k Ck ];    −1 ⊤ ⊤  Rk+1 ≡ Rk ⊕ Rk − Rk B⊤ Ck Bk Rk k Ck Tk I + Ck Gk Ck Tk     −1 ≡ Rk ⊕ Rk − I + Rk B⊤ Rk B⊤ k Hk Bk k Hk Bk Rk

(9) (10a) (10b)

≡ Rk ⊕ R˘ k ,

(10c)

and

   −1   Tk+1 ≡ Tk ⊕ Tk − Tk Ck⊤ Gk Ck Tk I + Ck⊤ Gk Ck Tk     −1 ≡ Tk ⊕ Tk − Tk Ck⊤ Bk I + Rk B⊤ Rk B⊤ k Hk Bk k Ck Tk

(11a) (11b)

≡ Tk ⊕ T˘k . Note that  Rk ,  Tk and Sk in (8) are symmetric, and the spans of  Bk and  Ck are Krylov subspaces [26, Section 3.1].

(11c)

2.2. Truncation of Bk and Ck Now we shall consider an important aspect of the SDA for large-scale DAREs (SDA_ls) — the growth of Bk and Ck . Obviously, as the SDA converges, components of low ranks, which are increasing smaller but with higher ranks, are added to Bk and Ck . Apparent from (9), the growth in the sizes and ranks of these iterates is potentially exponential. To reduce the dimensions (1) (2) of Bk and Ck , and control the dimensions of Dk and Dk , we truncate their columns by orthogonalization. Assume without loss of generality that Bk and Ck in (5) are orthogonal. Consider the QR decompositions with column pivoting: I

Γ1b,2

Γ1b,3

 [Bk , Ak Bk ](I ⊕ Πb ) = [Bk , Φk , Φτ ,k ] 0

Γ2b,2

Γ2b,3  ,



0

0

I

Γ1c,2 Γ2c,2



[Ck , Ak Ck ](I ⊕ Πc ) = [Ck , Ψk , Ψτ ,k ] 0 ⊤

0



Γτb

 Γ1c,3 c Γ 2 ,3  ,

0

(12)

Γτ

c

where Πb , Πc are permutation matrices, ∥Γτ ∥ 6 ξb,k and ∥Γτc ∥ 6 ξc ,k . Here, ξb,k and ξc ,k are some tolerances controlling the truncation process. We denote b

Γu =



Γ1u,2 Γ2u,2

Γ1u,3 Γ2u,3



Πu⊤ ,

  Γu,τ = 0, Γτu Πu⊤ ,

u = Γ



Γu Γu,τ



,

 Mu =

I 0

Γu



,

(13)

where u = b or c. Let Bk+1 = [Bk , Φk ], τ

Rk+1

Ck+1 = [Ck , Ψk ],

= Mb Rk+1 Mb ⊤ ,

Tkτ+1 = Mc Tk+1 Mc ⊤ .

(14)

From (9)–(14), we then have τ ⊤  Gk+1 =  Bk+1 Rk+1 B⊤ k+1 = Bk+1 Rk+1 Bk+1 + △Gk+1 ,

 Hk+1 =  Ck+1 Tk+1 Ck⊤+1 = Ck+1 Tkτ+1 Ck⊤+1 + △Hk+1 ,

(15)

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

119

where R˘ k and T˘k are defined respectively in (10c) and (11c), and

 △Gk+1 = [Bk+1 , Φτ ,k ]

0

Γb,τ R˘ k Γb⊤

 △Hk+1 = [Ck+1 , Ψτ ,k ]

0

Γc ,τ T˘k Γc⊤

  Γb R˘ k Γb⊤,τ B⊤ k+1 , Φτ⊤,k Γb,τ R˘ k Γb⊤,τ   Γc T˘k Γc⊤,τ Ck⊤+1 . Ψτ⊤,k Γc ,τ T˘k Γc⊤,τ

(16)

Let Gτk+1 = Bk+1 Rτk+1 B⊤ k +1 ,

Hkτ+1 = Ck+1 Tkτ+1 Ck⊤+1 .

(17)

We now have all the additional components required for the SDA_ls, with the truncation of Bk and Ck , and the recursive form of Ak . The superscripts of Tkτ+1 and Rτk+1 , indicating truncated quantities, will later be dropped in Algorithm 1 for simplicity. In other words, we replace the notations Tkτ+1 and Rτk+1 defined in (14) by Tk+1 and Rk+1 , respectively, in Algorithm 1. Theorem 2.2 shows that the sequences {Gk } and {Hk } generated by the SDA are monotonically increasing and bounded above. In the following theorem, we shall prove the monotonicity of the sequences {Gτk } and {Hkτ } generated by the SDA_ls, where τ ⊤ Gτk = Bk Rk B⊤ k and Hk = Ck Tk Ck . To this end, we need the following useful results. τ ⊤ ˘ ˘ Lemma 2.3. Let Gτk = Bk Rk B⊤ k and Hk = Ck Tk Ck with Rk and Tk being positive definite. Then Rk and Tk in (10c) and (11c),

respectively, are positive definite and ∥R˘ k ∥ 6 ∥Rk ∥ = ∥Gτk ∥ and ∥T˘k ∥ 6 ∥Tk ∥ = ∥Hkτ ∥. Proof. From (10)–(11), we have ⊤ −1 ⊤ ⊤ −1 R˘ k = Rk − Rk B⊤ Ck Bk Rk = Rk (I + B⊤ , k Ck Tk (I + Ck Gk Ck Tk ) k Ck Tk Ck Bk Rk )

T˘k = Tk − Tk Ck⊤ Gk Ck Tk (I + Ck⊤ Gk Ck Tk )−1 = Tk (I + Ck⊤ Gk Ck Tk )−1 . Since Rk and Tk are positive definite, we deduce that R˘ k and T˘k are positive definite, and R˘ k 6 Rk and T˘k 6 Tk . As Bk and Ck are orthogonal, we have ∥R˘ k ∥ 6 ∥Rk ∥ = ∥Gτk ∥ and ∥T˘k ∥ 6 ∥Tk ∥ = ∥Hkτ ∥.  Theorem 2.4. For given tolerances η, ξb,k and ξc ,k , suppose that the sequences {Gτk } and {Hkτ } are generated by the SDA_ls, and τ ⊤ Gτk = Bk Rk B⊤ k and Hk = Ck Tk Ck . Then we have

(a) Gτk 6 Gτk+1 and Hkτ 6 Hkτ+1 ; (b) rankη (Gτk ) 6 rankη (Gτk+1 ) and rankη (Hkτ ) 6 rankη (Hkτ+1 ); and

(c) ∥Gτk+1 − Gτk ∥ = ∥Γb R˘ k Γb⊤ ∥ and ∥Hkτ+1 − Hkτ ∥ = ∥Γc T˘k Γc⊤ ∥, where R˘ k , T˘k and Γu (u = b or c) are given in (10), (11) and (13), respectively.

Proof. Since R0 and T0 are positive definite, from (10), (11), (14) and Lemma 2.3, it is easily seen that Rk and Tk are positive definite for each k ∈ N by induction. From (10), (11), (13), (14) and (17), we have Gτk+1 − Gτk = Φk Γb R˘ k Γb⊤ Φk⊤ ,

Hkτ+1 − Hkτ = Ψk Γc T˘k Γc⊤ Ψk⊤ .

(18)

It follows from Lemma 2.3 that Gτk 6 Gτk+1 and Hkτ 6 Hkτ+1 . This proves (a), and (b) follows directly. Since Bk+1 and Ck+1 are orthogonal, we can obtain (c) from (18).  For particular choices of ξb,k and ξc ,k , we shall show below that the truncation tolerances τg and τh are also the upper bounds of the relative truncation errors in Gτk+1 and Hkτ+1 , respectively. Theorem 2.5. For given truncation tolerances τg and τh , suppose that ξb,k = τg /(2∥Ak Bk ∥) and ξc ,k = τh /(2∥A⊤ k Ck ∥). Then

∥ Gk+1 − Gτk+1 ∥ 6 τg ∥Gτk+1 ∥,

∥ Hk+1 − Hkτ+1 ∥ 6 τh ∥Hkτ+1 ∥,

(19)

where Gτk+1 and Hkτ+1 are given in (17), and  Gk+1 and  Hk+1 in (15).

 Proof. With ∥Γb,τ ∥ = ∥Γτb ∥ 6 τg /(2∥Ak Bk ∥) and ∥Γc ,τ ∥ = ∥Γτc ∥ 6 τh /(2∥A⊤ k Ck ∥), (12) leads to ∥Ak Bk ∥ = ∥Γb ∥ > ∥Γb ∥  and ∥A⊤ C ∥ = ∥ Γ ∥ > ∥ Γ ∥ . From (15)–(16), Lemma 2.3 and Theorem 2.4(a), the truncation errors satisfy c c k k b ∥∥R˘ k ∥ 6 τg ∥Gτk ∥ 6 τg ∥Gτk+1 ∥, ∥ Gk+1 − Gτk+1 ∥ = ∥△Gk+1 ∥ 6 2∥Γb,τ ∥∥Γ c ∥∥T˘k ∥ 6 τh ∥Hkτ ∥ 6 τg ∥Hkτ+1 ∥.  ∥ Hk+1 − Hkτ+1 ∥ = ∥△Hk+1 ∥ 6 2∥Γc ,τ ∥∥Γ Remark 2.1. In the SDA_ls with given truncation tolerances τg and τh , we set ξb,k = τg /(2∥Ak Bk ∥) and ξc ,k = τh /(2∥A⊤ k Ck ∥). We can guarantee that τg and τh are the upper bounds of the relative truncation errors, ∥△Gk+1 ∥/∥Gτk+1 ∥ and ∥△Hk+1 ∥/ ∥Hkτ+1 ∥, respectively. As an additional safeguard, mk and lk (the widths of Bk and Ck respectively) can be restricted to be less than or equal to their respective maxima mmax and lmax in the SDA_ls. When mk+1 > mmax (or lk+1 > lmax ) for the given τg

120

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

(or τh ), we force mk+1 = mmax (or lk+1 = lmax ) and the corresponding error ∥△Gk+1 ∥ (or ∥△Hk+1 ∥) can still be computed by (16), and hence the resulting truncation tolerance (on the left-hand-side) and relative error (on the right-hand-side) are equal:

τg ,k+1 =



∥△Gk+1 ∥ ∥Rk+1 ∥

or τh,k+1 =

 ∥△Hk+1 ∥ . ∥Tk+1 ∥

(20)

2.3. Algorithm and residuals We now consider the difference of successive iterates, δ Gk+1 ≡ Gτk+1 − Gτk and δ Hk+1 ≡ Hkτ+1 − Hkτ . From Theorem 2.4(c), we have ∥δ Gk+1 ∥ = ∥Γb R˘ k Γb⊤ ∥ and ∥δ Hk+1 ∥ = ∥Γc T˘k Γc⊤ ∥, where R˘ k , T˘k and Γu (u = b or c) are given in (10), (11) and (13), respectively. Let Hkτ = Ck Tk Ck⊤ ,

τ −1 Tk+1 = Tk − Tk Ck⊤ B0 (I + R0 B⊤ R0 B⊤ 0 Hk B0 ) 0 Ck Tk ,

   Ck+1 ≡ Ck , A⊤ Ck .

Λk = A⊤ Ck Tk+1 Ck⊤ A, From the SDA_ls, we have

D (Hkτ ) = −Hkτ + Λk + H = −Ck Tk Ck⊤ + A⊤ Ck Tk+1 Ck⊤ A + CT0 C ⊤

= Ck+1

  T −T k + 0 0



0 0



⊕ Tk+1  Ck⊤+1 ≡  Ck+1 Tk+1 Ck⊤+1 .

Because Ck is orthogonal,  Ck+1 has the QR decomposition

 Ck+1 = [Ck , Θ ]



I 0

Υ



. ≡ [Ck , Θ ]Υ

For the residual rk ≡ ∥D (Hkτ )∥ and the relative residual  rk ≡ rk /(∥Hkτ ∥ + ∥Λk ∥ + ∥H ∥), we have the efficient formulae:

∥Hkτ ∥ = ∥Tk ∥,

∥H ∥ = ∥T0 ∥,

∥Λk ∥ = ∥Υ Tk+1 Υ ⊤ ∥,

  ⊤ ∥. ∥D (Hkτ )∥ = ∥Υ Tk+1 Υ

Algorithm 1 (SDA_ls) Input: Output:

A ∈ Rn×n , orthogonal B ∈ Rn×m and C ∈ Rn×l , R−1 = R−⊤ ∈ Rm×m , T −1 = T −⊤ ∈ Rl×l , positive tolerances τg , τh and ϵ , and mmax , lmax ; τg ,ϵ , τh,ϵ and Bϵ ∈ Rn×mϵ , Rϵ ∈ Rmϵ ×mϵ , Cϵ ∈ Rn×lϵ and Tϵ ∈ Rlϵ ×lϵ , with Cϵ Tϵ Cϵ⊤ (≈ X in (1)) and Bϵ Rϵ B⊤ ϵ (≈ Y in (4)); (1 )

(2)

Set A0 = A, D0 = D0 = 0 ∈ Rn×l , S0 = Il , B0 = B, R0 = R−1 , C0 = C , T0 = T −1 , δ0 = 2ϵ, τg ,ϵ = τg , τh,ϵ = τh and k = 0; Do until convergence: If δk < ϵ , Set Bϵ = Bk , Rϵ = Rk , Cϵ = Ck and Tϵ = Tk ; Exit; End If Compute  Bk+1 = [Bk , Ak Bk ],  Ck+1 = [Ck , A⊤ k Ck ]; 

 −1 ⊤ ⊤  Rk+1 = Rk ⊕ Rk − Rk B⊤ Ck Bk Rk , k Ck Tk I + Ck Gk Ck Tk     −1 ⊤ ⊤  Tk+1 = Tk ⊕ Tk − Tk Ck Bk I + Rk Bk Hk Bk Rk B⊤ k Ck Tk ;  ⊤ (2) (1 ) with Ak+1 = A2k − Dk+1 Sk+1 Dk+1 ,  −1 (1) (2) ⊤ ; Dk+1 = Ak Gk Ck , Dk+1 = Ak Ck , Sk+1 = Tk Il + Ck⊤ Gk Ck Tk Compute Bk+1 , Ck+1 by truncating  Bk+1 ,  Ck+1 with tolerances τg , τh ; modify Rk+1 ← Rτk+1 and Tk+1 ← Tkτ+1 ; If mk+1 > mmax (or lk+1 > lmax ), Set mk+1 ← mmax (or lk+1 ← lmax ); Compute τg ,k+1 (or τh,k+1 ) as in (20), and τg ,ϵ ← max{τg ,ϵ , τg ,k+1 } (or τh,ϵ ← max{τh,ϵ , τh,k+1 }); End If Compute k ← k + 1, δk = max{∥δ Hk ∥, ∥δ Gk ∥} and  rk .

End Do

Eqs. (7) (used recursively but not explicitly), (8a) (or (8b)), (9), (10a) (or (10b)), (11a) (or (11b)) and (17) (with the superscript (·)τ denoting truncation dropped for simpler notation), together with the corresponding initial condition (3), constitute the SDA_ls. We summarize our method in Algorithm 1, with the particular choice of (7), (9), (17), (10a), (11b) and (8a). We would like to emphasize again that care has to be exercised in Algorithm 1, where multiplications by Ak+1 and A⊤ k+1 are carried out recursively using (7) and (8a) or (8b). Otherwise, computation cannot be carried out in O(n) complexity. Similar care has to be taken in the computation of differences of iterates or residuals (used in Algorithm 1). 2.4. Operation and memory counts Recall that the maximum number of nonzero elements in any row or column in A is bounded by nA , for bounding the number of flops in the second column in Table 1. In the third column, the number of variables is recorded. Only O(n) operations

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

121

Table 1 Operation and memory counts for the kth Iteration in Algorithm 1 (SDA_ls). Computation

Flops

 Bk+1 ,  Ck+1  Rk+1 ,  Tk+1



Memory

2 nA (lk + mk ) + (6lk + 1)Nk n 2lk mk n

Nk+1 n O(l2k + m2k )

Dk+1 , Dk+1

2l2k n

–

Sk+1

O(l3k )

O(l2k )

Truncate Bk+1 , Ck+1

4(l2k + m2k )n

–

Modify Rk+1 , Tk+1

O(l3k + m3k )

–

r˜k+1

4l2k n

–

Total



(1 )

(2)



k+1

2k+1 nA (lk + mk ) + (6lk + 1)Nk + 2(5l2k + 2m2k + lk mk ) n



Nk+1 n

or memory requirement are included. Note that most of the work is done in the computation of  Bk+1 and  Ck+1 , for which Ak Bk k and A⊤ C have to be calculated recursively, as A is not available explicitly. We shall use the notation N k k ≡ k k j=1 (lj + mj ), and the operation count for the QR decomposition of an n × r matrix is 4nr 2 flops [29, p. 250]. With lk and mk controlled by the truncation in Section 2.2, the operation count will be dominated by the calculation of  Bk+1 and  Ck+1 . In our numerical examples in Section 4, the flop count at the end of Algorithm 1 dominates, with the work in one iteration approximately doubled that of the previous one. This corresponds to the 2k+1 factor in the total flop count. 3. Errors of SDA_ls Interesting links between the SDA and Krylov subspaces or Galerkin methods can be found in [26, Section 3] and will be ignored here. We shall perform an error analysis for the SDA_ls with truncation. From [14], the SDA is equivalent to a doubling process applied to the associated symplectic matrix pencils. The error analysis, conducted via these matrix pencils, produces some interesting and elegant results. We start from the equivalence of the DARE (1), for some stable matrix S = (I + GX )−1 A, to the symplectic eigenvalue problem



A −H

0 I

  I X



I = 0

G A⊤

 

I S. X

The kth iteration in the SDA is then equivalent to

  Mk

I X

  = Lk

k I S2 , X

 Mk ≡



Ak −H k



0 , I

Lk ≡

I 0



Gk . A⊤ k

(21)

The SDA is equivalent to the following doubling action, finding (Mk∗ , Lk∗ ) such that Mk∗ Lk = Lk∗ Mk and (Mk+1 , Lk+1 ) = (Mk∗ Mk , Lk∗ Lk ). The effect of this doubling action is the ‘‘doubling’’ (or squaring) of S in       k I I I Mk+1 = Mk∗ Mk = Mk∗ Lk S2 X X X       k k+1 k+1 I I I S 2 = Lk∗ Lk S2 = Lk+1 S2 . = Lk∗ Mk X

X

X

k

Note, interestingly and importantly, that the stable deflating subspace (corresponding to eigenvalues in S 2 ) stays the same, spanned by the columns of [I , X ⊤ ]⊤ . Let Ak , Gk and Hk be perturbed to Aτk , Gτk and Hkτ , respectively, by errors from the truncation in Section 2.2, round-off or whatever sources. Let  τ    I Gτk k ≡ Ak τ 0 , k ≡ M L . −Hk I 0 Aτk ⊤ Suppose that ϵk1 and ϵk2 are the errors such that

      I I ϵ 2k   Mk = Lk S + k1 . X X ϵk2

(22)

k∗ L k = L k∗ M k with Perform doubling using the perturbed quantities, i.e., M  k∗ = M

Aτk (I + Gτk Hkτ )−1 −(Aτk )⊤ Hkτ (I + Gτk Hkτ )−1



0 , I

 k∗ = L

I 0

Aτk (I + Gτk Hkτ )−1 Gτk (Aτk )⊤ (I + Hkτ Gτk )−1

and

 k+1 ≡ M

 Ak+1 − Hk+1



0 k∗ M k , =M I

 k+1 ≡ L

I 0

  Gk+1 k∗ L k . =L  A⊤ k +1



122

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

We then have

  I X

k+1 M

      k∗ ϵk1 k I S 2k + M k∗ L k I = M k∗ M =M ϵk2 X X     k k∗ ϵk1 k I S 2 + M k∗ M =L ϵk2 X       k∗ ϵk1 k∗ ϵk1 S 2k + M k I S 2k+1 + L k∗ L =L ϵk2 ϵk2 X     k+1 I S 2k+1 + ϵˆk+1,1 , =L ϵˆk+1,2 X

(23)

where

    k k ϵk1 S 2 + Aτk (I + Gτk Hkτ )−1 ϵk1 + Aτk (I + Gτk Hkτ )−1 Gτk ϵk2 S 2 ϵˆk+1,1 . = k ϵˆk+1,2 ϵk2 + (Aτ )⊤ (I + H τ Gτ )−1 ϵk2 S 2 − (Aτ )⊤ (I + H τ Gτ )−1 H τ ϵk1 k

k

k

k

k

k

k

Suppose that △Gk+1 and △Hk+1 in (16) are respectively truncated from  Gk+1 and  Hk+1 , i.e.,  Gk+1 = Gτk+1 +△Gk+1 and  Hk+1 = Hkτ+1 + △Hk+1 . From (22) and (23), we have       k+1 I = L k+1 I S 2k+1 + ϵk+1,1 , M X X ϵk+1,2 where

      2k+1 ϵk+1,1 ϵˆ = k+1,1 + △Gk+1 XS . ϵk+1,2 ϵˆk+1,2 △Hk+1

(24)

Let τg and τh be the small truncation tolerances for the relative truncation errors of  Gk and  Hk respectively, i.e., ∥△Gk ∥ 6 τg ∥Gτk ∥ and ∥△Hk ∥ 6 τh ∥Hkτ ∥ (c.f. (19)). Suppose that there is only truncation from the k0 th iteration and that the SDA_ls converges after kmax iterations (i.e., δkmax = max{∥δ Hkmax ∥, ∥δ Gkmax ∥} < ϵ , where ϵ is some accuracy tolerance). Let sk ≡ max{∥Aτk (I + Gτk Hkτ )−1 ∥, ∥Aτk ⊤ (I + Hkτ Gτk )−1 ∥}, h ≡ ∥Hkτmax ∥,

g ≡ ∥Gτkmax ∥,

(25)

ρ ≡ ρ(S ),

where ρ ∈ (0, 1) is the spectral radius of S. Note from Theorem 2.4 that ∥Gτk ∥ 6 g and ∥Hkτ ∥ 6 h for k = k0 , . . . , kmax . Since the sequence {Ak } generated by the SDA in (2) converges to zero, we expect {sk } to converge to a value near zero. Suppose ˘ S˘ being diagonal. Since that the stable matrices S and S˘ defined in Theorem 2.2 are diagonalizable, with XS−1 SXS and X ˘−1 SX S

ρ(S ) = ρ(S˘ ), we have ∥S ∥ 6 κs ρ and ∥S˘ ∥ 6 κs˘ ρ , where κs ≡ ∥XS ∥∥XS−1 ∥ > 1 and κs˘ ≡ ∥XS˘ ∥∥XS˘−1 ∥ > 1. Let κ ≡ max{κs , κs˘ } > 1,

(26)

2k

it is easy to see that ∥S ∥ 6 κρ

2k

2k

2k

and ∥S˘ ∥ 6 κρ . Since ρ < 1, there exists a smallest k∗ ∈ N such that κρ  k0 k∗ 2 2k0 +1 ) · · · (1 + κρ 2 ) if k0 6 k∗ , ξ∗ = (1 + κρ )(1 + κρ 1 if k∗ < k0 .

We have the following result. Lemma 3.1. Let ρ ∈ (0, 1) and κ be constants defined in (25) and (26), respectively. Then kmax −1



ξ∗

i

(1 + κρ 2 ) 6

,

k0

1 − ρ2

i=k0

where ξ∗ is given in (27). Proof. Let k∗ be the smallest positive integer such that κ∗ = κρ 2 kmax −1



i

(1 + κρ 2 ) 6

i=k0

∞ 

k∗

6 1. If k0 6 k∗ , then

i

(1 + κρ 2 )

i=k0 k +1

k

k∗

= (1 + κρ 2 0 )(1 + κρ 2 0 ) · · · (1 + κρ 2 )

∞  

i

1 + κ∗ ρ 2 −2

k∗



i=k∗ +1

6 ξ∗

∞  

1 + κ∗ ρ

2i

i=k∗

Similarly, if k∗ < k0 then

kmax −1 i=k0



6 ξ∗

∞  

1 + κ∗ ρ

i=k0 i

2i



6 ξ∗

∞ 

k

ρ j2 0 6

j =0 k

k

(1 + κρ 2 ) 6 1/(1 − ρ 2 0 ) = ξ∗ /(1 − ρ 2 0 ).



ξ∗ k0

1 − ρ2

.

2k∗

6 1. Denote

(27)

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

123

From (23) and (24), we obtain



  k ∥ϵk+1,1 ∥ sk + κρ 2 6 ∥ϵk+1,2 ∥ sk h

sk g κρ 2



k

1 + sk κρ 2

k

   2k+1 ∥ϵk1 ∥ + τg g ∥X ∥κρ ∥ϵk2 ∥ τh h

(28)

k0

with ∥ϵk0 ,1 ∥ ≈ τg g ∥X ∥κρ 2 and ∥ϵk0 ,2 ∥ ≈ τh h, for each k = k0 , . . . , kmax . In the following, we give upper bounds on the relative forward errors in the SDA_ls. Theorem 3.2. Let X and Y be the solutions of DARE (1) and its dual (4). Suppose that the SDA_ls computes kmax iterations and executes the truncation procedure of Gk and Hk for k = k0 , . . . , kmax with relative truncation tolerances τg and τh as in (19), respectively. Then, we have

 τ H

kmax

 − X

6

∥Hkτmax ∥  τ G

kmax

 − Y

kmax −1 k=k0

1 − ρ2

1 − ρ2



k

τh +

k0

c ξ∗ (kmax − k0 + 1)

6

∥Gτkmax ∥ where c ≡

c ξ∗ (kmax − k0 + 1)



k0

τg g ∥X ∥κρ 2 0

+

h k

τg +



τh h∥Y ∥κρ 2 0

kmax

h

∥Aτkmax Y ∥κρ 2

,

(29a)

,

(29b)

kmax

 +

g

⊤ ∥Aτkmax X ∥κρ 2

g

max{1, sk (h + 1), sk (g + 1)}, sk , g , h and ρ are defined in (25), and κ and ξ∗ in (26) and (27), respectively. k

Proof. For each k = k0 , . . . , kmax − 1, we denote Ek = Lk + ρ 2 κ Uk where





sk Lk = sk h



0 , 1

1 Uk = 0



sk g . sk

From (28), we have



     2k0 2k0 +1 ∥ϵkmax ,1 ∥ τ g ∥ X ∥κρ τ g ∥ X ∥κρ g g 6 Ek0 + Ek0 +1 + ··· ∥ϵkmax ,2 ∥ τh h τh h     2kmax −1 2kmax τ g ∥ X ∥κρ τ g ∥ X ∥κρ g g + Ekmax −1 + , τh h τh h

(30)

where Ek = Ekmax −1 Ekmax −2 · · · Ek , for k = k0 , . . . , kmax − 1. Let ℓk ≡ ∥Lk ∥1 = max{1, sk (1 + h)}, uk ≡ ∥Uk ∥1 = max{1, sk (1 + g )} and kmax −1

c=



kmax −1

max{ℓk , uk } =

k=k0



max{1, sk (h + 1), sk (g + 1)}.

k=k0

From the definition of Ek and Lemma 3.1, we obtain that for each k > k0 ,

     kmax −1  2k  i Ek τg g ||X ||κρ  6 τh h + τg g ∥X ∥κρ 2k (ℓi + ui κρ 2 )   τh h 1

i=k



6 c τh h + τg g ∥X ∥κρ 2

k0

 kmax −1

i

(1 + κρ 2 ) 6

i=k0

c ξ∗ 1−ρ

 2k0

 k τh h + τg g ∥X ∥κρ 2 0 .

It follows from (30) that

∥ϵkmax ,2 ∥ 6

c ξ∗ (kmax − k0 + 1)  1 − ρ2

 k τh h + τg g ∥X ∥κρ 2 0 .

k0

⊤ ⊤ Since ∥S 2 ∥ 6 κρ 2 , we have ∥Aτkmax XS 2 ∥ 6 ∥Aτkmax X ∥κρ 2 Similarly, the dual equation (4) can be rewritten as kmax



A⊤ −G

kmax

 

0 I

I Y

kmax

 =

I 0

H A

kmax

. From (22) and h = ∥Hkτmax ∥, we obtain (29a).

 

I ˘ S, Y

where S˘ = (I + HY )−1 A⊤ , and (31) implies (29b).

(31) 

124

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

Table 2 Example 4.1 (n = 20209, τg = 10−4 , τh = 10−9 , ϵ = 10−10 , mmax = lmax = 200). k

δk

∥Hkτ ∥

∥Gτk ∥

rk

 rk

mk

lk

δ tk

tk

1 2 3 4 5 6 7 8 9 10 11 12

1.079e−01 2.027e−01 3.880e−01 7.113e−01 1.197e+00 1.701e+00 1.749e+00 9.838e−01 1.879e−01 5.330e−03 4.153e−06 1.324e−11

0.2226 0.4197 0.7994 1.511 2.707 4.407 6.153 7.129 7.313 7.318 7.318 7.318

3.393e−5 3.393e−5 4.914e−5 6.505e−5 7.496e−5 7.862e−5 7.924e−5 7.927e−5 7.927e−5 7.927e−5 7.927e−5 7.927e−5

1.021e−01 9.915e−02 9.352e−02 8.322e−02 6.594e−02 4.151e−02 1.658e−02 2.700e−03 7.415e−05 5.763e−08 1.773e−12 1.729e−12

2.327e−01 1.335e−01 7.073e−02 3.463e−02 1.536e−02 5.793e−03 1.566e−03 2.031e−04 5.131e−06 3.936e−09 1.211e−13 1.181e−13

14 23 34 46 60 76 94 113 132 143 143 143

12 17 17 17 17 17 17 17 17 17 17 17

7.6e−2 1.7e−1 3.0e−1 5.8e−1 1.1e+0 2.7e+0 4.9e+0 1.1e+1 2.3e+1 4.9e+1 1.0e+2 2.0e+2

7.6e−2 2.4e−1 5.5e−1 1.1e+0 2.3e+0 5.0e+0 9.9e+0 2.0e+1 4.3e+1 9.3e+1 2.0e+2 4.0e+2

⊤ Remark 3.1. In fact, the last term in (29a) is the upper bound for ∥Aτkmax XS 2

kmax

⊤ ∥Aτkmax XS

2kmax

∥/h and

⊤ ∥ ≈ ∥Aτkmax Hkτmax (I + Gτkmax Hkτmax )−1 Aτkmax ∥ ≈ δkmax +1 .

Thus, δkmax +1 is an important quantity for controlling the relative forward error of X . In the SDA_ls, we use δk < ϵ with some accuracy tolerance ϵ as the stopping criterion. Remark 3.2. For the fixed relative truncation tolerances τg and τh , the upper bounds of relative forward errors are given on the right-hand-side of (29). In the SDA_ls, the effective truncation tolerance τg ,ϵ (or τh,ϵ ) may be larger than the original tolerance τg (or τh ), if we enforce mk 6 mmax and lk 6 lmax . The upper bounds of the relative forward errors can be obtained from (29) by replacing τg (or τh ) with τg ,ϵ (or τh,ϵ ). k0

For large enough k0 near convergence, we have ρ 2 ≈ 0. Essentially, the truncation tolerances τh and τg controlled the accuracy in Hkτmax and Gτkmax , respectively. Our numerical experiments illustrate that this is indeed the case, with the accuracy of Hkτmax usually slightly better than the asymptotic bound c ξ∗ (kmax − k0 + 1)τh in (29a), possibly due to some cancellation of errors, or the approximation effect of the corresponding Krylov subspaces. Observe the expected ill-effect of ρ ≈ 1 in (29). Finally, we would like to remind readers that τg and τh control the relative errors in Gτk and Hkτ from any source, not just the truncation process. Notice the interesting fact that the forward errors on the left-hand-sides of (29) have been obtained directly from the SDA via the associated symplectic pencil (Mk , Lk ) in (21). This contrasts with most numerical algorithms that forward errors are difficult to derive, leading to the celebrated approach of backward error analysis from Wilkinson. Of course, ill-conditioning still occurs in DAREs, associated with the eigenvalue problem of (Mk , Lk ) and the closeness from unity of the spectral radius ρ of S. 4. Numerical experiments We shall only present one selected example, with two values of n. The suite of problems involves continuous-time systems originated from the boundary control problem modelling the cooling of rail sections. The PDE model was semidiscretized using 2D finite elements to a continuous-time linear system with n variables, where m = 7, l = 6 and n = 20,209, 79,841. We transformed these continuous-time systems (discretized with ∆t = 0.01) to discrete-time models [18], yielding Examples 4.1 and 4.2. All the numerical experiments were conducted using MATLAB Version 2010a on an iMac with a 2.97 GHz Intel i7 processor, 8 Gigabyte RAM and machine accuracy eps = 2.22 × 10−16 . For the numerical results, we present δk ≡ {∥δ Hk ∥, ∥δ Gk ∥}. In Algorithm 1, the stopping criterion is δk < ϵ with some accuracy tolerance ϵ ; please also consult the convergence results in Theorem 3.2 and Remark 3.1. Note that δ ti is the execution k time for the ith iteration and the sub-total execution time tk = i=1 δ ti . In Examples 4.1 and 4.2, the iterations in the SDA_ls are reported for a corresponding set of (τg , τh ). In Tables 2 and 3, k, δk , ∥Hkτ ∥, ∥Gτk ∥, rk , rk , mk , lk , δ tk and tk are displayed. Example 4.1 ([19], With n = 20209). From Table 2, the (near-)machine accuracy of O(10−13 ) (in relative residual) is τ τ achieved within 11 iterations, in 200 s (execution time). In the 12th iteration, the relative modification of H11 is δ12 /∥H11 ∥≈ −12 −9 1.81 × 10 , which is smaller than the given truncation tolerance τh = τh,ϵ = 10 . Hence, the relative forward error of τ the computed solution H11 is about 10−9 by (29a). Notice, from columns 3 and 4 in Table 2, the monotonicity of Gk and Hk , as summarized in Theorems 2.2(b, c) and 2.4(a). Observe also the near-quadratic convergence of δk and rk , and the doubling of δ tk and tk , with respect to increasing k.

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

125

Table 3 Example 4.2 (n = 79841, τg = 10−4 , τh = 10−9 , ϵ = 10−10 , mmax = lmax = 200). k

δk

∥Hkτ ∥

∥Gτk ∥

rk

 rk

mk

lk

δ tk

tk

1 2 3 4 5 6 7 8 9 10

1.510e−01 2.770e−01 4.669e−01 6.660e−01 6.894e−01 3.924e−01 7.651e−02 2.254e−03 1.890e−06 6.725e−12

0.3132 0.5880 1.055 1.720 2.409 2.798 2.873 2.875 2.875 2.875

5.552e−5 7.633e−5 9.039e−5 9.590e−5 9.689e−5 9.693e−5 9.693e−5 9.693e−5 9.693e−5 9.693e−5

1.425e−01 1.270e−01 1.009e−01 6.381e−02 2.573e−02 4.269e−03 1.216e−04 1.017e−07 3.122e−13 1.382e−14

2.218e−01 1.218e−01 5.737e−02 2.236e−02 6.172e−03 8.169e−04 2.139e−05 1.765e−08 5.415e−14 2.397e−15

14 26 40 55 74 92 106 113 113 113

12 17 17 17 17 17 17 17 17 17

1.7e−1 4.7e−1 1.3e+0 2.8e+0 5.6e+0 1.2e+1 2.2e+1 4.4e+1 8.3e+1 1.6e+2

1.7e−1 6.4e−1 1.9e+0 4.7e+0 1.0e+1 2.2e+1 4.4e+1 8.8e+1 1.7e+2 3.3e+2

Example 4.2 ([19], With n = 79841). From Table 3, a (near-)machine accuracy of O(10−15 ) (in relative residual) is achieved within 10 iterations, in 330 s (execution time). In the 10th iteration, the relative modification of H9τ is δ10 /∥H9τ ∥ ≈ 2.34 × 10−12 , which is smaller than the given truncation tolerance τh = τh,ϵ = 10−9 . The relative forward error of the computed solution H9τ is about 10−9 by (29a). The ratio between t10 for different values of n is 3.548, against the expected value of 79841/20209 ≈ 3.951. Similar ratios in [26, Section 4] illustrate the O(n) computational complexity of the SDA_ls. More numerical results can be found in [26]. Numerical examples for CAREs, which are transformed to DAREs by Cayley transform before the SDA is applied, can be found in [25]. 5. Conclusions We have proposed a structure-preserving doubling algorithm for the large-scale discrete-time algebraic Riccati equation (1), the SDA_ls, with A being large and sparse, and B and C have low ranks. Similar continuous-time algebraic Riccati equations can be treated after an application of a Cayley transform. We apply the Sherman–Morrison–Woodbury formula when appropriate and do not form Ak (the iterate for A) explicitly. For well-behaved DAREs (or CAREs), with eigenvalues of the corresponding symplectic pencil (or Hamiltonian matrix) not on or near the unit circle (or the imaginary axis) and I + Gk Hk being invertible for all k, approximations of low ranks to the solutions X and Y can be obtained efficiently. The convergence of the SDA_ls is quadratic (ignoring the truncation of Bk and Ck ), as shown in [13,14]. The computational complexity and memory requirement are both O(n), provided that the growth of Bk and Ck is controlled. We have also shown that the truncation process in Section 2.2 controls the relative errors in the approximate solutions. Similar to the methods in [18], our technique can be applied when A is large and sparse, or is a product (inverse) of such matrices (a matrix). The feasibility of the SDA_ls depends on whether Av and A⊤ v (or A−1 v and A−⊤ v for CAREs) can be formed efficiently, for an arbitrary v . Contrasting other methods, ours does not require any difficult initialization nor any inner ADI or Smith iterations, for the Stein equations arisen from the inexact Newton’s iteration. As illustrated by our numerical results, the SDA_ls solves large-scale DAREs efficiently. The main results in Section 3 concern the effects of truncation on the SDA_ls for large-scale DAREs. It will be interesting to extend this to other applications of the SDA, with their individual structures, assumptions and requirements. Acknowledgements We thank Profs. Y.-C. Kuo, T. Li and W.-W. Lin for their participation in the earlier part of the research project. Part of the work was completed when the first author visited the Shanghai Key Laboratory of Contemporary Applied Mathematics at Fudan University and the National Taiwan Normal University. The second author was supported by Monash Graduate and International Postgraduate Research Scholarships, and a Monash Postgraduate Publication Award. References [1] [2] [3] [4] [5] [6] [7]

B.D.O. Anderson, J.B. Moore, Linear Optimal Control, Prentice Hall, N.J, 1971. M. Athan, P.L. Falb, Optimal Control: An Introduction to the Theory and Its Applications, McGraw-Hill, New York, 1965. W.L. Brogan, Modern Control Theory, third ed., Prentice-Hall, London, 1991. K. Ogata, Modern Control Engineering, fifth ed., Prentice-Hall, London, 2010. W.M. Wonham, Linear Multivariable Control, Springer-Verlag, Berlin, 1970&1979. K. Zhou, J.C. Doyle, Essentials of Robust Control, Prentice Hall, N.J., 1998. V.L. Mehrmann, The Autonomous Linear Quadratic Control Problem, in: Lecture Notes in Control and Information Sciences, vol. 163, Springer Verlag, Berlin, 1991. [8] B.D.O. Anderson, J.B. Moore, Optimal Filtering, Prentice Hall, N.J, 1979. [9] L.S. Pontryagin, V.G. Boltyanskii, R.V. Gamkrelidze, E.F. Mishchenko, The Mathematical Theory of Optimal Processes, John Wiley & Sons, New York, 1962. [10] R.E. Kalman, Contributions to the theory of optimal control, Bol. Soc. Mat. Mex. 5 (1960) 102–119.

126

E.K.-w. Chu, P.C.-Y. Weng / Journal of Computational and Applied Mathematics 277 (2015) 115–126

[11] P. Lancaster, L. Rodman, Algebraic Riccati Equations, Clarendon Press, Oxford, 1995. [12] R.E. Kalman, A new approach to linear filtering and prediction problem, J. Basic Eng. (ASME) 82D (1960) 35–45. [13] E.K.-W. Chu, H.-Y. Fan, W.-W. Lin, A structure-preserving doubling algorithm for continuous-time algebraic Riccati equations, Linear Algebra Appl. 396 (2005) 55–80. [14] E.K.-W. Chu, H.-Y. Fan, W.-W. Lin, C.-S. Wang, A structure-preserving doubling algorithm for periodic discrete-time algebraic Riccati equations, Internat. J. Control 77 (2004) 767–788. [15] C.-H. Guo, A.J. Laub, A Schur method for solving algebraic Riccati equations, IEEE Trans. Automat. Control (1979) 913–921. [16] Mathworks, MATLAB User’s Guide, 2010. [17] D. Kleinman, On an iterative technique for Riccati equation computations, IEEE Trans. Automat. Control 13 (1968) 114–115. [18] P. Benner, H. Fassbender, On the numerical solution of large-scale sparse discrete-time Riccati equations, Adv. Comput. Math. 35 (2011) 119–147. [19] P. Benner, J. Saak, A semi-discretized heat transfer model for optimal cooling of steel profiles, in: P. Benner, V. Mehrmann, D.C. Sorensen (Eds.), Dimension Reduction of Large-Scale Systems, in: Lecture Notes in Computational Science and Engineering, 45, Springer-Verlag, Berlin/Heidelberg, Germany, 2005, pp. 353–356. [20] M. Heyouni, K. Jbilou, An extended block Arnoldi algorithm for large-scale solutions of the continuous-time algebraic Riccati equations, Electron. Trans. Numer. Anal. 33 (2009) 53–62. [21] I.M. Jaimoukha, E.M. Kasenally, Krylov subspace methods for solving large Lyapunov equations, SIAM J. Numer. Anal. 31 (1994) 227–251. [22] K. Jbilou, Block Krylov subspace methods for large continuous-time algebraic Riccati equations, Numer. Algorithms 34 (2003) 339–353. [23] K. Jbilou, An Arnoldi based algorithm for large algebraic Riccati equations, Appl. Math. Lett. 19 (2006) 437–444. [24] k. jbilou, Low rank approximate solutions to large Sylvester matrix equations, Appl. Math. Comput. 177 (2006) 365–376. [25] T. Li, E.K.-W. Chu, W.-W. Lin, C.-Y. Weng, Solving large-scale continuous-time algebraic Riccati equations by doubling, J. Comput. Appl. Math. 237 (2013) 373–383. [26] T. Li, Y.-C. Kuo, E.K.-W. Chu, W.-W. Lin, Solving Large-scale Discrete-time Algebraic Riccati Equations by Doubling, Technical Report, NCTS, National Chiao Tung University, Hsinchu, Taiwan, 2012. [27] P.C.-Y. Weng, H.-Y. Fan, E.K.-W. Chu, Low-rank approximation to the solution of a nonsymmetric algebraic Riccati equation from transport theory, Appl. Math. Comput. 219 (2012) 729–740. [28] L. Zadeh, C. Desoer, Linear System Theory, The State Space Approach, McGraw Hill, New York, 1963. [29] G.H. Golub, C.F. Van Loan, Matrix Computations, third ed., Johns Hopkins University Press, Baltimore, MD, 1996. [30] J.S. Respondek, Numerical recipes for the high efficient inverse of the confluent Vandermonde matrices, Appl. Math. Comput. 218 (2011) 2044–2054. [31] J.S. Respondek, Recursive numerical recipes for the high efficient inverse of the confluent Vandermonde matrices, Appl. Math. Comput. 225 (2013) 718–730. [32] P. Benner, H. Fassbender, The symplectic eigenvalue problem, the butterfly form, the SR algorithm, and the Lanczos method, Linear Algebra Appl. 275–276 (1998) 19–47. [33] P. Benner, Solving large-scale control problems, IEEE Control Syst. Mag, 14 (2004) 44–59. [34] P. Benner, Editorial of special issue on large-scale matrix equations of special type, Numer. Linear Algebra Appl. 15 (2008) 747–754. [35] A. Antoulas, Approximation of Large-Scale Dynamical Systems, SIAM Publications, Philadelphia, PA, 2005. [36] K. Jbilou, ADI preconditioned Krylov methods for large Lyapunov matrix equations, Linear Algebra Appl. 432 (2010) 2473–2485. [37] K. Jbilou, A. Riquet, Projection methods for large Lyapunov matrix equations, Linear Algebra Appl. 415 (2006) 344–358. [38] T. Li, C.-Y. Weng, E.K.-W. Chu, W.-W. Lin, Large-scale Stein and Lyapunov equations, Smith method, and applications, Numer. Algorithms 63 (2013) 727–752. [39] J.-R. Li, J. White, Low-rank solution of Lyapunov equations, SIAM Rev. 46 (2004) 693–713. [40] J. Saak, H. Mena, P. Benner, Matrix Equation Sparse Solvers (MESS): A Matlab Toolbox for the Solution of Sparse Large-Scale Matrix Equations, Chemnitz University of Technology, Germany, 2010. [41] P. Benner, J.-R. Li, T. Penzl, Numerical solution of large Lyapunov equations, Riccati equations and linear-quadratic control problems, Numer. Linear Algebra Appl. 15 (2008) 755–777. [42] L. Amodei, J.-M. Buchot, An invariant subspace method for large-scale algebraic Riccati equation, Appl. Numer. Math. 60 (2010) 1067–1082. [43] T. Li, E.K.-W. Chu, Y.-C. Kuo, W.-W. Lin, Solving large-scale nonsymmetric algebraic Riccati equations by doubling, SIAM J. Matrix Anal. Appl. 34 (2013) 1129–1147. [44] C.-Y. Weng, E.K.-W. Chu, Y.-C. Kuo, W.-W. Lin, Solving large-scale nonlinear matrix equations by doubling, Linear Algebra Appl. 439 (2012) 914–932. [45] P. Benner, J. Saak, A Galerkin–Newton–ADI method for solving large-scale algebraic Riccati equations, DFG Priority Programme 1253 Optimization with Partial Differential Equations, Preprint SPP1253-090. [46] T. Pappas, A.J. Laub, N.R. Sandell, On the numerical solution of the discrete-time algebraic Riccati equation, IEEE Trans. Automat. Control 25 (1980) 631–641.