Applied Mathematics and Computation 242 (2014) 716–728
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Condition number and backward errors of nonsymmetric algebraic Riccati equation q Lan-dong Liu a,⇑, Shu-fang Xu b a b
Department of Mathematics, China University of Mining and Technology, Beijing 100083, China School of Mathematics Sciences, Peking University, Beijing 100871, China
a r t i c l e
i n f o
Keywords: Nonsymmetric algebraic Riccati equation Perturbation analysis Condition number Backward error M-matrix Minimal nonnegative solution
a b s t r a c t The normwise condition number and two kinds of backward errors are considered when nonsymmetric algebraic Riccati equation (NARE) has the minimal nonnegative solution. Based on the techniques for the symmetric case, we apply the condition number theory developed by Rice to define condition number for NARE. The explicit expression is derived in a uniform manner. Meanwhile, two kinds of backward errors are defined and evaluated by the explicit formulas. Numerical experiments are listed to illustrate and compare the practical performance. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction We consider nonsymmetric algebraic Riccati equation (NARE)
XCX XD AX þ B ¼ 0;
ð1:1Þ mn
where A; B; C; D are real matrices of sizes m m; m n; n m; n n, respectively, and solution X 2 R . NARE (1.1) has many applications in transport theory and Markov models [4,13]. The solution of practical interest is the minimal nonnegative solution. So we assume that
K¼
D
C
B
A
ð1:2Þ
is a nonsingular M-matrix to guarantee NARE (1.1) has the minimal nonnegative solution. The numerical algorithm of NARE (1.1) has attracted some authors’ attention. Theory and numerical methods of nonsymmetric algebraic Riccati equation are well developed [3,1,2,4–6,8,10,11,26]. In particular, structure-preserving doubling algorithm [11,9,14,16] is a very efficient algorithm. For these algorithms details, please refer to the references therein. Our interest here is to discuss the perturbation analysis of NARE (1.1). As is known perturbation analysis contains forward perturbation analysis and backward perturbation analysis. The purpose of forward perturbation is to ascertain the stability of an equation or a problem itself. The result of the forward perturbation analysis may be a perturbation bound, or a condition number, or a perturbation expansion. The back perturbation analysis is to test the stability of a computation or an algorithm,
q This project is supported by National Natural Science Foundation of China under Grant (11371364) and the Fundamental Research Funds for the Central Universities (2009QS09). ⇑ Corresponding author. E-mail address:
[email protected] (L.-d. Liu).
http://dx.doi.org/10.1016/j.amc.2014.06.036 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728
717
and ascertain the accuracy of an approximate solution. The result of backward perturbation analysis may be a backward error or a residual. For more details about perturbation theory, please refer to Higham [12] and Sun [19]. For algebraic Riccati equation some perturbation analysis results have been achieved, such as symmetric algebraic Riccati equation [26,29], discrete-time Riccati equation [20–22,15,23,24]. The study of perturbation analysis of NARE (1.1) was developed by Xu [26]. He presented the perturbation bound and the residual bound of NARE (1.1). Guo and Bai [10] discussed the perturbation bound and structured condition number about the minimal nonnegative solution of NARE (1.1). Recently some authors researched on the mixed and componentwise condition numbers of algebraic Riccati equation. For example, Zhou et al. [29] discussed symmetric algebraic Riccati equation, and Liu [17] discussed nonsymmetric algebraic Riccati equation. Xue et al. [28,27] considered the relative perturbation theory and its entrywise relatively accurate numerical solutions of M-matrix algebraic Riccati equations and M-matrix Sylvester equations. The purpose of this paper is to evaluate the condition number and backward error of a minimal nonnegative solution to NARE (1.1). Condition number of NARE (1.1), as the measure of the sensitivity for the minimal nonnegative solution to small changes in the coefficient matrices, plays a key role in the perturbation theory. We will apply the condition number theory developed by Rice [18] to define condition number of NARE (1.1), and derive an explicit expression of condition number in a uniform manner. Backward error is an important concept proposed by Wilkinson [25] in the 1960s. Now it becomes one of the basic tools to evaluate the quality of the calculated solution. This paper will discuss two kinds of backward errors of NARE (1.1) and derive their corresponding explicit expressions. This research is greatly influenced by Sun [23]. Throughout this paper we will use the following notations: k kF kk k k2 vec I
X
rðA; B; C; DÞ
Frobenius norm the operator norm induced by the Frobenius norm in an associated matrix space Euclidean vector norm vectorization operator, which stacks the columns of a matrix one under another Kronecker product Hadamard product identity matrix with the size determined by the context Rmm Rmn Rnm Rnn qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kAk2F þ kBk2F þ kCk2F þ kDk2F
This paper is organized as follows. In Section 2 some preliminary results are given. In Section 3 we define the condition number for the minimal nonnegative solution of NARE (1.1) and obtain a calculable explicit expression. In Section 4 two kinds of backward errors are defined, meanwhile the computable explicit expressions are given. In Section 5 the results are illustrated by some numerical examples. 2. Preliminaries In this section we give some preliminary results. Lemma 1 [7] and Lemma 2 [26] are introduced. Lemma 1. Let K be a nonsingular M-matrix, then NARE (1.1) has a minimal nonnegative solution X, such that
DC ¼ D CX
ð2:1Þ
AC ¼ A XC
ð2:2Þ
and
are all nonsingular M-matrices. Lemma 2. Let the map f : Rn ! Rm be resolved into
f ðxÞ ¼ uðxÞ þ gðxÞ; where x 2 Rn . Assume that (1) f ðxÞ; uðxÞ; gðxÞ are continuous, (2) uðaxÞ ¼ auðxÞ; a 2 Rþ ; 8x 2 Rn , (3) limkxk!0 kgðxÞk ¼ 0. kxk Then for an arbitrary given vector p ¼ ðe1 ; e2 ; . . . ; en ÞT , here
ei > 0 ði ¼ 1; 2; . . . ; nÞ, it holds that
718
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728
c ¼ lim sup
d!0 kpxk6d
kf ðxÞk kuðxÞk ¼ max ; x–0 d px
ð2:3Þ
where p x ¼ ðe1 x1 ; e2 x2 ; . . . ; en xn ÞT ; x ¼ ðx1 ; x2 ; . . . ; xn ÞT 2 Rn . Define the linear operators £ : Rmn ! Rmn as follows
£W ¼ AC W þ WDC ;
W 2 Rmn :
ð2:4Þ
From Lemma 1 we know that AC and DC are all nonsingular M-matrices. Consequently it follows that the operator £ is invertible. 3. Condition number In this section we discuss the condition number of NARE (1.1). Let the coefficient matrices A; B; C; D of NARE (1.1) be slightly perturbed by DA; DB; DC; DD, and let
e ¼ A þ DA; A
e ¼ B þ DB; B
e ¼ C þ DC; C
e ¼ D þ DD: D
Assume that the matrix
" e¼ K
e D e B
e C e A
# ð3:1Þ
is still a nonsingular M-matrix. So nonsymmetric algebraic Riccati equation with small perturbation
eX eX eC eX eD eA eþB e¼0 X
ð3:2Þ
e . Then still has a minimal nonnegative solution, denoted as X
e X: DX ¼ X Applying the condition number theory developed by Rice [18], we define the normwise condition number CðXÞ as follows. Definition 3.1. The condition number CðXÞ of X is defined as
CðXÞ ¼ lim sup e!0 De
kDXkF ; ne
ð3:3Þ
where
De ¼
DA DB DC DD 6e ; ðDA; DB; DC; DDÞ 2 X : . ; ; ; a b c d
here
DA DB DC DD ¼ . ; ; ; a b c d
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 kDAkF kDBkF kDCkF kDDkF þ þ þ a b c d
and n; a; b; c; d are all positive parameters. Taking
n¼a¼b¼c¼d¼1 in (3.3) the absolute condition number C abs ðXÞ is given, and taking
n ¼ kXkF ;
a ¼ kAkF ; b ¼ kBkF ; c ¼ kCkF ; d ¼ kDkF
in (3.3) the relative condition number C rel ðXÞ is given. Now we derive an explicit expression of condition number CðXÞ defined as in (3.3). Theorem 3.1. Assume that the coefficient matrices A; B; C; D of NARE (1.1) satisfy ðA; B; C; DÞ 2 X. Let X be a minimal nonnegative solution to NARE (1.1). And let CðXÞ be the condition number defined as in (3.3). Then CðXÞ has the explicit expression
CðXÞ ¼
1 k Cr k 2 ; n
ð3:4Þ
where
Cr ¼ ½ aU r ; bSr ; cT r ; dV r in which
ð3:5Þ
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728
U r ¼ L1 ðX T Im Þ;
Sr ¼ L1 ;
T r ¼ L1 ðX T XÞ;
V r ¼ L1 ðIn XÞ
719
ð3:6Þ
here the matrix L is shown as follows
L ¼ In AC þ DTC Im :
ð3:7Þ
Proof. Subtracting (1.1) from (3.2), we get
£ðDXÞ ¼ AC DX þ DXDC ¼ E þ uðDXÞ;
ð3:8Þ
E ¼ DB DAX X DD þ X DCX;
ð3:9Þ
where
uðDXÞ ¼ DXðC þ DCÞDX ðDA X DCÞDX DXðDD DCXÞ:
ð3:10Þ
Since AC and DC are the nonsingular M-matrices, the linear operator £ is invertible, it follows from (3.8) that
DX ¼ £1 ðDB DAX X DD þ X DCXÞ þ Oð.ðDA; DB; DC; DDÞÞ2 :
ð3:11Þ
We substitute (3.11) into (3.3) and use Lemma 2 to omit the higher-order, it holds that
CðXÞ ¼
1 k£1 ðDB DAX X DD þ X DCXÞk
max : n 0–ðDA;DB;DC;DDÞ2X . DA ; DB ; DC ; DD a
b
c
ð3:12Þ
d
Let
DA
a
DB ¼ N; b
¼ M;
DC
c
¼ P;
DD ¼ Q; d
ð3:13Þ
so from (3.12) it follows that
CðXÞ ¼
1 k£1 ðaMX þ bN þ cXPX dXQ Þk max : n 0–ðM;N;P;Q Þ2X .ðM; N; P; Q Þ
ð3:14Þ
Let T
2
h ¼ ½vecðMÞT ; vecðNÞT ; vecðPÞT ; vecðQ ÞT 2 RðmþnÞ : Notice that the linear operator £ has the matrix representation (3.7). Hence, associating with (3.5) and (3.6), from (3.14) we get
CðXÞ ¼
1 kCr hk2 kCr k2 max ¼ : n h–0 khk2 n
ð3:15Þ
Remark 3.1. From (3.4)–(3.7), we have the relative condition number
1 k ½ kAkF U r ; kBkF Sr ; kCkF T r ; kDkF V r k2 kXkF
C rel ðXÞ ¼
ð3:16Þ
and the absolute condition number
C abs ðXÞ ¼ k ½ U r ; Sr ; T r ; V r k2 :
ð3:17Þ
Remark 3.2. From (3.11)–(3.15), we get the following estimate about the relative perturbation with respect to solution X
e Xk kX F 6 C rel ðXÞ.rel þ Oð.2rel Þ; kXkF
ð3:18Þ
where
.rel ¼ .
DA DB DC DD : ; ; ; kAkF kBkF kCkF kDkF
ð3:19Þ
720
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728
Remark 3.3. In order to study the sensitivity of the coefficient matrices perturbation, we define the partial condition numbers C A ðXÞ; C B ðXÞ; C C ðXÞ and C D ðXÞ, respectively as
C A ðXÞ ¼ lim sup e!0
kDAkF 6ae DB¼DC¼DD¼0
C B ðXÞ ¼ lim sup e!0
kDBkF 6be DA¼DC¼DD¼0
C C ðXÞ ¼ lim sup e!0
kDCkF 6ce DA¼DB¼DD¼0
kDXkF ; ne
ð3:20Þ
kDXkF ; ne
ð3:21Þ
kDXkF ne
ð3:22Þ
kDXkF ; ne
ð3:23Þ
and
C D ðXÞ ¼ lim sup e!0
kDDkF 6de DA¼DB¼DC¼0
where DX is expressed by (3.11). Using the technique described in Theorem 3.1 we get the following expressions:
C A ðXÞ ¼
a n
b C B ðXÞ ¼ kSr k; n
kU r k;
c
C C ðXÞ ¼ kT r k; n
d C D ðXÞ ¼ kV r k n
ð3:24Þ
here the matrices U r ; Sr ; T r ; V r are the same as (3.6). Taking
n¼a¼b¼c¼d¼1 abs abs abs in (3.24) the partial absolute condition numbers are given, denote them as C abs A ðXÞ; C B ðXÞ; C C ðXÞ and C D ðXÞ, respectively. And taking
n ¼ kXkF ;
a ¼ kAkF ; b ¼ kBkF ; c ¼ kCkF ; d ¼ kDkF
rel rel rel in (3.24) the partial relative condition numbers are given, denote them as C rel A ðXÞ; C B ðXÞ; C C ðXÞ and C D ðXÞ, respectively.
4. Backward errors 4.1. The residual bound The result of backward perturbation analysis may be a backward error or a residual bound. The residual bound of NARE (1.1) has been discussed by Xu [26]. In Section 5, we will use the residual bound to demonstrate the validity of numerical examples. So the main results about the residual bound of NARE (1.1) are listed below. For the reduction details, please refer to Xu [26]. e be an approximate minimal nonnegative solution of NARE (1.1), the residual matrix R is Let X
eCX eX e D AX e þ B: R¼X
ð4:1Þ
Then we introduce the following symbols
bC ¼ X e C A; A
bC ¼ CX e D: D
ð4:2Þ
b be the linear operator defined as follows And let £
bcW þ W D b ¼A b c; £W
W 2 Rmn :
ð4:3Þ
Furthermore we have
^l ¼ k £ bCW þ W D b 1 k1 ¼ min k A bCk kWk¼1 W2Rmn
ð4:4Þ
and denote
e ¼ k b£ 1 Rk; c ¼ kCk:
ð4:5Þ
b c be the linear operator defined as follows Next let £
bc Z ¼ D b TZ þ ZD b c; £ c
Z 2 Rnn ;
ð4:6Þ
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728
721
we have
^lc ¼ k £ b 1 k1 ¼ min k D b TZ þ ZD b C k: c c
ð4:7Þ
kZk¼1 Z2Rnn
With these preparations, we have the following main result about the residual bound. Theorem 4.1 (Xu [26]). Assumptions and symbols are the same as before. If
e<
1 min f^l; ^lc g; 4c
ð4:8Þ
then the residual bound is
e Xk kX 2^l e qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 , bResB : ek ek kX ^l þ ^l2 4c^le k X
ð4:9Þ
Next we mainly discuss the normwise backward error of the minimal nonnegative solution of NARE (1.1). There are some e In the following subsection we propose two kinds different ways to define backward errors of NARE (1.1) with respect to X. of normwise backward errors. eÞ 4.2. The backward error gð X e does not satisfy NARE (1.1), but there exists an infinite number of In general, the approximate solution X e is a solution of the following equation ðDA; DB; DC; DDÞ 2 X such that X
e ðC þ DCÞ X eX e ðD þ DDÞ ðA þ DAÞ X e þ B þ DB ¼ 0: X
ð4:10Þ
Based on the back stability of an algorithm, we focus on the smallest perturbations (DA; DB; DC; DD) which satisfy Eq. (4.10), and we hope to know which one is sufficiently close to the original one in all of NARE (4.10) and what are the distances? Specifically, let all of (DA; DB; DC; DD) which satisfy Eq. (4.10) be in set S, that is
S ¼ fðDA; DB; DC; DDÞ 2 X : DA; DB; DC; DD satisfy NARE ð4:10Þg:
ð4:11Þ
e is given by the following definition. Now the normwise backward error of X e Þ of an approximate solution X e for NARE (4.10) is defined as Definition 4.1. The normwise backward error gð X
gð Xe Þ ¼
1
min
s ðDA;DB;DC;DDÞ2S
rðDA; DB; DC; DDÞ;
ð4:12Þ
where s is a positive parameter and set S is shown as (4.11). e Þ is called the normwise absolute backward error, denoted as g ð X e Þ. If s ¼ 1, then gð X abs e Þ is called the normwise relative backward error, denoted as g ð X e Þ. If s ¼ rðA; B; C; DÞ, then gð X rel e if gð XÞ e is small, then we conclude that the computing process for X e is For the computed minimal nonnegative solution X, e and optimal backward perturbabackward stable. The following result gives a calculable explicit expression between gð XÞ tions ðDAopt ; DBopt ; DC opt ; DDopt Þ. Theorem 4.2. Let
eCX eX e D AX e þB R¼X
ð4:13Þ
e , and denote that be the residual matrix of NARE (1.1) with respect to the approximate solution X
eX e T Þ1 ; W 1 ¼ ðIm þ X
eT X e Þ1 ; W 2 ¼ ðIn þ X
H ¼ W 1 RW 2 :
ð4:14Þ
Then
gð Xe Þ ¼
1
s
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e T k2 þ kHk2 þ k X e HX e T k2 þ k X e T Hk2 : kH X F F F F
ð4:15Þ
And the corresponding optimal backward perturbations DAopt ; DBopt ; DC opt ; DDopt are expressed by
eT; DAopt ¼ H X
DBopt ¼ H;
e HX eT; DC opt ¼ X
e T H: DDopt ¼ X
Proof. For any ðDA; DB; DC; DDÞ 2 S, from NARE (4.10) the following equation is obtained
ð4:16Þ
722
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728
e ðDCÞ X eX e ðDDÞ ðDAÞ X e þ DB ¼ R; X
ð4:17Þ
where R is shown as (4.13). In addition, (4.17) is equivalent to
2
vecðDAÞ
3
6 vecðDBÞ 7 7 6 M6 7 ¼ vecðRÞ; 4 vecðDCÞ 5
ð4:18Þ
vecðDDÞ where the matrix M is
h i e T Im ; Imn ; X eT X e ; In X e 2 RmnðmþnÞ2 : M ¼ X It holds that (4.18) has a unique minimal 2-norm solution
3 e T Im X 7 6 6 vecðDB Þ 7
1
1 7 6 opt 7 Imn 6 T T 1 eT X e eX eT 7 In þ X Im þ X vecðRÞ: 7 ¼ M ðMM Þ vecðRÞ ¼ 6 6 7 6 T e X e 5 4 vecðDC opt Þ 5 4 X e vecðDDopt Þ In X 2
vecðDAopt Þ
2
3
Thus from the notation of (4.14), it follows from above formula that
2
vecðDAopt Þ
3
2
eTÞ vecðH X
3
7 6 6 vecðDB Þ 7 6 vecðHÞ 7 opt 7 6 7: 7 ¼ 6 6 6 e X eTÞ 7 4 vecðDC opt Þ 5 5 4 vecð XH vecðDDopt Þ e T HÞ vecð X
ð4:19Þ
From (4.19) we obtain the explicit expressions of DAopt ; DBopt ; DC opt ; DDopt shown as (4.16). And by using Definition 4.1, formula (4.15) is obtained. h e Þ of an approximate solution X e as follows Remark 4.1. We also define the normwise backward error gð X
gð Xe Þ ¼
DA s ðDA;DB;DC;DDÞ2S DC
1
DB : DD F
min
e and taking s ¼ A Taking s ¼ 1 in (4.20) yields the absolute backward error gabs ð XÞ, C e Þ. error grel ð X
ð4:20Þ
B yields the relative backward D F
This definition is exactly the same as Definition 4.1. Therefore, the calculated formula is also equivalent to (4.15) obtained by Theorem 4.2. e 4.3. The backward error g ð XÞ e Þ and its computable expression. In this subsection we give another definition of the backward error g ð X e Þ of the approximate solution X e is defined as Definition 4.2. For NARE (4.10), the normwise backward error g ð X
g ð Xe Þ ¼
" DA a min ðDA;DB;DC;DDÞ2S DC c
# ; DD d DB b
ð4:21Þ
F
where a; b; c; d are all positive parameters and set S is shown as (4.11). e Þ is called the normwise absolute backward error, denoted as g ð XÞ. e If a ¼ b ¼ c ¼ d ¼ 1, then g ð X abs e e Þ. If a ¼ kAkF ; b ¼ kBkF ; c ¼ kCkF ; d ¼ kDkF , then g ð XÞ is called the normwise relative backward error, denoted as g rel ð X e The following result gives a computable formula of g ð X Þ. e for NARE (1.1). And denote Theorem 4.3. Let R be the residual matrix shown as (4.13) with respect to the approximate solution X the matrices Q ; G as follows
h i e T Im ; bImn ; c X eT X e ; dIn X e 2 RmnðmþnÞ2 ; Q ¼ a X G ¼ ðQQ T Þ
1
1 e T Þ2 Im þ b2 Imn þ c2 ð X e T Þ2 ð X e Þ2 þ d2 In ð X e Þ2 ¼ a2 ð X :
e Þ is expressed by Then the backward error g ð X
ð4:22Þ
ð4:23Þ
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728
e ¼ kQ T GvecðRÞk : g ð XÞ 2
723
ð4:24Þ
Proof. For any ðDA; DB; DC; DDÞ 2 S, from NARE (4.10) we have the following equation
e ðDCÞ X eX e ðDDÞ ðDAÞ X e þ DB ¼ R X
ð4:25Þ
or equivalently,
3 vec DaA 7 6 6 vec DB 7 7 6 b 7 Q6 7 ¼ vecðRÞ; 6 6 vec DcC 7 5 4 vec DdD 2
ð4:26Þ
where the matrix Q is shown as (4.22). cA; D cC ; d cB; D Let ð D DDÞ 2 S be the optimal backward perturbation, i.e.,
2 DbA a e g ð X Þ ¼ 6 4 DbC c
3 7 5 : c DD d DbB b
ð4:27Þ
F
Then from (4.26) we get
3 vec DaA 7 6 6 vec DB 7 7 6 b y T T 1 6 7 7 ¼ Q vecðRÞ ¼ Q ðQQ Þ vecðRÞ: 6 6 vec DcC 7 5 4 vec DdD 2
ð4:28Þ
Utilizing Kronecker product property ðA BÞðC DÞ ¼ AC BD, we have
e T Im Þð X e T Im Þ þ b2 Imn þ c2 ð X eT X e Þð X eT X e Þ þ d2 ðIn X e ÞðIn X eÞ QQ T ¼ a2 ð X e T Þ2 Im þ b2 Imn þ c2 ð X e T Þ2 ð X e Þ2 þ d2 In ð X e Þ2 : ¼ a2 ð X 1
We denote G ¼ ðQQ T Þ . Consequently, from (4.28) it follows that
3 vec DaA
7 6 6 vec DB 7 7 6 b T 6 7 7 ¼ Q GvecðRÞ; 6 6 vec DcC 7 5 4 vec DdD 2
ð4:29Þ
where G is shown as (4.23), So (4.29) implies formula (4.24). h
5. Numerical experiments In this section, some simple examples are given. We compute the absolute, relative condition number and two kinds of backward errors of NARE (1.1). They illustrate the validity of Theorems 3.1, 4.2, 4.3, respectively. All examples are carried out using MATLAB 7.10. Example 5.1. Consider NARE (1.1), in which
A¼
3 þ 10m 2
; 4 þ 10m 1
B¼C¼
1 1 1 1
;
and D ¼
4 þ 10m
2
1
3 þ 10m
;
where m is a positive integer. Note that the matrix K (shown as (1.2)) is a nonsingular M-matrix. Suppose that the perturbations in the coefficient matrices are
DA ¼ and
0:3 0:1 10j ; 0:1 0:2
DB ¼
0:1 0:2 10j 0:2 0:3
724
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728
DC ¼
0:3 0:1 0:1
0:2
10j ;
DD ¼
0:3
0:2
0:1
0:2
10j ;
where j is an integer. e ¼ A þ DA ; B e ¼ C þ DC ; D e ¼ B þ DB ; C e ¼ D þ DC be the coefficient matrices of the perturbed NARE (3.2). Assume Let A kAk kBk kCk kCk e are the minimal nonnegative solutions of NARE (1.1) and NARE (3.2), respectively. that X and X e in advance. In this paper we deal with this problem in the In general, we can not know the exact solutions X and X b of the unique following way. By using structure-preserving doubling algorithm (SDA [11]), we compute an approximate X b e of the unique minimal nonnegative solution X e to the minimal nonnegative solution X to NARE (1.1) and an approximate X b e have very high precision (see Table 3). Note that the approximate b and X perturbed NARE (3.2). The approximate solutions X b ke X b X kF ke X Xk of kXk F have at least 2 correct significant digits (A number’s significant digits are the digits of the mantissa values F kb X kF b eXkF e b can be substituted by k X X kF . For which do not count leading zeros). We deduce from these approximate solutions that k XkXk F kb X kF b b e k X Xk ke X X kF . the details, please see Sun [23]. So in all tables listed below, kXk F is replaced by the approximate computing values F kb X kF The numerical results in Table 1 and 2 display the relative perturbation bounds C rel ðXÞ .rel (.rel shown as (3.18)) and the condition numbers shown as Theorem 3.1 and Remark 3.3. we observe from Table 1 that the condition numbers of NARE (1.1) are increased while m is increased and j ¼ 12 is fixed. This indicates that Example 5.1 is an gradually ill-conditioned equation. The results listed in Table 1 in the relative sense display solution X is increasingly sensitive to the perturbations of matrices A; B; C; D. The fact is intuitively expected because matrix K is close to singular while m approaches infinity. Especially we notice that matrix K is a singular M-matrix while m ¼ þ1. That is
2 K¼
D
C
C
A
4
2 1 1
3
6 1 3 1 1 7 7 6 ¼6 7 4 1 1 3 1 5 1 1 2
4
is singular M-matrix. Table 2 exhibits the compared results in which m is fixed and j is increased. Both Tables 1 and 2 show that the relative perturbation bounds are fairly sharp. e Þ; g ð XÞ e are The backward errors and residual bound of Example 5.1 are listed in Table 3. The backward errors gabs ð X rel e g ð XÞ e are computed by formula (4.24). And the residual bound computed by formula (4.15). The backward errors g abs ð XÞ; rel e is the exact minimal bResB is computed by (4.9). The backward errors listed in Table 3 show that the computed solution X nonnegative solution to the slightly perturbed NARE (4.10). In other words, the computing process of the approximate solution is quite stable. Moreover, from the results listed in Table 3 we observe that in the case of m ¼ 1; 2; . . . ; 7 we have 10 e Xk 6 k X ek b kX ; ResB < 10 F F
ð5:1Þ
e approximate the exact solutions X up to about 10 significant digits in the sense which implies that the computed solutions X of the Frobenius norm. Example 5.2 [11]. We construct NARE (1.1) according to the rules given in [7]. We first generate and save a random 2n 2n nonzero matrix W by using randð2n; 2nÞ in MATLAB. Let the matrix K = diagðWeÞ W, with e ¼ ð1; 1; . . . ; 1ÞT 2 R2n . For a given nonnegative constant a, we define
A ¼ Kðn þ 1 : 2n; n þ 1 : 2nÞ þ aI;
B ¼ Kðn þ 1 : 2n; 1 : nÞ;
C ¼ Kð1 : n; n þ 1 : 2nÞ;
D ¼ Kð1 : n; 1 : nÞ þ aI:
Table 1 Relative perturbation bounds and condition numbers for Example 5.1 with different values of m (j ¼ 12).
.rel
ke X XkF kXkF
C rel A ðXÞ
C rel B ðXÞ
C rel C ðXÞ
C rel D ðXÞ
C rel ðXÞ
C rel ðXÞ
2
8:2184
3:5082
2:4902
8:1719
12:3611
1:8265 1012
1:3051 1013
4
82:0776
30:4695
29:4519
82:0196
123:5301
1:7901 1011
1:0657 1012
6
820:9522
300:2675
299:2500
820:8928
1:2359 10
10
1:7907 10
1:0464 1011
8
8:2097 103
2:9983 103
2:9972 103
8:2097 103
1:2360 104
1:7908 109
1:0414 1010
10
4
8:2097 10
4
2:9978 10
2:9977 10
4
4
5
8
12
8:2099 105
2:9979 105
2:9978 105
m
3
8:2097 10
1:2360 10
1:7908 10
1:0594 109
8:2099 105
1:2360 106
1:7908 107
1:0311 108
725
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728 Table 2 Relative perturbation bounds and condition numbers for Example 5.1 with different values of j (m ¼ 8). C rel ðXÞ .rel
ke X XkF kXkF
j
C abs ðXÞ
C rel ðXÞ
8
1:0461 106
1:7908 105
3:0566 103
1:2360 104
9
7
6
3
1:2359 103
1:0433 10
1:7908 10
3:0566 10
10
1:0429 108
1:7908 107
3:0566 103
1:2360 104
11
1:0414 10
9
8
3
1:2359 104
12
1:0414 1010
3:0566 103
1:2359 104
1:7908 10
3:0566 10
1:7908 109
Table 3 Backward errors and residual bounds for Example 5.1 with different values of m (j ¼ 12). e gabs ð XÞ
m 1
grel ð Xe Þ 13
2:08 10
e g abs ð XÞ
2:50 10
14
2:02 10
14
e g rel ð XÞ 13
8:84 1014
1:31 10
0.85
1:31 1013
1:67 10
6:54 10
1:58 1013
6:24 1014
3:47 1013
0.96
3:49 1013
13
14
12
1.00
1:07 1012
6:20 1014
3:29 1012
1.01
3:24 1012
14
12
1.02
1:05 1011
3:29 1011
1.02
3:32 1011
1:58 1013
1:92 1014
4
1:55 10
13
14
1:55 10
6:12 10
5
1:55 1013
1:88 1014
1:55 1013
6
1:54 10
13
14
13
7
1:54 1013
13
1:87 1014
0.57
13
ke X kF
14
1:67 10
1:87 10
8:88 10
14
8:69 10
2:08 10
3
1:87 10
kXe X kF
14
2
13
ek kX F
bResB
1:54 10
6:18 10
1:54 1013
6:14 1014
1:05 10
9:95 10
Note that K is a nonsingular M-matrix, A; B; C and D are n n matrices with A; D being nonsingular M-matrices and B; C being nonnegative matrices. The perturbations in the coefficient matrices are randomly generated ð0; 1Þ matrices by using MATLAB function rand, and
DA ¼ ðrandðsizeðAÞÞ=kAkÞ 1012 ;
DB ¼ ðrandðsizeðBÞÞ=kBkÞ 1012 ;
DC ¼ ðrandðsizeðCÞÞ=kCkÞ 1012 ;
DD ¼ ðrandðsizeðDÞÞ=kDkÞ 1012 :
ð5:2Þ
While a ¼ 0; a ¼ 5 and size of matrices n is taken from 5 to 50, we compare the relative perturbation bounds C rel ðXÞ .rel e e X kF X kF with kX . Fig. 1 and Fig. 2 display the compared results of the nature logarithm of C rel ðXÞ .rel and kX . The results show kXk kXk F
F
that the estimate C rel ðXÞ .rel gives an upper perturbation bound of X. And from Table 4 and Table 5 we see that the
-21
|| X||/||X|| Relative perturbation bounds
-22
Logarithm of compared results
-23
-24
-25
-26
-27
-28
5
10
15
20
25
30
35
40
45
50
n
Fig. 1. Logarithm of
kDXkF kXkF
and relative perturbation bounds C rel ðXÞ .rel with different values of n for Example 5.2 ða ¼ 0Þ.
726
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728 -27.5
-28 || X||/||X|| Relative perturbation bounds Logarithm of compared results
-28.5
-29
-29.5
-30
-30.5
-31
5
10
15
20
25
30
35
40
45
50
n
Fig. 2. Logarithm of
kDXkF kXkF
and relative perturbation bounds C rel ðXÞ .rel with different values of n for Example 5.2 ða ¼ 5Þ.
Table 4 Backward errors and residual bounds for Example 5.2 with different values of n (a ¼ 0). e gabs ð XÞ
n 5
e grel ð XÞ 13
e g abs ð XÞ 13
ek kX F
kXe X kF
12
1.04
2:34 1012
bResB 13
ke X kF
7:39 10
2:62 10
2:34 10
15
8:43 1013
1:02 1014
8:43 1013
1:86 1013
7:51 1012
1.02
7:50 1012
25
8:99 10
13
15
13
8:99 10
13
1:41 10
1:78 10
12
1.02
1:79 1012
35
9:22 1013
3:12 1015
9:23 1013
1:12 1013
1:81 1012
0.97
1:80 1012
45
13
15
13
14
12
1.01
8:92 1012
e k Xk F
kXe X kF
13
0.28
1:93 1013
14
9:94 1014
7:35 10
9:32 10
4:73 10
14
e g rel ð XÞ
5:13 10 2:19 10
9:33 10
9:79 10
8:89 10
Table 5 Backward errors and residual bounds for Example 5.2 with different values of n (a ¼ 5). n 10
e gabs ð XÞ
grel ð Xe Þ 13
9:81 10
e g abs ð XÞ 14
1:45 10
g rel ð Xe Þ 13
9:81 10
3:02 10
bResB 13
1:92 10
ke X kF
20
8:94 10
5:66 10
8:94 10
1:47 10
9:94 10
0.40
30
9:29 1013
3:47 1015
9:29 1013
1:25 1013
7:83 1014
0.47
7:83 1014
40
13
9:29 10
15
2:31 10
13
9:29 10
14
9:52 10
14
6:45 10
0.52
6:48 1014
50
9:22 1013
1:68 1015
9:22 1013
8:05 1014
5:52 1014
0.55
5:57 1014
13
15
13
13
backward errors are small enough to reveal the computing process of the approximate solution is backward stable. Moreover, from the results listed in Table 4 we observe that in the case of a ¼ 0; n ¼ 5; 15; 25; 35; 45, we have 11 e Xk 6 k X ek b kX : ResB < 10 F F
ð5:3Þ
e approximate the exact solutions X up to about 11 significant digits in the (5.3) implies that the computed solutions X sense of the Frobenius norm. Example 5.3 [11]. Consider NARE (1.1), in which
A ¼ D ¼ tridiagðIm ; T; Im Þ 2 Rnn
727
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728 Table 6 e Þ; g ð XÞ e and residual bound b Relative perturbation bound C rel ðXÞ .rel , backward errors grel ð X ResB with different values of n for Example 5.3 (n ¼ 0:2). rel kXe X kF
C rel ðXÞ .rel
e grel ð XÞ
e g rel ð XÞ
bResB
e k Xk F
4
2:9682 1011
1:9479 1010
4:5950 1014
2:6778 1012
2:9682 1011
0.0024
9
3:7273 1011
3:6172 1010
8:6259 1014
4:8496 1012
3:7273 1010
0.0059
16
4:9768 1011
6:1757 1010
1:4567 1013
8:0389 1012
4:9768 1011
0.0112
25
7:3295 1011
9:7197 1010
2:4319 1013
1:3253 1011
7:3294 1011
0.0184
36
8:9407 1011
1:4576 109
3:2104 1013
1:7227 1011
8:9407 1011
0.0274
n
ke X kF
are block tridiagonal matrices with n ¼ m2 , where
T ¼ tridiag 1; 4 þ
150 ðm þ 1Þ2
! ; 1
2 Rmm :
And matrices B and C are taken as following tridiagonal matrices
B¼
1 tridiagð1; 2; 1Þ 2 Rnn ; 50
C ¼ nB;
where n is a positive constant. The perturbations in the coefficient matrices are randomly generated ð0; 1Þ matrices by using MATLAB function randðnÞ, and DA; DB; DC; DD are taken the same as (5.2). In Table 6, while n ¼ 0:2 and size of matrices n is taken to be 4,9,16,25,36, we kXe Xk compare kXk F with the relative perturbation bound C rel ðXÞ .rel , the residual bound bResB , the backward relative errors F
e and g ð X e Þ, respectively. The results show that the relative perturbation bound C rel ðXÞ . is still sharp. We also see grel ð XÞ rel rel from Table 6 that the backward errors are small enough to reveal the computing process of the approximate solution is backward stable. Moreover, from the results listed in Table 6 we observe that in the case of n ¼ 0:2; n ¼ 4; 9; 16; 25; 36, we have 11 e Xk 6 k X ek b kX ; ResB < 10 F F
ð5:4Þ
e approximate the exact solutions X up to about 11 significant digits in the sense which implies that the computed solutions X of the Frobenius norm.
6. Conclusions In this paper, we discuss the condition number and two kinds of backward errors of nonsymmetric algebraic Riccati equation. The explicit expression of condition number is derived in a uniform manner. Two kinds of backward errors of the minimal nonnegative solution are evaluated. Numerical examples show that our estimate is fairly sharp. Condition number can reflect the sensitivity of the approximate solution when the coefficient matrices A; B; C; D have small perturbations. And the backward errors can test the stability of the computing process of the approximate solution. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2014.06.036. References [1] Z.-Z. Bai, X.-X. Guo, S.-F. Xu, Alternately linearized implicit iteration methods for the minimal nonnegative solutions of nonsymmetric algebraic Riccati equations, Numer. Linear Algebra Appl. 13 (2006) 655–674. [2] Z.-Z. Bai, Y.-H. Gao, L.-Z. Lu, Fast iterative schemes for nonsymmetric algebraic Riccati equations arising from transport theory, SIAM J. Sci. Comput. 30 (2008) 804–818. [3] D.A. Bini, B. Iannazzo, F. Polon, A fast Newton’s method for a nonsymmetric algebraic Riccati equation, SIAM J. Matrix Anal. Appl. 30 (2008) 276–290. [4] G. Freiling, A survey of nonsymmetric Riccati equations, Linear Algebra Appl. 351–352 (2002) 243–270. [5] Y.-H. Gao, Z.-Z. Bai, On inexact Newton methods based on doubling iteration scheme for nonsymmetric algebraic Riccati equations, Numer. Linear Algebra Appl. 18 (2011) 325–341. [6] C.-H. Guo, A.J. Laub, On the iterative solution of a class of nonsymmetric algebraic Riccati equations, SIAM J. Matrix Anal. Appl. 22 (2000) 376–391. [7] C.-H. Guo, Nonsymmetric algebraic Riccati equations and Wiener–Hopf factorization for M-matrices, SIAM J. Matrix Anal. Appl. 23 (2001) 225–242. [8] C.-H. Guo, A new class of nonsymmetric algebraic Riccati equations, Linear Algebra Appl. 426 (2007) 636–649. [9] C.-H. Guo, B. Iannazzo, B. Meini, On the doubling algorithm for a (shifted) nonsymmetric algebraic Riccati equation, SIAM J. Matrix Anal. Appl. 29 (2007) 1083–1100.
728
L.-d. Liu, S.-f. Xu / Applied Mathematics and Computation 242 (2014) 716–728
[10] X.-X. Guo, Z.-Z. Bai, On the minimal nonnegative solution of nonsymmetric algebraic Riccati equation, J. Comput. Math. 23 (2005) 305–320. [11] X.-X. Guo, W.-W. Lin, S.-F. Xu, A structure-preserving doubling algorithm for nonsymmetric algebraic Riccati equation, Numer. Math. 103 (2006) 393– 412. [12] N.J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM Society for Industrial and Applied Mathematics, Philadelphia, 2002. [13] J. Juang, W.-W. Lin, Nonsymmetric algebraic Riccati equations and Hamiltonian-like matrices, SIAM J. Matrix Anal. Appl. 20 (1999) 228–243. [14] T.-X. Li, E.K. 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. [15] W.-W. Lin, J.-G. Sun, Perturbation analysis of the periodic discrete-time algebraic Riccati equation, SIAM J. Matrix Anal. Appl. 24 (2002) 411–438. [16] W.-W. Lin, S.-F. Xu, Convergence analysis of structure-preserving doubling algorithms for Riccati-type matrix equations, SIAM.J. Matrix Anal. Appl. 28 (2006) 26–39. [17] L.-D. Liu, Mixed and componentwise condition numbers of nonsymmetric algebraic Riccati equation, Appl. Math. Comput. 218 (2012) 7595–7601. [18] J.R. Rice, A theory of condition, SIAM J. Numer. Anal. 3 (1966) 287–310. [19] J.-G. Sun, Perturbation Analysis of Matrix, second ed., Science Press, Beijing, 2001 (in Chinese). [20] J.-G. Sun, Backward error for the discrete-time algebraic Riccati equation, Linear Algebra Appl. 259 (1997) 183–208. [21] J.-G. Sun, Sensitivity analysis of the discrete-time algebraic Riccati equation, Linear Algebra Appl. 275–276 (1998) 595–615. [22] J.-G. Sun, Condition numbers of algebraic Riccati equations in the Frobenius norm, Linear Algebra Appl. 350 (2002) 237–261. [23] J.-G. Sun, Perturbation analysis of algebraic Riccati equations, Technical Report UMINF 02.03, Department of Computing Science, Umeå University, 2002. [24] J.-G. Sun, Backward perturbation analysis of the periodic discrete-time algebraic Riccati equation, SIAM J. Matrix Anal. Appl. 26 (2004) 1–19. [25] J.H. Wilkinson, The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, England, 1965. [26] S.-F. Xu, Matrix Computation in Control System, Higher Education Press, Beijing, 2011 (in Chinese). [27] J.-G. Xue, S.-F. Xu, R.-C. Li, Accurate solutions of M-matrix Sylvester equations, Numer. Math. (120) (2012) 639–670. [28] J.-G. Xue, S.-F. Xu, R.-C. Li, Accurate solutions of M-matrix algebraic Riccati equations, Numer. Math. (120) (2012) 671–700. [29] L.-M. Zhou, Y.-Q. Lin, Y.-M. Wei, S.-Z. Qiao, Perturbation analysis and condition numbers of symmetric algebraic Riccati Equations, Automatica 45 (2009) 1005–1011.