Condition number and backward errors of nonsymmetric algebraic Riccati equation

Condition number and backward errors of nonsymmetric algebraic Riccati equation

Applied Mathematics and Computation 242 (2014) 716–728 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepag...

530KB Sizes 0 Downloads 50 Views

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





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





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



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.