# Iterative Algorithm for Solving a System of Nonlinear Matrix Equations

Authors: Asmaa M. Al-Dubiban
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

We discuss the positive definite solutions for the system of nonlinear matrix equationsX−A∗Y−nA=I and Y−B∗X−mB=I, where n, m are two positive integers. Some properties of solutions are studied, and the necessary and sufficient conditions for the existence of positive definite solutions are given. An iterative algorithm for obtaining positive definite solutions of the system is proposed. Moreover, the error estimations are found. Finally, some numerical examples are given to show the efficiency of the proposed iterative algorithm.


---

## Body

## 1. Introduction

In this paper, we consider the system of nonlinear matrix equations that can be expressed in the form(1.1)X-A*Y-nA=I,Y-B*X-mB=I,
where n,m are two positive integers, X,Y are r×r unknown matrices, I  is the r×r identity matrix, and A,B are given nonsingular matrices. All matrices are defined over the complex field.System of nonlinear matrix equations of the form of (1.1) is a special case of the system of algebraic discrete-type Riccati equations of the form
(1.2)Xi=Vi*XiVi+Qi-(Vi*XiBi+Ai)(Ri+Bi*XiBi)-1(Bi*XiVi+Ai*),
where i=1,2,…,k [1, 2], when m=n=1. It is well known that the algebraic Riccati equations often arise in control theory, stability theory, communication system, dynamic programming, signal processing, statistics, and so forth [1–3]. In the recent years, some special case of the system (1.2) has been studied in many papers [4–14]. For example, Costa and Marques [5] have studied the maximal and stabilizing hermitian solutions for discrete-time-coupled algebraic Riccati equations. Czornik and Świerniak [7, 8] have studied the lower and the upper bounds on the solution of coupled algebraic Riccati equation. In [13] Mukaidani et al. proposed a numerical algorithm for finding solution of cross-coupled algebraic Riccati equations. In [4] Aldubiban has studied the properties of special case of Sys. (1.2) and obtained the sufficient conditions for existence of a positive definite solution and proposed an iterative algorithm to calculate the solutions. In [10] Ivanov proposed a method to solve the discrete-time-coupled algebraic Riccati equations.This paper is organized as following: in Section2, we derive the necessary and sufficient conditions of existence solutions for the Sys. (1.1). In Section 3, we introduce an iterative algorithm to obtain the positive definite solutions of Sys. (1.1). We discuss the convergence of the proposed iterative algorithm and study the convergence of the algorithm. Some numerical examples are given to illustrate the efficiency for suggested algorithm in Section 4.The following notations are used throughout the rest of the paper. The notationA≥0(A>0) means that A is a positive semidefinite (positive definite), A⋆ denotes the complex conjugate transpose of A, and I is the identity matrix. Moreover, A≥B(A>B) is used as a different notation for A-B≥0(A-B>0). We denote by ρ(A) the spectral radius of A, λr(X), and μr(Y) mean the eigenvalues of X and Y, respectively. The norm used in this paper is the spectral norm of the matrix A, that is, ∥A∥=ρ(AA⋆) unless otherwise noted.


## 2. Conditions for Existence of the Solutions

