On the Iterative Method for the System of Nonlinear Matrix Equations

Abstract

The positive definite solutions for the system of nonlinear matrix equations 𝑋+𝐴^∗𝑌^−𝑛𝐴=𝐼,𝑌+𝐵^∗𝑋^−𝑚𝐵=𝐼 are considered, where 𝑛,𝑚 are two positive integers and A, B are nonsingular complex matrices. Some sufficient conditions for the existence of positive definite solutions for the system are derived. Under some conditions, an iterative algorithm for computing the positive definite solutions for the system is proposed. Also, the estimation of the error is obtained. Finally, some numerical examples are given to show the efficiency of the proposed iterative algorithm.

1. Introduction

Linear and nonlinear matrix equations have been widely used for solving many problems in several areas such as control theory, optimal control, optimization control, stability theory, communication system, dynamic programming, signal processing, and stochastic filtering and statistics, [1–3]. Many authors studied the existence of solutions for several classes of the matrix equations (see, e.g., [4–14]), in particular, Lyapunov matrix equation [15], Sylvester matrix equations [11, 14], algebraic Riccati equations [3], some special case of linear and nonlinear matrix equations [16–21], and coupled matrix equations [22–24].

In recent years, many types of algebraic Riccati equations have been the subject of great activity, the aim being to achieve a fast and reliable algorithm that generates numerical positive definite solutions.

In this paper, we will consider the system (Sys.) of nonlinear matrix equations that can be expressed in the form: 𝑋+𝐴^∗𝑌^−𝑛𝐴=𝐼,𝑌+𝐵^∗𝑋^−𝑚𝐵=𝐼, (1) where 𝑛,𝑚 are two positive integers,  𝑋,𝑌 are 𝑟×𝑟 unknown matrices,  𝐼 is the 𝑟×𝑟 identity matrix, and  𝐴,𝐵 are nonsingular matrices. All matrices are defined over the complex field. The system of nonlinear matrix equations with the form of (1) is a special case of the system of algebraic discrete-type Riccati equations of the form: 𝑋_𝑖=𝑉^∗_𝑖𝑋_𝑖𝑉_𝑖+𝑄_𝑖−𝑉^∗_𝑖𝑋_𝑖𝐵_𝑖+𝐴_𝑖𝑅_𝑖+𝐵^∗_𝑖𝑋_𝑖𝐵_𝑖⁻¹𝐵^∗_𝑖𝑋_𝑖𝑉_𝑖+𝐴^∗_𝑖, (2) where 𝑖=1,2,…,𝑘, [2, 3]. The efficient numerical solutions for some special case of the system (2) have been extensively studied by several authors [4–10, 22–26]. For example, Mukaidani [22] proposed a new algorithm for solving cross-coupled sign-indefinite algebraic Riccati equations for weakly coupled large-scale systems, while in [4, 5] Al-Dubiban has studied special cases of Sys. (2) by obtained sufficient conditions for the existence of positive definite solutions for the systems and proposed iterative algorithms to calculate the solutions. In [10], Davies proposed upper bounds for the sum of the maximal eigenvalues of the solutions of the continuous and discrete coupled algebraic Riccati equations. In [25], Ivanov has studied a set of discrete-time coupled algebraic Riccati equations which arise in quadratic optimal control and proposed two iterations for computing a symmetric solution of this system.

In this paper, we derive the sufficient conditions of the existence of solutions for the Sys. (1). We introduce an iterative algorithm to obtain the positive definite solutions of Sys. (1). We discuss the convergence of this iterative algorithm. Finally, some numerical examples are given to illustrate the efficiency for suggested algorithm.

The following notations are used throughout the rest of the paper. The notation 𝐴≥0  (𝐴>0) means that 𝐴 is positive semidefinite (positive definite), 𝐴^⋆ denotes the complex conjugate transpose of 𝐴, and 𝐼 is the identity matrix. Moreover, 𝐴≥𝐵  (𝐴>𝐵) is used as a different notation for 𝐴−𝐵≥0  (𝐴−𝐵>0). We denote by 𝜌(𝐴) the spectral radius of  𝐴;  𝜆_𝑟(𝑋),𝜇_𝑟(𝑌) represent the eigenvalues of 𝑋 and 𝑌, respectively. The norm used in this paper is the spectral norm of the matrix 𝐴; that is,  ‖𝐴‖=𝜌(𝐴𝐴^⋆) unless otherwise noted.

