Nonrecursive solution for the discrete algebraic Riccati equation and X + A*X-1A=L

Abstract In this paper, we present two new algebraic algorithms for the solution of the discrete algebraic Riccati equation. The first algorithm requires the nonsingularity of the transition matrix and is based on the solution of a standard eigenvalue problem for a new symplectic matrix; the proposed algorithm computes the extreme solutions of the discrete algebraic Riccati equation. The second algorithm solves the Riccati equation without the assumption of the nonsingularity of the transition matrix; the proposed algorithm is based on the solution of the matrix equation X + A*X-1A=L, where A is a singular matrix and L a positive definite matrix.


Introduction
The well known discrete algebraic Riccati equation (DARE) is formed where A; Q and X are n n complex matrices, R is an m m matrix, B is an n m; and it is assumed that Q; R are Hermitian and positive definite matrices; here A denotes the complex conjugate A of the transpose A T of the transition matrix A: The discrete algebraic Riccati equation (DARE) has attracted enormous attention. One of the more important motivations for the study of DARE was the recognition of the work of Kalman [15] on filtering and prediction problem of discrete type. The discrete time Kalman filter [4,15] is the well known algorithm that solves the filtering problem, that solves the corresponding DARE emanating from Kalman filter. The solution is required to be a unique Hermitian positive semidefinite matrix (or positive definite, if there exists). The existence of this solution, X R , depends on the modulus of all the eigenvalues of matrix A B.RCB X R B/ 1 B X R A: In view of the importance of the Riccati equation, there exists considerable literature on its algebraic solutions [4,16,22] and recursive solutions [4,6,11,18] concerning per step or doubling algorithms.
A detailed discussion of the solution theory for the case where A is a nonsingular complex matrix is given in [16,20], whereas the case of singular A is treated in [16] (see and references therein). The algebraic solution method of (1) is closely connected to the standard eigenvalue problem [22] or the generalized eigenvalue problem [9,12,17,19] for a 2n 2n real matrix, which depends on the coefficients matrices A; B; Q; R of (1). The usual solution method is the classical eigenvector approach and the case of the nonsingularity of A described in [22]; the Schur vector method is established in [9,17,19] without the assumption of nonsingularity of A.
The nonlinear matrix equation, which arises in various research areas including control theory, stochastic filtering and statistics, (see [1,8] and the references therein), is formed where A is an n n matrix and L is an n n positive definite matrix. It is well known [1,3,8] that, the existence of the positive definite solution of the equation (2) depends on the numerical radius of matrix L 1=2 AL 1=2 : In particular, it is required where F .A/ D fx Ax W for every x 2 C n with x x D 1g is defined the numerical range (also known as the field of values) of A; and r.A/ D maxfjzj W z 2 F .A/g is defined the numerical radius of A; [8,10]. Also, if the nonlinear matrix equation (2) has a unique positive definite solution X; then there exist minimal and maximal solutions X mi n and X max ; respectively, such that 0 < X mi n Ä X Ä X max for any positive definite solution X: We will refer to X mi n and X max ; as the extreme solutions of (2). Here, when X and Y are Hermitian matrices, X > Y means that X Y is a positive definite matrix denoted X Y > 0 and X Y means that X Y is a positive semidefinite matrix denoted X Y 0: The theoretical properties of the positive definite solution X max of (2) have been investigated by many authors (see [8,14,21] and the references therein) and the available methods for computing X max are algebraic (nonrecursive) [1][2][3] and recursive algorithms concerning per step or doubling [5,7,8,13,14].
The positive definite solutions X R ; X max of the nonlinear equations (1) and (2) are related and the solutions methods are combined [1,2,7,13,14].
Although for applications the real case (i.e., when all the matrices A; B; Q; R are real and real solution matrices X are to be found) is especially important, we consider here the general (i.e., complex) case as well as the real case. The objective of this paper is to develop a new algebraic method for the computation of the solution of the DARE independent on the nonsingularity of the transition matrix A and link the method to the solution of the matrix equation X C A X 1 A D L in (2). The proposed algebraic method for the computation of the solution of the DARE (1) is based on the standard eigenvalue problem for the 2n 2n matrix which depends only on the coefficients A; B; Q; R of (1), (see the developed theory in terms of Lagrangian subspaces in [16, chapters 12-15]). The key idea of this paper is to construct the solution of the DARE in (1) using a basis for the eigenspace of T: The paper is organized as follows. In Section 2, the classical algebraic Riccati equation method for nonsingular transition matrix is presented using a new symplectic matrix in order to be computed the extreme solutions of the DARE. Section 3 contains the relationship between the discrete algebraic Riccati equations and the matrix equation In Section 4, the new algebraic Riccati equation method for singular transition matrix is developed and the proposed algorithm is implemented in two examples. Finally, concluding remarks are given in Section 5.