In this section, we will discuss some properties of the solutions for the Sys. (1.1), and we obtain the necessary and sufficient conditions for the existence of the solutions of Sys. (1.1).Theorem 2.1.
Ifλ-,λ+ are the smallest and the largest eigenvalues of a solution X of Sys. (1.1), respectively, and μ-,μ+  are the smallest and the largest eigenvalues of a solution Y of Sys. (1.1), respectively, η,ξ are eigenvalues of A,B, then
(2.1)(λ--1)μ-n≤|η|≤(λ+-1)μ+n,(μ--1)λ-m≤|ξ|≤(μ+-1)λ+m.Proof.
Letν be an eigenvector corresponding to an eigenvalue η of the matrix  A  and |ν|=1, and let  ω be   an eigenvector corresponding to an eigenvalue ξ of the matrix B  and |ω|=1. Since the solution (X,Y) of Sys. (1.1) is a positive definite solution then (λ--1)I≤X-I≤(λ+-1)I and (μ--1)I≤Y-I≤(μ+-1)I.
From the Sys. (1.1), we have
(2.2)〈Xν,ν〉-〈(A*Y-nA)ν,ν〉=〈Iν,ν〉,〈Xν,ν〉-〈Iν,ν〉=〈Y-nAν,Aν〉,〈(X-I)ν,ν〉=|η|2〈Y-nν,ν〉,
that is
(2.3)(λ--1)μ-n≤|η|2≤(λ+-1)μ+n.
Hence
(2.4)(λ--1)μ-n≤|η|≤(λ+-1)μ+n.
Also, from the Sys. (1.1), we have
(2.5)〈Yω,ω〉-〈(B*X-mB)ω,ω〉=〈Iω,ω〉,〈Yω,ω〉-〈Iω,ω〉=〈X-mBω,Bω〉,〈(Y-I)ω,ω〉=|ξ|2〈X-mω,ω〉,
that is
(2.6)(μ--1)λ-m≤|ξ|2≤(μ+-1)λ+m.
Hence
(2.7)(μ--1)λ-m≤|ξ|≤(μ+-1)λ+m.Theorem 2.2.
If Sys. (1.1) has a positive definite solution (X,Y), then
(2.8)I<X<I+A*A,I<Y<I+B*B.Proof.
Since(X,Y)  is a positive definite solution of Sys. (1.1), then
(2.9)X>I,A*Y-nA>0,Y>I,B*X-mB>0.
From the inequality Y-1<I, we have
(2.10)A*Y-nA<A*A,
that is
(2.11)X=I+A*Y-nA,<I+A*A.
Hence
(2.12)I<X<I+A*A.
From the inequality X-1<I, we have
(2.13)B*X-mB<B*B,
that is
(2.14)Y=I+B*X-mB,<I+B*B.
Hence
(2.15)I<Y<I+B*B.Theorem 2.3.
Sys. (1.1) has a positive definite solution (X,Y) if and only if the matrices A,B have the factorization
(2.16)A=(P*P)n/2N,B=(Q*Q)m/2M,
where P,Q are nonsingular matrices satisfying the following system:
(2.17)Q*Q-N*N=I,P*P-M*M=I.
In this case the solution is (Q*Q,P*P).Proof.
Let Sys. (1.1) have a positive definite solution (X,Y);  then X=Q*Q,Y=P*P, where Q,P are nonsingular matrices. Then Sys. (1.1) can be rewritten as
(2.18)Q*Q-A*(P*P)-nA=I,P*P-B*(Q*Q)-mB=I,Q*Q-A*(P*P)-n/2(P*P)-n/2A=I,P*P-B*(Q*Q)-m/2(Q*Q)-m/2B=I.
Letting N=(P*P)-n/2A, M=(Q*Q)-m/2B, then A=(P*P)n/2N, B=(Q*Q)m/2M, then the Sys. (1.1) is an equivalent to Sys. (2.17).
Conversely, ifA,B have the factorization (2.16) and satisfy Sys. (2.17), let X=Q*Q,Y=P*P,  then X,Y are positive definite matrices, and we have
(2.19)X-A*Y-nA=Q*Q-A*(P*P)-nA=Q*Q-A*(P*P)-n/2(P*P)-n/2A=Q*Q-N*N=I,(2.20)Y-B*X-mB=P*P-B*(Q*Q)-mB=P*P-B*(Q*Q)-m/2(Q*Q)-m/2B=P*P-M*M=I.
Hence Sys. (1.1) has a positive definite solution.


## 3. Iterative Algorithm for Solving the System