2. Main Theorems

In this section, we will introduce an iterative algorithm which is applicable for computing the positive definite solutions of the Sys. (1). We start with some results which will be used throughout this paper.

Lemma 1 (see [27, 28]). If 𝑃>𝑄>0 (or 𝑃≥𝑄>0), then 𝑃^𝛼>𝑄^𝛼 (or 𝑃^𝛼≥𝑄^𝛼>0) for all 𝛼∈(0,1], and 𝑃^𝛼<𝑄^𝛼 (or 0<𝑃^𝛼≤𝑄^𝛼) for all 𝛼∈[−1,0).

Theorem 2 (see [29]). Let the matrices 𝑃,𝑄, and 𝑅 be positive definite 𝑟×𝑟 matrices, such that the integral ^∞₀𝑒^𝑃𝑡𝑅𝑒^𝑄𝑡𝑑𝑡 (3) exists and lim_𝑡→∞𝑒^𝑃𝑡𝑅𝑒^𝑄𝑡=0; (4) then the matrix 𝑍=−^∞₀𝑒^𝑃𝑡𝑅𝑒^𝑄𝑡𝑑𝑡 (5) is the solution of the matrix equation: 𝑃𝑍+𝑍𝑄=𝑅. (6)

The solution of Sys. (1) can be found by the following iterative algorithm.

Algorithm 3. 𝑋₀=𝑌₀=𝛽𝐼,𝑋_𝑠+1=𝐵𝐼−𝑌_𝑠⁻¹𝐵^∗1/𝑚,𝑌_𝑠+1=𝐴𝐼−𝑋_𝑠⁻¹𝐴^∗1/𝑛,𝑠=0,1,2,…. (7)

Theorem 4. If there exist numbers 𝛼,𝛽 satisfying 0<𝛼<𝛽≤min{𝑛/(𝑛+1),𝑚/(𝑚+1)}, and the following conditions hold:(i)𝛼^𝑛(1−𝛼)𝐼<𝐴𝐴^∗<𝛽^𝑛(1−𝛽)𝐼, (ii)𝛼^𝑚(1−𝛼)𝐼<𝐵𝐵^∗<𝛽^𝑚(1−𝛽)𝐼, then the sequences {𝑋_𝑠},{𝑌_𝑠} defined by Algorithm 3 converge to a positive definite solution (𝑋,𝑌) of Sys. (1).