Nonrecursive algebraic solution for the DARE with nonsingular transition matrix
Here, motivated by [2,22], an algebraic Riccati equation solution method is presented, which is based on the use of a basis for the eigenspace of T in (4). The algebraic Riccati equation solution method requires the transition matrix A to be nonsingular, so that the matrix T to be well defined. At this point we make the assumption that A is a nonsingular matrix, which will be removed in the next section.
Using the 2n 2n Hermitian matrix where I n denotes the n n identity matrix, a straightforward calculation shows that A matrix T is called symplectic, whenever equality (6) is verified. It is obvious by (6) that a symplectic matrix is a nonsingular matrix, thus all its eigenvalues are non-zero, that is 0 … .T /; where .T / denotes the spectrum of T: Lemma 2.1. Let A be a nonsingular matrix, Q; R; be positive definite matrices and T a symplectic matrix as in (4).
If is an eigenvalue of T; then . / 1 2 .T /: Proof. Since T is nonsingular matrix, it is easy to verify that its inverse matrix is formulated as Assume that is an eigenvalue of T; and h x 1 x 2 i T is the corresponding eigenvector, thus Since is an eigenvalue of T T ; it follows that 1 2 ..T 1 / T /I hence from (7) we obtain x 2 x 1 # and the conjugate of the above equality yields Since the linear systems in (8) and (9) coincide for and . / 1 ; the proof is completed.
According to Lemma 2.1 and supposing that the eigenvalues of the symplectic matrix T in (4) are distinct, T can be written as where is the diagonal matrix, which contains the eigenvalues of matrix T that lie outside the unit disk and W is the matrix of the corresponding eigenvectors of T ; its partition into four n n blocks is denoted as A matrix A is said to be stable if all eigenvalues of A lie in the open unit disk, i.e. j i .A/j < 1 for i D 1; : : : ; n; the pair .A; B/ is called stabilizable, if there exists an m n matrix K such that all the eigenvalues of A C BK lie in the open unit disk, and a Hermitian solution X of the DARE is called stabilizing if all the eigenvalues of A B.R C B XB/ 1 B XA lie in the open unit disk.
Since the existence of a Hermitian solution of (1) guarantees the uniqueness of the maximal (Hermitian) solution X R , [16, Theorem 13.1.1], for which holds X R X for all Hermitian solutions X of (1), firstly we are interested to find one.
Assuming that .A; B/ is a stabilizable pair, the DARE (1) is known to have a unique X R positive semidefinite solution and all the eigenvalues of the matrix O [16,Corollary 13.1.2]). If, additional, .A ; Q/ is a stabilizable pair, then O A 1 is stable; X R is positive definite, when .A ; Q/ is a controllable pair, (see, [16,Theorem 13.1.3]). There are, of course, many other solutions of (1) but we are interested in computing the positive semidefinite one, (or the positive definite one, if it exists).
In the following, the maximal solution of (1) is given in terms of the eigenvectors of matrix T; which is achieved by the property of stability.
is the unique maximal Hermitian solution of (1) and X R is positive definite matrix.
Proof. Rewriting the equality in (10) as where T; W are given by (4), (11), respectively, the following equalities are derived: The nonsingularity of W 12 allows us to write the equality (15) in the form and substituting the above equality in (16) arises Multiplying the above equality by W 1 12 (on the right) and A (on the left) we obtain The nonsingularity of A allows us to write the above equality in the form We claim that I n C W 22 W 1 12 BR 1 B is nonsingular; indeed, let the row vector y be such that Multiplying (17) Multiplying (17) Hence from (19) we obtain and so X D W 22 W 1 12 satisfies (1). Since Q; R are positive definite matrices and the pairs .A; B/ and .A ; Q/ are stabilizable, all the hypotheses of Theorem 13.1.3 in [16] are satisfied, hence the existence of the maximal Hermitian solution i) It is obvious from (4) that the algebraic method for the computation of the maximal solution (12) of the DARE (1) cannot be used, when the matrix A is singular. ii) The assumptions of Theorem 2.2 guarantee that T has no eigenvalues on the unit disk [19,Theorem 3]. iii) If the block matrices W 11 ; W 21 are nonsingular, using (13)- (14) and following the same statements as in the proof of Theorem 2.2, the minimal Hermitian solution X of (1) is given by In the real case, the DARE (1) is formed where the coefficients A; B; Q; R are real matrices of sizes n n; n m; n n; and m m; respectively, with Q; R symmetric and Q; R > 0; X is an n n real symmetric matrix to be found. The 2n 2n matrix requires the nonsingularity of the matrix A in (22) in order to be well defined and using the 2n 2n real matrix it is proved that T is a symplectic matrix. Thus, T can be written as where is the diagonal matrix, which contains the eigenvalues of matrix T that lie outside the unit disk and W is the matrix, which contains the corresponding eigenvectors of T in (22). Consider the associated partition of W as in (11), then, the unique maximal real symmetric solution of (21) is related to the eigenvectors of matrix T and formulated as in (12), i.e., X R D W 22 W 1 12 : i) It is obvious from (22) that the above algebraic method cannot be used, when the matrix A is singular.
ii) The maximal real solution of DARE in (21) is formulated through the eigenvectors of a symplectic matrix˚in [2,22], which has the same form as in (12), (see, [2, subsection 2.2.1]). It is expected, since T in (22) is the inverse matrix of˚in [2,22], hence using the factorization of T in (24), we can write where W is given by (11) and is the diagonal matrix with the eigenvalues of matrix˚, that lie outside the unit disk.