In this section, we will investigate the iterative solution of the Sys. (1.1). From this section to the end of the paper we will consider that the matrices A,B are normal satisyfing A-1B=BA-1 and A-1B*=B*A-1.Let us consider the following iterative algorithm.Algorithm 3.1.
TakeX0=I,Y0=I.
Fors=0,1,2,… compute
(3.1)Xs+1=I+A*Ys-nA,Ys+1=I+B*Xs-mB.Lemma 3.2.
For the Sys. (1.1), we have
(3.2)AXs=XsA,BYs=YsB,AYs-1=Ys-1A,BXs-1=Xs-1B,
where {Xs},{Ys},s=0,1,2,…, are determined by Algorithm 3.1.Proof.
SinceX0=Y0=I,  then
(3.3)AX0=X0A,BY0=Y0B,A-1Y0=Y0A-1,B-1X0=X0B-1.
Using the conditions AA*=A*A,BB*=B*B, we obtain
(3.4)AX1=A(I+A*A)=A+AA*A=A+A*AA=(I+A*A)A=X1A.
Also, we have
(3.5)BY1=Y1B.
Using the conditions A-1B=BA-1,A-1B*=B*A-1, we obtain
(3.6)A-1Y1=A-1(I+B*B)=A-1+A-1B*B=A-1+B*A-1B=A-1+B*BA-1=(I+B*B)A-1=Y1A-1.
By the same manner, we get
(3.7)B-1X1=X1B-1.
Further, assume that for each k it is satisfied  that
(3.8)AXk-1=Xk-1A,BYk-1=Yk-1B,A-1Yk-1=Yk-1A-1,B-1Xk-1=Xk-1B-1.
Now, by induction, we will prove
(3.9)AXk=XkA,BYk=YkB,A-1Yk=YkA-1,B-1Xk=XkB-1.
Since the two matrices A,B are normal and using the equalities (3.8), therefor
(3.10)AXk=A(I+A*Yk-1-nA)=A+AA*Yk-1-nA=A+A*AYk-1-nA=A+A*Yk-1-nAA=(I+A*Yk-1-nA)A=XkA.
Similarly
(3.11)BYk=YkB.
By using the conditions A-1B=BA-1,A-1B*=B*A-1 and the equalities (3.8), we have
(3.12)A-1Yk=A-1(I+B*Xk-1-mB)=A-1+A-1B*Xk-1-mB=A-1+B*A-1Xk-1-mB=A-1+B*Xk-1-mA-1B=A-1+B*Xk-1-mBA-1=(I+B*Xk-1-mB)A-1=YkA-1.
Also, we can prove
(3.13)B-1Xk=XkB-1.
Therefore, the equalities (3.2) are true for all s=0,1,2,….Corollary 3.3.
From Lemma3.2, we have
(3.14)AYs-n=Ys-nA,BXs-m=Xs-mB,
where {Xs},{Ys},s=0,1,2,…, are determined by Algorithm 3.1.Lemma 3.4.
For the Sys. (1.1), we have
(3.15)XsXs+1=Xs+1Xs,YsYs+1=Ys+1Ys,
where {Xs},{Ys},s=0,1,2,…, are determined by Algorithm 3.1.Proof.
SinceX0=Y0=I, then X0X1=X1X0,Y0Y1=Y1Y0.
By using the equalities (3.14), we have
(3.16)X1X2=(I+A*A)(I+A*Y1-nA)=I+A*A+A*Y1-nA+A*AA*Y1-nA=I+A*A+A*Y1-nA+A*AY1-nA*A=I+A*A+A*Y1-nA+A*Y1-nAA*A=(I+A*Y1-nA)(I+A*A)=X2X1.
Similarly we get
(3.17)Y1Y2=Y2Y1.
Further, assume that for each k it is satisfied that
(3.18)Xk-1Xk=XkXk-1,Yk-1Yk=YkYk-1.
Now, we will prove
(3.19)XkXk+1=Xk+1Xk,YkYk+1=Yk+1Yk.
From the equalities (3.18), we have
(3.20)Xk-1-mXk-m=Xk-mXk-1-m,Yk-1-nYk-n=Yk-nYk-1-n.
By using the equalities (3.14) and (3.20), we have
(3.21)XkXk+1=(I+A*Yk-1-nA)(I+A*Yk-nA)=I+A*Yk-1-nA+A*Yk-nA+A*Yk-1-nAA*Yk-nA=I+A*Yk-1-nA+A*Yk-nA+A*AYk-1-nYk-nA*A=I+A*Yk-1-nA+A*Yk-nA+A*AYk-nYk-1-nA*A=I+A*Yk-1-nA+A*Yk-nA+A*Yk-nAA*Yk-1-nA=(I+A*Yk-nA)(I+A*Yk-1-nA)=Xk+1Xk.
By the same manner,  we can prove
(3.22)YkYk+1=Yk+1Yk.
Therefore, the equalities (3.15) are true for all s=0,1,2,….Theorem 3.5.
IfA,B are satisfying the following conditions: (i)
∥A∥2(1+∥A∥2)m-1<1/m,(ii)
∥B∥2(1+∥B∥2)n-1<1/n, then the Sys. (1.1) has a positive definite solution (X,Y), which satisfy
(3.23)X2s<X<X2s+1,Y2s<Y<Y2s+1,s=0,1,2,…,max(∥X-X2s∥,∥X2s+1-X∥)<qs∥A∥2,max(∥Y-Y2s∥,∥Y2s+1-Y∥)<qs∥B∥2,
where q=nm∥A∥2∥B∥2(1+∥A∥2)m-1(1+∥B∥2)n-1<1, {Xs},{Ys},s=0,1,2,…, are determined by Algorithm 3.1.Proof.
ForX1,Y1 we have X1=I+A*A>X0  and  Y1=I+B*B>Y0.
SinceX1>X0,Y1>Y0 then X1-m<X0-m,Y1-n<Y0-n and B*X1-mB<B*X0-mB,A*Y1-nA<A*Y0-nA, hence I=X0<X2=I+A*Y1-nA<I+A*Y0-nA=X1, I=Y0<Y2=I+B*X1-mB<I+B*X0-mB=Y1, that is,
(3.24)X0<X2<X1,Y0<Y2<Y1.
We find the relation between X2,X3,X4,  and X5 and the relation between Y2,Y3,Y4, and Y5.
SinceX0<X2,Y0<Y2, then X3=I+A*Y2-nA<I+A*Y0-nA=X1, X3=I+A*Y2-nA>I+A*Y1-nA=X2, Y3=I+B*X2-mB<I+B*X0-mB=Y1, Y3=I+B*X2-mB>I+B*X1-mB=Y2.
SinceX2<X3<X1,Y2<Y3<Y1, then X4=I+A*Y3-nA>I+A*Y1-nA=X2, X4=I+A*Y3-nA<I+A*Y2-nA=X3, Y4=I+B*X3-mB>I+B*X1-mB=Y2, Y4=I+B*X3-mB<I+B*X2-mB=Y3.
Also sinceX2<X4<X3,Y2<Y4<Y3, then X5=I+A*Y4-nA<I+A*Y2-nA=X3, X5=I+A*Y4-nA>I+A*Y3-nA=X4, Y5=I+B*X4-mB<I+B*X2-mB=Y3, Y5=I+B*X4-mB>I+B*X3-mB=Y4.
Thus we get(3.25)X0<X2<X4<X5<X3<X1,Y0<Y2<Y4<Y5<Y3<Y1.
So, assume that for each k it is satisfied that
(3.26)X0<X2k<X2k+2<X2k+3<X2k+1<X1,Y0<Y2k<Y2k+2<Y2k+3<Y2k+1<Y1.
Now, we will prove X2k+2<X2k+4<X2k+5<X2k+3 and Y2k+2<Y2k+4<Y2k+5<Y2k+3.
By using the inequalities (3.26) we have
(3.27)X2k+4=I+A*Y2k+3-nA>I+A*Y2k+1-nA=X2k+2,X2k+4=I+A*Y2k+3-nA<I+A*Y2k+2-nA=X2k+3.
Also we have
(3.28)Y2k+4=I+B*X2k+3-mB>I+B*X2k+1-mB=Y2k+2,Y2k+4=I+B*X2k+3-mB<I+B*X2k+2-mB=Y2k+3.
Similarly
(3.29)X2k+5=I+A*Y2k+4-nA<I+A*Y2k+2-nA=X2k+3,X2k+5=I+A*Y2k+4-nA>I+A*Y2k+3-nA=X2k+4.
Also we have
(3.30)Y2k+5=I+B*X2k+4-mB<I+B*X2k+2-mB=Y2k+3,Y2k+5=I+B*X2k+4-mB>I+B*X2k+3-mB=Y2k+4.
Therefore, the inequalities (3.26) are true for all s=0,1,2,…; consequently the subsequences {X2s},{X2s+1},{Y2s}, and {Y2s+1} are monotonic and bounded. therefore they are convergent to positive definite matrices. To prove that  the sequences {X2s},{X2s+1} have a common limit, we have
(3.31)∥X2s+1-X2s∥=∥A*Y2s-nA-A*Y2s-1-nA∥=∥A*(Y2s-n-Y2s-1-n)A∥≤∥A∥2∥Y2s-n(Y2s-1n-Y2sn)Y2s-1-n∥≤∥A∥2∥Y2s-n∥∥Y2s-1-n∥∥Y2s-1n-Y2sn∥=∥A∥2∥Y2s-n∥∥Y2s-1-n∥∥(Y2s-1-Y2s)(∑i=1n‍Y2s-1n-iY2si-1)∥.
Since I<Ys<I+B*B, then we have
(3.32)∥Ys-n∥<1,∥Ys∥<1+∥B∥2.
Consequently
(3.33)∥X2s+1-X2s∥<∥A∥2∥Y2s-1-Y2s∥∑i=1n‍∥Y2s-1∥n-i∥Y2s∥i-1<n∥A∥2(1+∥B∥2)n-1∥Y2s-1-Y2s∥.
Also, to prove that the sequences {Y2s},{Y2s+1} have a common limit, we have
(3.34)∥Y2s+1-Y2s∥=∥B*X2s-mB-B*X2s-1-mB∥=∥B*(X2s-m-X2s-1-m)B∥≤∥B∥2∥X2s-m(X2s-1m-X2sm)X2s-1-m∥≤∥B∥2∥X2s-m∥∥X2s-1-m∥∥X2s-1m-X2sm∥=∥B∥2∥X2s-m∥∥X2s-1-m∥∥(X2s-1-X2s)(∑i=1m‍X2s-1m-iX2si-1)∥.
Since I<Xs<I+A*A, then we have
(3.35)∥Xs-m∥<1,∥Xs∥<1+∥A∥2.
Consequently
(3.36)∥Y2s+1-Y2s∥<∥B∥2∥X2s-1-X2s∥∑i=1m‍∥X2s-1∥m-i∥X2s∥i-1<m∥B∥2(1+∥A∥2)m-1∥X2s-1-X2s∥.
By using (3.36) in (3.33) and (3.33) in (3.36), we have
(3.37)∥X2s+1-X2s∥<nm∥A∥2∥B∥2(1+∥A∥2)m-1(1+∥B∥2)n-1∥X2s-1-X2s-2∥,∥Y2s+1-Y2s∥<nm∥A∥2∥B∥2(1+∥A∥2)m-1(1+∥B∥2)n-1∥Y2s-1-Y2s-2∥.
Therefore
(3.38)∥X2s+1-X2s∥<q∥X2s-1-X2s-2∥<⋯<qs∥X1-X0∥,∥Y2s+1-Y2s∥<q∥Y2s-1-Y2s-2∥<⋯<qs∥Y1-Y0∥.
Consequently the subsequences {X2s},{X2s+1} are convergent and have a common positive definite limit X. Also, the subsequences {Y2s},{Y2s+1} are convergent and have a common positive definite limit Y. Therefore (X,Y) is a positive definite solution of Sys. (1.1).