Proof. From Algorithm 3, we get 𝑋₁=𝐵𝐼−𝛽𝐼⁻¹𝐵^∗1/𝑚=11−𝛽𝐵𝐵^∗1/𝑚<𝛽^𝑚1−𝛽𝐼1−𝛽^1/𝑚=𝛽𝐼=𝑋₀. (8) Also, we have 𝑋₁=[𝐵(𝐼−𝛽𝐼)⁻¹𝐵^∗]^1/𝑚=11−𝛽𝐵𝐵^∗1/𝑚>1(1−𝛼)𝐵𝐵^∗1/𝑚>𝛼^𝑚(1−𝛼)()𝐼1−𝛼^1/𝑚=𝛼𝐼. (9) That is,  𝑋₀>𝑋₁>𝛼𝐼, similarly we get 𝑌₁=𝐴𝐼−𝛽𝐼⁻¹𝐴^∗1/𝑛=11−𝛽𝐴𝐴^∗1/𝑛<𝛽^𝑛1−𝛽𝐼1−𝛽^1/𝑛=𝛽𝐼=𝑌₀. (10) Also, we have 𝑌₁=𝐴𝐼−𝛽𝐼⁻¹𝐴^∗1/𝑛=11−𝛽𝐴𝐴^∗1/𝑛>1(1−𝛼)𝐴𝐴^∗1/𝑛>𝛼^𝑛(1−𝛼)()𝐼1−𝛼^1/𝑛=𝛼𝐼. (11) That is,  𝑌₀>𝑌₁>𝛼𝐼.
Suppose that 𝑋_𝑠−1>𝑋_𝑠>𝛼𝐼,𝑌_𝑠−1>𝑌_𝑠>𝛼𝐼. (12) Now, we will prove that 𝑋_𝑠>𝑋_𝑠+1>𝛼𝐼 and 𝑌_𝑠>𝑌_𝑠+1>𝛼𝐼.
By using the inequalities (12), we have 𝑋_𝑠+1=𝐵𝐼−𝑌_𝑠⁻¹𝐵^∗1/𝑚<𝐵𝐼−𝑌_𝑠−1⁻¹𝐵^∗1/𝑚=𝑋_𝑠. (13) Also, we have 𝑋_𝑠+1=𝐵𝐼−𝑌_𝑠⁻¹𝐵^∗1/𝑚>𝐵(𝐼−𝛼𝐼)⁻¹𝐵^∗1/𝑚=1(1−𝛼)𝐵𝐵^∗1/𝑚>𝛼^𝑚(1−𝛼)()𝐼1−𝛼^1/𝑚=𝛼𝐼. (14) Similarly, we get 𝑌_𝑠+1=𝐴𝐼−𝑋_𝑠⁻¹𝐴^∗1/𝑛<𝐴𝐼−𝑋_𝑠−1⁻¹𝐴^∗1/𝑛=𝑌_𝑠. (15) Also, we have 𝑌_𝑠+1=𝐴𝐼−𝑋_𝑠⁻¹𝐴^∗1/𝑛>𝐴(𝐼−𝛼𝐼)⁻¹𝐴^∗1/𝑛=1(1−𝛼)𝐴𝐴^∗1/𝑛>𝛼^𝑛(1−𝛼)()𝐼1−𝛼^1/𝑛=𝛼𝐼. (16) Therefore, the inequalities (12) are true for all 𝑠=1,2,3,….
Hence, the sequences {𝑋_𝑠},{𝑌_𝑠} are monotonically decreasing and bounded from below by the matrix 𝛼𝐼. Consequently, the sequences converge to a positive definite limit (𝑋,𝑌) which is a solution of Sys. (1).

Theorem 5. If there exist numbers 𝛼,𝛽 satisfying 0<𝛼<𝛽≤min{𝑛/(𝑛+1),𝑚/(𝑚+1)}, and the following conditions hold: (i)𝛼^𝑛(1−𝛼)𝐼<𝐴𝐴^∗<𝛽^𝑛(1−𝛽)𝐼, (ii)𝛼^𝑚(1−𝛼)𝐼<𝐵𝐵^∗<𝛽^𝑚(1−𝛽)𝐼, (iii)𝑞=𝛽^{(2𝑘+2𝑡)}/(2𝛼)^(𝑘+𝑡)(1−𝛽)²<1, where 𝑚=2^𝑘,𝑛=2^𝑡, then Sys. (1) has a positive definite solution (𝑋,𝑌) which satisfies 𝑋_𝑠−𝑋≤𝑞𝑋_𝑠−2−𝑋≤⋯≤𝑞^𝑠/2,𝑌𝛽−𝛼_𝑠−𝑌≤𝑞𝑌_𝑠−2−𝑌≤⋯≤𝑞^s/2𝛽−𝛼. (17)