DARE solution via the solution of the matrix equation
Assuming that X is a nonsingular matrix and since R is positive definite, applying the matrix inversion lemma in the DARE (1), the following equivalent DARE is derived Under the assumptions of Theorem 2.2 it is obvious that X R in (12) is a nonsingular matrix, hence, replacing X by X R , the DARE in (1) is equivalent to the DARE in (25). Now, we are able to derive a nonlinear matrix equation as in (2), that follows from the DARE in (25). Indeed, setting in (25) The inverse matrix of Q C A ˘C BR 1 B 1 A is derived by the matrix inversion lemma in the above equation Setting the equation (27) can be written as which is formulated as (2), where the coefficients A; L are Thus, we are able to compute the solution of the DARE (25) via the solution of the matrix equation (2) with the coefficients A; L in (30), which are calculated in terms of the parameters A; B; Q; R of the DARE.
Recall that the solvability of the matrix equation (2) is related to the numerical radius of L 1=2 AL 1=2 and the relation is formulated in (3). Considering that the matrix equation (2) is solvable and X max D X is the unique maximal solution, [8], the unique maximal solution X R of the DARE (25) can be computed in terms of X max . In fact, using (26) and (28) we have

Matrix equation solution via the solution of the DARE
Consider that A in the matrix equation (2) is nonsingular and r.L 1=2 AL 1=2 / Ä 1 2 : Working as in [1,2], we are able to derive a DARE, which is implied by (2). Indeed, since the matrix equation (2) can be written as X D L A X 1 A; the substitution of the matrix X with L A X 1 A and the nonsingularity of A yield the following equality The above related DARE is formed as (25) and the associated coefficients are The assumption for the numerical radius guarantees the existence of the solutions of the matrix equation (2), thus now we are able to compute the unique maximal positive definite solution X max of (2), which coincides with the unique positive definite solution X R of the above related DARE, i.e., Since the DARE (1) is equivalent to DARE (25), the solution of the above related DARE can be derived using the algebraic proposed method in Section 2. More specifically, using the following symplectic matrix and (10)- (11), the maximal solution X R of the related DARE is given in terms of the eigenvectors of matrix T as in (12). Finally, the maximal positive definite solution of (2) is computed from (12) and (32). The boundary of the numerical range L 1=2 AL 1=2 is illustrated in Figure 1, from which the numerical radius appears as r.L 1=2 AL 1=2 / 0:39: Since r.L 1=2 AL 1=2 / Ä 1 2 ; the given matrix equation has a unique positive definite solution X max , which is computed by (32). The symplectic matrix T; formulated in (33) of the paragraph 3.2, is T D : Considering the partition of above W as in (11) easily it is verified that the matrices W 11 ; W 12 ; W 21 and W 22 are nonsingular, hence the maximal solution X max follows from (32) and (12) and it is equal to