## 4. Numerical Examples

In this section the numerical examples are given to display the flexibility of the method. The solutions are computed for some different matricesA,B with different orders. In the following examples we denote by X,Y the solutions which are obtained by iterative Algorithm 3.1, and ϵ1(X)=∥X-Xs∥, ϵ2(X)=∥Xs-A*Ys-nA-I∥, ϵ1(Y)=∥Y-Ys∥,  and ϵ2(Y)=∥Ys-B*Xs-mB-I∥.Example 4.1.
Consider Sys. (1.1) with n=5, m=5 and matrices
(4.1)A=(0410-1801-2),B=(0-1-1-1-2-683-1).
By computation, we get
(4.2)X=(100019-5.999840-5.9998452.0729),Y=(1-2.49178×10-7-1.30736×10-6-2.49178×10-721-1.30736×10-612.00002),λr(X)={53.1277,17.9452,1},μr(Y)={3.00002,1.00001,1}.
The results are given in the Table 1.Table 1
s
ϵ
1
(
X
)
ϵ
1
(
Y
)
ϵ
2
(
X
)
ϵ
2
(
Y
)
0
51.0729
1.00002
69.0000
65.0000
1
17.9271
65.0000
59.7925
65.0000
2
41.8654
1.91479
E
-
06
41.8656
35.8160
3
2.71389
E
-
04
35.8160
42.8061
35.8160
4
42.8058
8.21136
E
-
11
42.8058
5.46612
E
-
02
5
1.14466
E
-
08
5.46612
E
-
02
6.06199
5.46612
E
-
02
6
6.06199
4.79780
E
-
15
6.06199
4.90522
E
-
05
7
7.60281
E
-
13
4.90522
E
-
05
4.26692
E
-
03
4.90522
E
-
05
8
4.26692
E
-
03
2.70627
E
-
19
4.26692
E
-
03
1.10251
E
-
09
9
0
1.10251
E
-
09
1.51071
E
-
07
1.10251
E
-
09
10
1.51071
E
-
07
0
1.51071
E
-
07
4.69266
E
-
14
11
0
4.69266
E
-
14
7.73070
E
-
12
4.69266
E
-
14Example 4.2.
Consider Sys. (1.1) with n=22, m=12 and matrices
(4.3)A=-0.1(0211240010421020),B=-0.1(1212200010012010).
By computation, we get
(4.4)X=(1.017170.02510930.0089920.004810470.02510931.11444-0.0619579-0.01276460.008992-0.06195791.124610.05232160.00481047-0.01276460.05232161.02878),Y=(1.0460.007308320.0167610.008766690.007308321.034040.01642780.03227750.0167610.01642781.016050.01284310.008766690.03227750.01284311.0343),λr(X)={1.19574,1.07941,1.00868,1.00118},μr(Y)={1.08177,1.04193,1.0064,1.0003}.
The results are given in the Table 2.Table 2
s
ϵ
1
(
X
)
ϵ
1
(
Y
)
ϵ
2
(
X
)
ϵ
2
(
Y
)
0
1.24608
E
-
01
4.60012
E
-
02
2.10000
E
-
01
1.00000
E
-
01
1
8.55563
E
-
02
5.39988
E
-
02
9.70059
E
-
02
7.95578
E
-
02
2
1.14496
E
-
02
2.55590
E
-
02
2.74686
E
-
02
3.15318
E
-
02
3
1.60190
E
-
02
0.59728
E
-
02
1.68545
E
-
02
1.36886
E
-
02
4
8.35485
E
-
04
7.71578
E
-
03
0.28393
E
-
02
8.14852
E
-
03
5
2.00382
E
-
03
4.32737
E
-
04
2.01636
E
-
03
1.42354
E
-
03
6
4.29946
E
-
05
9.90806
E
-
04
1.08242
E
-
04
9.89517
E
-
04
7
9.56981
E
-
05
0.19967
E
-
04
1.00788
E
-
04
4.57126
E
-
05
8
1.02495
E
-
05
4.08207
E
-
05
0.22025
E
-
04
4.13549
E
-
05
9
1.17755
E
-
05
4.68667
E
-
06
0.12947
E
-
04
9.61810
E
-
06
10
3.60760
E
-
06
4.93143
E
-
06
6.08098
E
-
06
5.52865
E
-
06
11
2.47338
E
-
06
1.44559
E
-
06
3.22091
E
-
06
2.24954
E
-
06
12
7.47536
E
-
07
1.06400
E
-
06
1.70353
E
-
06
1.15972
E
-
06
13
9.55996
E
-
07
3.31722
E
-
07
1.01939
E
-
06
7.40346
E
-
07
14
6.33981
E
-
08
4.08624
E
-
07
2.03521
E
-
07
4.37652
E
-
07
15
1.40123
E
-
07
2.90281
E
-
08
1.38180
E
-
07
9.25059
E
-
08
16
3.56253
E
-
09
6.34778
E
-
08
8.08867
E
-
09
6.24061
E
-
08


## 5. Conclusion

In this paper we considered the system of nonlinear matrix equations (1.1) where n,m are two positive integers. We achieved the general conditions for the existence of a positive definite solution. Moreover, we discussed an iterative algorithm from which a solution can always be calculated numerically, whenever the system is solvable. The numerical examples included in this paper showed the efficiency of the iterative algorithm which is described.



---

*Source: 461407-2012-12-04.xml*