Proof. From Theorem 4, the two sequences {𝑋_𝑠},{𝑌_𝑠} defined by Algorithm 3 are convergent to a positive definite solution (𝑋,𝑌) of Sys. (1). We compute the spectral norm of the matrices 𝑋_𝑠−𝑋,𝑌_𝑠−𝑌. For that, we have 𝑋_𝑠−𝑋=^𝑚𝐵(𝐼−𝑌_𝑠−1)⁻¹𝐵^∗−^𝑚𝐵(𝐼−𝑌)⁻¹𝐵^∗. (18) We denote 𝑃₁=𝐵𝐼−𝑌_𝑠−1⁻¹𝐵^∗,𝑄₁=𝐵(𝐼−𝑌)⁻¹𝐵^∗. (19) We use the following equality: ^𝑚𝑃₁^𝑚𝑃₁−^𝑚𝑄₁+^𝑚𝑃₁−^𝑚𝑄₁^𝑚𝑄₁=^𝑚/2𝑃₁−^𝑚/2𝑄₁. (20)Since 𝑌_𝑠>𝑌>0 for each 𝑠=0,1,2,…, then by using Lemma 1 we have the matrix 𝐷₁=^𝑚𝑃₁−^𝑚𝑄₁ being a positive definite solution of the matrix equation: ^𝑚𝑃₁𝐷₁+𝐷₁^𝑚𝑄₁=^𝑚/2𝑃₁−^𝑚/2𝑄₁. (21) According to Theorem 2, we have 𝐷₁=^∞₀𝑒^{−𝑚𝑃1𝑡𝑚/2}𝑃₁−^𝑚/2𝑄₁𝑒^{−𝑚𝑄1𝑡}𝑑𝑡. (22) Since ^𝑚𝑃₁,^𝑚𝑄₁ are positive definite matrices, then the integral (22) exists, and 𝑒^{−𝑚𝑃1𝑡𝑚/2}𝑃₁−^𝑚/2𝑄₁𝑒^{−𝑚𝑄1𝑡}→0,as𝑡→∞. (23) By using (18) and (22), we have 𝑋_𝑠−𝑋≤^∞₀^𝑚/2𝑃₁−^𝑚/2𝑄₁𝑒^{−𝑚𝑃1𝑡}e^{−𝑚𝑄1𝑡}𝑑𝑡. (24) However, 𝑌_𝑠>𝑌≥𝛼𝐼; hence, ^𝑚𝑃₁>𝛼𝐼,^𝑚𝑄₁≥𝛼𝐼. (25) Then, we have 𝑋_𝑠−𝑋≤^𝑚/2𝑃₁−^𝑚/2𝑄₁^∞₀𝑒^{−𝑚𝑃1𝑡}𝑒^{−𝑚𝑄1𝑡}≤𝑑𝑡_𝑚/2𝑃₁−^𝑚/2𝑄₁^∞₀𝑒^{−𝛼𝐼𝑡2}=𝑑𝑡_𝑚/2𝑃₁−^𝑚/2𝑄₁^∞₀𝑒^−2𝛼𝑡=1𝑑𝑡2𝛼^𝑚/2𝑃₁−^𝑚/2𝑄₁≤12𝛼^22𝑚/2𝑃₁−^2𝑚/2𝑄₁. (26) After 𝑘 times as above, we get 𝑋_𝑠≤1−𝑋2𝛼^{𝑘𝑘𝑚/2}𝑃₁−^𝑘𝑚/2𝑄₁. (27) Let 𝑚=2^𝑘 be special case, then we have 𝑋_𝑠≤1−𝑋2𝛼^𝑘𝑃₁−𝑄₁=12𝛼^𝑘𝐵𝐼−𝑌_𝑠−1⁻¹𝐵^∗−𝐵(𝐼−𝑌)⁻¹𝐵^∗=12𝛼^𝑘𝐵𝐼−𝑌_𝑠−1⁻¹𝑌_𝑠−1−𝑌(𝐼−𝑌)⁻¹𝐵^∗≤12𝛼^𝑘‖𝐵‖²𝐼−𝑌_𝑠−1⁻¹(I−𝑌)⁻¹𝑌_𝑠−1−𝑌. (28) Since 𝑌<𝑌_𝑠−1≤𝛽𝐼, then (𝐼−𝑌)⁻¹<(𝐼−𝑌_𝑠−1)⁻¹≤(1/(1−𝛽))𝐼; that is, (𝐼−𝑌)⁻¹<𝐼−𝑌_𝑠−1⁻¹≤1.1−𝛽 (29) Therefore, we get 𝑋_𝑠≤1−𝑋2𝛼^𝑘𝛽^2𝑘𝑌1−𝛽^𝑠−1−𝑌. (30) Also, we have 𝑌_𝑠−𝑌=^𝑛𝐴𝐼−𝑋_𝑠−1⁻¹𝐴^∗−^𝑛𝐴(I−𝑋)⁻¹𝐴^∗. (31) We denote 𝑃₂=𝐴𝐼−𝑋_𝑠−1⁻¹𝐴^∗,𝑄₂=𝐴(𝐼−𝑋)⁻¹𝐴^∗. (32) We use the following equality: ^𝑛𝑃₂^𝑛𝑃₂−^𝑛𝑄₂+^𝑛𝑃₂−^𝑛𝑄₂^𝑛𝑄₂=^𝑛/2P₂−^𝑛/2Q₂. (33)
Since 𝑋_𝑠>𝑋>0 for each 𝑠=0,1,2,…, then by using Lemma 1 we have the matrix 𝐷₂=^𝑛𝑃₂−^𝑛𝑄₂ being a positive definite solution of the matrix equation: ^𝑛𝑃₂𝐷₂+𝐷₂^𝑛𝑄₂=^𝑛/2𝑃₂−^𝑛/2𝑄₂. (34) According to Theorem 2, we have 𝐷₂=^∞₀𝑒^{−𝑛𝑃2𝑡𝑛/2}𝑃₂−^𝑛/2𝑄₂𝑒^{−𝑛𝑄2𝑡}𝑑𝑡. (35) Since ^𝑛𝑃₂,^𝑛𝑄₂ are positive definite matrices, then the integral (35) exists, and 𝑒^{−𝑛𝑃2𝑡𝑛/2}𝑃₂−^𝑛/2𝑄₂𝑒^{−𝑛𝑄2𝑡}→0,as𝑡→∞. (36) By using (31) and (35), we have 𝑌_𝑠−𝑌≤^∞₀^𝑛/2𝑃₂−^𝑛/2𝑄₂𝑒^{−𝑛𝑃2𝑡}𝑒^{−𝑛𝑄2𝑡}𝑑𝑡. (37) However, 𝑋_𝑠>𝑋≥𝛼𝐼; hence, ^𝑛𝑃₂>𝛼𝐼,^𝑛𝑄₂≥𝛼𝐼. (38) Then, we have 𝑌_𝑠−𝑌≤^𝑛/2𝑃₂−^𝑛/2𝑄₂^∞₀𝑒^{−𝑛𝑃2𝑡}𝑒^{−𝑛𝑄2𝑡}≤𝑑𝑡_𝑛/2𝑃₂−^𝑛/2𝑄₂^∞₀𝑒^{−𝛼𝐼𝑡2}=𝑑𝑡_𝑛/2𝑃₂−^𝑛/2𝑄₂^∞₀𝑒^−2𝛼𝑡=1𝑑𝑡2𝛼^𝑛/2𝑃₂−^𝑛/2𝑄₂≤12𝛼^22𝑛/2𝑃₂−^2𝑛/2𝑄₂. (39) After 𝑡 times as above, we get 𝑌_𝑠≤1−𝑌2𝛼^{𝑡𝑡𝑛/2}𝑃₂−^𝑡𝑛/2𝑄₂. (40)Let 𝑛=2^𝑡 be special case, then we have 𝑌_𝑠≤1−𝑌2𝛼^𝑡𝑃₂−𝑄₂=12𝛼^𝑡𝐴𝐼−𝑋_𝑠−1⁻¹𝐴^∗−𝐴(𝐼−𝑋)⁻¹𝐴^∗=12𝛼^𝑡𝐴𝐼−𝑋_𝑠−1⁻¹𝑋_𝑠−1−𝑋(𝐼−𝑋)⁻¹𝐴^∗≤12𝛼^𝑡‖𝐴‖²𝐼−𝑋_𝑠−1⁻¹(𝐼−𝑋)⁻¹𝑋_𝑠−1−𝑋. (41)
Since 𝑋<𝑋_𝑠−1≤𝛽𝐼, then (𝐼−𝑋)⁻¹<(𝐼−𝑋_𝑠−1)⁻¹≤(1/(1−𝛽))𝐼; that is, (𝐼−𝑋)⁻¹<𝐼−𝑋_𝑠−1⁻¹≤1.1−𝛽 (42) Therefore, we get 𝑌_𝑠≤1−𝑌2𝛼^𝑡𝛽^2𝑡𝑋1−𝛽^𝑠−1−𝑋. (43) By using (43) in (30) and (30) in (43), we have 𝑋_𝑠−𝑋≤𝑞𝑋_𝑠−2−𝑋≤⋯≤𝑞^𝑠/2,𝑌𝛽−𝛼_𝑠−𝑌≤𝑞𝑌_𝑠−2−𝑌≤⋯≤𝑞^𝑠/2𝛽−𝛼, (44) which completes the proof.