DARE solution method for singular transition matrix
In this section, we are going to develop a new algebraic method for the solution of the DARE (1) without the restriction of the nonsingularity of the transition matrix A as it has been required in the previous section; the basic idea is to solve the DARE via the maximal solution of the matrix equation (2). In particular, the maximal solution X max of (2) is computed by the implementation of the below algorithm and in the following the maximal solution X R is computed by the methodology of the paragraph 3.1. Firstly, the given DARE is transformed into the equivalent DARE (25), from which the matrix equation (29) follows. From (30), it is clear that the singularity of A implies the singularity of A; consequently it is needed to be developed a methodology in order to solve the matrix equation (2) with A singular.
In the following, working as in [1,2] multiplying the both sides of (2) with the nonsingular matrix L 1=2 arises the equation L 1=2 XL 1=2 C L 1=2 A X 1 AL 1=2 D I n ; in which setting Thus, it is obvious that the computation of the solution of (2) is achieved through the solving of the special matrix equation (36) and the solutions are related as in (34). Further, by (35) it is evident that the nonsingularity of A is equivalent to the nonsingularity of C . Based on the ideas of the algorithm mentioned in [1,7], we proposed an algorithm, which achieves the solution of (36) for singular matrix C: Using a Schur factorization of C in (35), there exists an n n unitary matrix U such that where C is an .n 1/ .n 1/ matrix, 0 is the .n 1/ 1 zero-vector and c is an 1 .n 1/ row vector. Consider that there exists an .n 1/ .n 1/ nonsingular matrix Z such that Substituting the matrices C; Y by (37), (38) in (36) and using the properties of the matrices U; Z; we derive: The above matrix equation leads to where L D I n 1 c c: Note that, c c is a Hermitian matrix, hence the matrix equation (39) is formed as (2) if and only if L > 0: Thus, we distinguish the following cases: i) If the matrix L is not positive definite, then the algorithm is finished. ii) If L > 0 and r.L 1=2 CL 1=2 / Ä 1 2 the following cases are examined: a. If C is the zero matrix, then the matrix equation (39) yields Z D L D I n 1 c c: The maximal solution of (2) follows from (38) and (34).
b. If C ¤ O; and detC ¤ 0; then the matrix equation (39) is formed as (2). Thus, the solving of (39) is achieved using the methodology in the paragraph 3.2 for the computation of the maximal solution Z: The solution of (2) follows from (38) and (34). c. If C ¤ O; and det C D 0; then multiplying both sides of (39) with the nonsingular matrix L 1=2 ; the matrix equation (39) is transformed in the equation of the type (36), and the algorithm is repeated.
Finally, applying the above algebraic algorithm the unique maximal solution X max of the matrix equation in (2) with A singular can be computed. Using the solution X max the maximal solution X R of the DARE (1), whose A is singular, arises from (31). In the following the above algorithm is implemented in two examples with A singular.
The boundary of the numerical range L 1=2 AL 1=2 is illustrated in Figure 2, from which the numerical radius appears as r.L 1=2 AL 1=2 / 0:45:  The DARE (1) with real coefficients founds an application of the associated Riccati equation emanating from Kalman filter, which arises from time invariant system with standard coefficients F; H; Q; R and formed as the DARE in (21). For time invariant system, it is well known [4] that if the signal process model is asymptotically stable (i.e. all eigenvalues of the transition matrix lie in the unit disk), then there exists a steady state value P of the prediction error covariance matrix; P can be calculated by solving the corresponding discrete time Riccati equation emanating from Kalman filter where F is the n n transition matrix and H is the m n output matrix. Matrix Q is the n n covariance matrix for the white plant noise process and R is the m m covariance matrix for the white measurement noise process and it is convenient to assume throughout that Q; R > 0: It is evident that setting F D A T and H D B T , the above DARE follows the DARE in (21)   The boundary of the numerical range L 1=2 AL 1=2 is illustrated in Figure 3, from which the numerical radius appears as r.L 1=2 AL 1=2 / 0:17: Finally, by (31) the maximal solution of the given DARE is computed

Conclusions
The standard algebraic solution of the discrete algebraic Riccati equation (DARE) requires the nonsingularity of the transition matrix. In this paper two new algebraic algorithms for solution of the discrete algebraic Riccati equation were proposed: one for nonsingular transition matrix and the other for singular transition matrix.
Concerning the nonsingular transition matrix case, the proposed algorithm requires the nonsingularity of the transition matrix and is based on the solution of a standard eigenvalue problem for a new symplectic matrix. The maximal and minimal solutions of the DARE are derived. Concerning the singular transition matrix case, the proposed algorithm solves the discrete algebraic Riccati equation through the solutions of the matrix equation X C A X 1 A D L for suitable A; L matrices, which are related to the matrices of the DARE. Note that the singularity of the transition matrix is equivalent to the singularity of A: The proposed method is especially important for solving the DARE emanating from Kalman filter, which arises from a time invariant system with standard real coefficients. The proposed algebraic method faces successfully the singular transition matrix case and computes the maximal solution, which is equivalent to the steady state prediction error covariance matrix.