3. Numerical Examples

We will give some numerical examples for computing the positive definite solution of the Sys. (1). The solution is computed for some different matrices 𝐴,𝐵 with different orders. Denote by 𝑋,𝑌 the solutions which are obtained by Algorithm 3 and 𝜖₁(𝑋)=‖𝑋−𝑋_𝑠‖,𝜖₁(𝑌)=‖𝑌−𝑌_𝑠‖.

For computing 𝑍^1/𝑛 for all 1/𝑛∈(0,1], we use the iterative algorithm.

Algorithm 6. 𝑊₀𝑍=𝑍𝑊₀,𝑊_𝑠+1=1𝑛(𝑛−1)𝑊_𝑠+𝑊^1−𝑛_s𝑍,𝑠=0,1,2,…. (45)
See [30].

Example 7. Consider Sys. (1) with 𝑛=3,𝑚=2,𝛽=0.5, and normal matrices [],[].𝐴=diag0.2,0.2,0.5𝐵=diag0.4,0.3,0.1 (46) By using Algorithms 3 and 6, we have [],[].𝑋=diag0.535247,0.388225,0.174485𝑌=diag0.441516,0.402862,0.67154 (47)
The results are given in Table 1.

Table 1 

Error analysis for Example 7.

Example 8. Consider Sys. (1) with 𝑛=5,𝑚=4,𝛽=0.4 and normal matrices 𝐴=diag−72𝑟²,−62𝑟²,…,𝑟−82𝑟²,1𝐵=diag,28+𝑟𝑟8+𝑟,…,.8+𝑟 (48) By using Algorithms 3 and 6, we have;
When 𝑟=4, [],[].𝑋=diag0.362018,0.510342,0.620206,0.707524𝑌=diag0.595686,0.590503,0.577588,0.556604 (49)
When 𝑟=7, [][].𝑋=diag0.290121,0.408797,0.497892,0.570656,0.631828,0.68313,0.723835𝑌=diag0.372659,0.363432,0.349091,0.329437,0.302795,0.265305,0.206669 (50)(51)
The results are given in Table 2.

Table 2 

Error analysis for Example 8.

Example 9. Consider Sys. (1) with 𝑛=2,𝑚=2,𝛽=0.4, and matrices 𝐴=0.2−0.10.10.2,𝐵=−0.40.10.20.3. (52) By using Algorithms 3 and 6, we have ,,𝜆𝑋=0.497814−0.0839355−0.08393550.420533𝑌=0.333326−0.00556924−0.005569240.281592_𝑟(𝑋)={},𝜇0.551576,0.36677_𝑟(𝑌)={0.333919,0.280999}. (53)
The results are given in Table 3.

Table 3 

Error analysis for Example 9.

Example 10. Consider Sys. (1) with 𝑛=2,𝑚=1,𝛽=0.65, and matrices 𝐴=−0.2−0.4−0.4−0.6−0.20.60.2−0.40.1,𝐵=−0.03−0.080000.0200.10. (54) By using Algorithms 3 and 6, we have ,,𝜆𝑋=0.0573078−0.000371142−0.00702446−0.0003711420.000752950.0000206791−0.007024460.00002067910.000896794𝑌=0.596078−0.02336790.0748523−0.02336790.8870690.01045920.07485230.01045920.454549_𝑟(𝑋)={},𝜇0.0581718,0.000751403,0.0000343988_𝑟(𝑌)={0.888982,0.627173,0.42154}. (55)
The results are given in Table 4.

Table 4 

Error analysis for Example 10.

4. Conclusion

In this paper, the positive definite solutions for Sys. (1) have been tackled. We presented sufficient conditions for the existence of positive definite solutions for Sys. (1). Moreover, we discussed an iterative algorithm from which solutions can always be calculated numerically whenever the system is solvable. Finally, we gave numerical examples that illustrated the behavior of the proposed algorithm.