# Iterative Solutions of a Set of Matrix Equations by Using the Hierarchical Identification Principle

Authors: Huamin Zhang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This paper is concerned with iterative solution to a class of the real coupled matrix equations. By using the
hierarchical identification principle, a gradient-based iterative algorithm is constructed to solve the real coupled
matrix equationsA
1
X
B
1
+
A
2
X
B
2
=
F
1 and C
1
X
D
1
+
C
2
X
D
2
=
F
2. The range of the convergence factor is derived to guarantee that the iterative algorithm is convergent for any initial value. The analysis indicates that
if the coupled matrix equations have a unique solution, then the iterative solution converges fast to the exact one
for any initial value under proper conditions. A numerical example is provided to illustrate the effectiveness of
the proposed algorithm.


---

## Body

## 1. Introduction

For systems with certain parameters, the controllability and stability are the important topics which are worth studying [1–3]. If the parameters of systems are uncertain, then how to identify these parameters is to put on the agenda. To identify the parameters of the large-scale systems, the hierarchical identification principle was proposed in [4–6]. The hierarchical gradient-based and least squares based identifications were presented for multivariable system [7]. The fruits of these effective strategies include identifications and adaptive control for dual-rate systems [8] and Hammerstein nonlinear systems [9].Many publications have studied the solutions to matrix equations from the different points of view [10–14]. Zhou et al. proposed the positive definite solutions of the nonlinear matrix equation X
+
A
H
X
¯
-
1
A
=
I [15]; Li et al. discussed a class of iterative methods for the generalized Sylvester equation [16]; the Riccati equation and a class of the coupled transpose matrix equations are investigated in [17, 18].Unlike the above methods and just like the Jacobi and Gauss-Seidel iterations, Ding and Chen proposed the gradient-based and least squares based iterations for solvingA
x
=
b and A
X
B
=
F [19, 20] and a large family of iterations. These iterations include the gradient iteration, the least squares iteration, and some classical iterations as their special cases [21, 22]. By using the hierarchical identification principle, the gradient-based iterative algorithms were derived for solving different real matrix equations, such as the generalized Sylvester matrix equations and general linear matrix equations [23–25].Ding’s strategy has received much attention. With the real representations of complex matrices as tools, Wu et al. applied the Ding’s strategy to solve the extended Sylvester-conjugate matrix equationsA
X
B
+
C
X
¯
D
=
F [26], the complex conjugate and transpose matrix equations [27], and the extended coupled Sylvester conjugate matrix equations [28]; Song and Chen presented the gradient based iterative algorithm for solving the extended Sylvester-conjugate transpose matrix equations [29].Different from the various single real matrix equations in [21–25] and the complex matrix equations in [26–29], the iterative algorithm for the real coupled matrix equations A
1
X
B
1
+
A
2
X
B
2
=
F
1 and C
1
X
D
1
+
C
2
X
D
2
=
F
2 is discussed by using the hierarchical identification principle, and the gradient-based iterative algorithm is proposed. Moreover, the gradient-based iteration is reported for solving the general real coupled matrix equations. We prove that the iterative solution always converges to the exact one for any initial value, if there exists a unique solution of these kinds of real coupled matrix equations.The iterative methods can be applied to nonlinear system identification [30, 31]. The proposed methods of this paper can combine the iterative identification methods [32–35], the auxiliary model identification methods [36–39], the multi-innovation identification methods [40–44], and the two-stage or multistage identification methods [45] to study identification problems for other linear systems [46–50] or nonlinear systems [51–54] and other systems with colored noises [55–58].This paper is organized as follows. Section2 offers some notation and basic lemmas. Section 3 derives the gradient-based iterative algorithm for solving the matrix equations A
1
X
B
1
+
A
2
X
B
2
=
F
1 and C
1
X
D
1
+
C
2
X
D
2
=
F
2. The iteration of real general coupled matrix equations is discussed in Section 4. Section 5 presents a numerical example to illustrate the effectiveness of the proposed algorithm. Finally, we offer some concluding remarks in Section 6.


## 2. Notations and Basic Lemmas

Some notation and lemmas are introduced first. The symbolA
T represents the transpose of A. For a square matrix A, λ
max
⁡
[
A
] indicates the maximum eigenvalue of A. ∥
A
∥ denotes the norm of A and is defined as ∥
A
∥
2
:
=
tr
⁡
[
A
T
A
]. For an m
×
n matrix
(1)
X
=
[
x
1
,
x
2
,
…
,
x
n
]
∈
ℝ
m
×
n
,
x
i
∈
ℝ
m
col
[
X
] is defined as
(2)
col
[
X
]
=
[
x
1
T
,
x
2
T
,
…
,
x
n
T
]
T
∈
ℝ
m
n
.
If A
=
[
a
i
j
]
∈
ℝ
m
×
n and B
∈
ℝ
p
×
q, then their Kronecker product is defined as
(3)
A
⊗
B
=
[
a
i
j
B
]
=
[
a
11
B
a
12
B
⋯
a
1
n
B
a
21
B
a
22
B
⋯
a
2
n
B
⋮
⋮
⋮
a
m
1
B
a
m
2
B
⋯
a
m
n
B
]
∈
ℝ
(
m
p
)
×
(
n
q
)
.The relationship between the vec-operatorcol and the Kronecker product can be expressed as the following lemma [59].Lemma 1.
IfA
∈
ℝ
m
×
n, B
∈
ℝ
p
×
q, and X
∈
ℝ
n
×
p, then we have
(4)
col
[
A
X
B
]
=
(
B
T
⊗
A
)
col
[
X
]
.The gradient-based iterative algorithm for solving matrix equationA
X
B
=
F is listed as follows [22].Lemma 2.
ForA
X
B
=
F, if A is a full-column rank matrix and B is a full-row rank matrix, then the iterative solution X
(
k
) given by the following gradient-based iterative algorithm converges to the exact solution X (i.e., lim
⁡
k
→
∞
X
(
k
)
=
X) for any initial values X
(
0
) [22]:
(5)
X
(
k
)
=
X
(
k
-
1
)
+
μ
A
T
[
F
-
A
X
(
k
-
1
)
B
]
B
T
,
0
<
μ
<
2
λ
max
⁡
[
A
A
T
]
λ
max
⁡
[
B
T
B
]
o
r
0
<
μ
⩽
2
∥
A
∥
2
∥
B
∥
2
.


## 3. The Coupled Matrix Equations

In this section, we consider the following real coupled matrix equations:(6)
A
1
X
B
1
+
A
2
X
B
2
=
F
1
,
C
1
X
D
1
+
C
2
X
D
2
=
F
2
,
where A
i
,
C
i
∈
ℝ
m
×
m, B
i
,
D
i
∈
ℝ
n
×
n, and F
i
∈
ℝ
m
×
n, i
=
1,2, are the given matrices, and X
∈
ℝ
m
×
n is the unknown matrix to be determined.


### 3.1. The Exact Solution

According to Lemma1, we rewrite (6) as
(7)
[
B
1
T
⊗
A
1
+
B
2
T
⊗
A
2
D
1
T
⊗
C
1
+
D
2
T
⊗
C
2
]
col
[
X
]
=
col
[
F
1
,
F
2
]
,
where [
F
1
,
F
2
] is a block matrix. The exact solution of (6) can be given by the following theorem.Theorem 3.
Equation (6) has a unique solution if and only if S
1 is a full-column rank matrix; in this case, the unique solution is
(8)
col
[
X
]
=
(
S
1
T
S
1
)
-
1
S
1
T
col
[
F
1
,
F
2
]
,
where
(9)
S
1
:
=
[
B
1
T
⊗
A
1
+
B
2
T
⊗
A
2
D
1
T
⊗
C
1
+
D
2
T
⊗
C
2
]
∈
ℝ
(
2
m
n
)
×
(
m
n
)
.
If F
1
=
F
2
=
0, then (6) has a unique solution X
=
0.


### 3.2. The Gradient Iterative Algorithm

The hierarchical identification principle implies that the related system can be decomposed into several subsystems, so we introduce four intermediate matrices:(10)
M
1
≔
F
1
-
A
2
X
B
2
,
M
2
≔
F
1
-
A
1
X
B
1
,
N
1
≔
F
2
-
C
2
X
D
2
,
N
2
≔
F
2
-
C
1
X
D
1
.
Using these four intermediate matrices, we decompose (6) into four subsystems
(11)
S
1
:
A
1
X
B
1
=
M
1
,
S
2
:
A
2
X
B
2
=
M
2
,
S
3
:
C
1
X
D
1
=
N
1
,
S
4
:
C
2
X
D
2
=
N
2
.
According to Lemma 2, it is not hard to get the iterative solutions X
1
(
k
), X
2
(
k
), X
3
(
k
), and X
4
(
k
) of these four subsystems S
1, S
2, S
3, and S
4 as follows:
(12)
X
1
(
k
)
=
X
1
(
k
-
1
)
+
μ
A
1
T
[
M
1
-
A
1
X
1
(
k
-
1
)
B
1
]
B
1
T
,
X
2
(
k
)
=
X
2
(
k
-
1
)
+
μ
A
2
T
[
M
2
-
A
2
X
2
(
k
-
1
)
B
2
]
B
2
T
,
X
3
(
k
)
=
X
3
(
k
-
1
)
+
μ
C
1
T
[
N
1
-
C
1
X
3
(
k
-
1
)
D
1
]
D
1
T
,
X
4
(
k
)
=
X
4
(
k
-
1
)
+
μ
C
2
T
[
N
2
-
C
2
X
4
(
k
-
1
)
D
2
]
D
2
T
.
We will determine the convergence factor μ later. Using (10), we have
(13)
X
1
(
k
)
=
X
1
(
k
-
1
)
+
μ
A
1
T
[
F
1
-
A
2
X
B
2
-
A
1
X
1
(
k
-
1
)
B
1
]
B
1
T
,
X
2
(
k
)
=
X
2
(
k
-
1
)
+
μ
A
2
T
[
F
1
-
A
1
X
B
1
-
A
2
X
2
(
k
-
1
)
B
2
]
B
2
T
,
X
3
(
k
)
=
X
3
(
k
-
1
)
+
μ
C
1
T
[
F
2
-
C
2
X
D
2
-
C
1
X
3
(
k
-
1
)
D
1
]
D
1
T
,
X
4
(
k
)
=
X
4
(
k
-
1
)
+
μ
C
2
T
[
F
2
-
C
1
X
D
1
-
C
2
X
4
(
k
-
1
)
D
2
]
D
2
T
.
Because the unknown matrix X appears in the right-hand side, this algorithm is impossible to realize. By using the hierarchical identification principle, we replace the unknown matrix X with its estimate at iteration k
-
1. Hence, we obtain
(14)
X
1
(
k
)
=
X
1
(
k
-
1
)
+
μ
A
1
T
×
[
F
1
-
A
2
X
1
(
k
-
1
)
B
2
-
A
1
X
1
(
k
-
1
)
B
1
]
B
1
T
,
X
2
(
k
)
=
X
2
(
k
-
1
)
+
μ
A
2
T
×
[
F
1
-
A
1
X
2
(
k
-
1
)
B
1
-
A
2
X
2
(
k
-
1
)
B
2
]
B
2
T
,
X
3
(
k
)
=
X
3
(
k
-
1
)
+
μ
C
1
T
×
[
F
2
-
C
2
X
3
(
k
-
1
)
D
2
-
C
1
X
3
(
k
-
1
)
D
1
]
D
1
T
,
X
4
(
k
)
=
X
4
(
k
-
1
)
+
μ
C
2
T
×
[
F
2
-
C
1
X
4
(
k
-
1
)
D
1
-
C
2
X
4
(
k
-
1
)
D
2
]
D
2
T
.
In fact, only an iterative solution X
(
k
) is needed in this algorithm; we take the average of X
1
(
k
), X
2
(
k
), X
3
(
k
), and X
4
(
k
) as X
(
k
) and obtain the gradient-based iterative algorithm:
(15)
X
(
k
)
=
X
1
(
k
)
+
X
2
(
k
)
+
X
3
(
k
)
+
X
4
(
k
)
4
,
(16)
X
1
(
k
)
=
X
(
k
-
1
)
+
μ
A
1
T
×
[
F
1
-
A
2
X
(
k
-
1
)
B
2
-
A
1
X
(
k
-
1
)
B
1
]
B
1
T
,
(17)
X
2
(
k
)
=
X
(
k
-
1
)
+
μ
A
2
T
×
[
F
1
-
A
1
X
(
k
-
1
)
B
1
-
A
2
X
(
k
-
1
)
B
2
]
B
2
T
,
(18)
X
3
(
k
)
=
X
(
k
-
1
)
+
μ
C
1
T
×
[
F
2
-
C
2
X
(
k
-
1
)
D
2
-
C
1
X
(
k
-
1
)
D
1
]
D
1
T
,
(19)
X
4
(
k
)
=
X
(
k
-
1
)
+
μ
C
2
T
×
[
F
2
-
C
1
X
(
k
-
1
)
D
1
-
C
2
X
(
k
-
1
)
D
2
]
D
2
T
,
(20)
0
<
μ
<
2
μ
0
,
μ
0
:
=
max
⁡
{
∑
i
=
1
2
‍
∥
A
i
∥
2
∥
B
i
∥
2
,
∑
i
=
1
2
‍
∥
C
i
∥
2
∥
D
i
∥
2
}
.Theorem 4.
If (6) has a unique solution X, then the iterative solution X
(
k
) given by (15)–(20) converges to X, that is, lim
⁡
k
→
∞
X
(
k
)
=
X, or the error matrix X
(
k
)
-
X converges to zero for any initial value X
(
0
).Proof.
We define estimation error matrices:(21)
X
~
i
(
k
)
:
=
X
i
(
k
)
-
X
,
i
=
1,2
,
3,4
,
X
~
(
k
)
:
=
X
(
k
)
-
X
.
From this definition and (15)–(19), we have the following error formulas:
(22)
X
~
1
(
k
)
=
X
~
(
k
-
1
)
-
μ
A
1
T
[
A
1
X
~
(
k
-
1
)
B
1
+
A
2
X
~
(
k
-
1
)
B
2
]
B
1
T
,
X
~
2
(
k
)
=
X
~
(
k
-
1
)
-
μ
A
2
T
[
A
1
X
~
(
k
-
1
)
B
1
+
A
2
X
~
(
k
-
1
)
B
2
]
B
2
T
,
X
~
3
(
k
)
=
X
~
(
k
-
1
)
-
μ
C
1
T
[
C
1
X
~
(
k
-
1
)
D
1
+
C
2
X
~
(
k
-
1
)
D
2
]
D
1
T
,
X
~
4
(
k
)
=
X
~
(
k
-
1
)
-
μ
C
2
T
[
C
1
X
~
(
k
-
1
)
D
1
+
C
2
X
~
(
k
-
1
)
D
2
]
D
2
T
,
(23)
X
~
(
k
)
=
[
X
~
1
(
k
)
+
X
~
2
(
k
)
+
X
~
3
(
k
)
+
X
~
4
(
k
)
]
4
.
Set
(24)
ξ
~
(
k
-
1
)
:
=
A
1
X
~
(
k
-
1
)
B
1
+
A
2
X
~
(
k
-
1
)
B
2
,
η
~
(
k
-
1
)
:
=
C
1
X
~
(
k
-
1
)
D
1
+
C
2
X
~
(
k
-
1
)
D
2
.
By the trace formula ∥
A
∥
2
=
tr
⁡
[
A
T
A
] and from (22), we obtain
(25)
∥
X
~
1
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
A
1
T
ξ
~
(
k
-
1
)
B
1
T
]
+
μ
2
∥
A
1
T
ξ
~
(
k
-
1
)
B
1
T
∥
2
,
(26)
∥
X
~
2
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
A
2
T
ξ
~
(
k
-
1
)
B
2
T
]
+
μ
2
∥
A
2
T
ξ
~
(
k
-
1
)
B
2
T
∥
2
,
(27)
∥
X
~
3
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
C
1
T
η
~
(
k
-
1
)
D
1
T
]
+
μ
2
∥
C
1
T
ξ
~
(
k
-
1
)
D
1
T
∥
2
,
(28)
∥
X
~
4
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
C
2
T
η
~
(
k
-
1
)
D
2
T
]
+
μ
2
∥
C
2
T
ξ
~
(
k
-
1
)
D
2
T
∥
2
.
Adding (25) to (26) gives
(29)
∥
X
~
1
(
k
)
∥
2
+
∥
X
~
2
(
k
)
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
A
1
T
ξ
~
(
k
-
1
)
B
1
T
]
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
A
2
T
ξ
~
(
k
-
1
)
B
2
T
]
+
μ
2
∥
A
1
T
ξ
~
(
k
-
1
)
B
1
T
∥
2
+
μ
2
∥
A
2
T
ξ
~
(
k
-
1
)
B
2
T
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
{
[
B
1
T
X
~
T
(
k
-
1
)
A
1
T
+
B
2
T
X
~
T
(
k
-
1
)
A
2
T
]
×
ξ
~
(
k
-
1
)
[
B
1
T
X
~
T
(
k
-
1
)
A
1
T
}
+
μ
2
∥
A
1
T
ξ
~
(
k
-
1
)
B
1
T
∥
2
+
μ
2
∥
A
2
T
ξ
~
(
k
-
1
)
B
2
T
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
{
[
A
1
X
~
(
k
-
1
)
B
1
+
A
2
X
~
(
k
-
1
)
B
2
]
T
×
ξ
~
(
k
-
1
)
}
+
μ
2
∥
A
1
T
ξ
~
(
k
-
1
)
B
1
T
∥
2
+
μ
2
∥
A
2
T
ξ
~
(
k
-
1
)
B
2
T
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
∥
ξ
~
(
k
-
1
)
∥
2
+
μ
2
∥
A
1
T
ξ
~
(
k
-
1
)
B
1
T
∥
2
+
μ
2
∥
A
2
T
ξ
~
(
k
-
1
)
B
2
T
∥
2
⩽
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
∥
ξ
~
(
k
-
1
)
∥
2
+
μ
2
(
∥
A
1
∥
2
∥
B
1
∥
2
+
∥
A
2
∥
2
∥
B
2
∥
2
)
×
∥
ξ
~
(
k
-
1
)
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
μ
[
2
-
μ
(
∥
A
1
∥
2
∥
B
1
∥
2
+
∥
A
2
∥
2
∥
B
2
∥
2
)
]
×
∥
ξ
~
(
k
-
1
)
∥
2
.
Similarly, adding (27) to (28) yields
(30)
∥
X
~
3
(
k
)
∥
2
+
∥
X
~
4
(
k
)
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
C
1
T
η
~
(
k
-
1
)
D
1
T
]
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
C
2
T
η
~
(
k
-
1
)
D
2
T
]
+
μ
2
∥
C
1
T
ξ
~
(
k
-
1
)
D
1
T
∥
2
+
μ
2
∥
C
2
T
ξ
~
(
k
-
1
)
D
2
T
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
∥
η
~
(
k
-
1
)
∥
2
+
μ
2
∥
C
1
T
ξ
~
(
k
-
1
)
D
1
T
∥
2
+
μ
2
∥
C
2
T
ξ
~
(
k
-
1
)
D
2
T
∥
2
⩽
2
∥
X
~
(
k
-
1
)
∥
2
-
μ
[
2
-
μ
(
∥
C
1
∥
2
∥
D
1
∥
2
+
∥
C
2
∥
2
∥
D
2
∥
2
)
]
×
∥
η
~
(
k
-
1
)
∥
2
.
Using (29) and (30), we have
(31)
∥
X
~
1
(
k
)
∥
2
+
∥
X
~
2
(
k
)
∥
2
+
∥
X
~
3
(
k
)
∥
2
+
∥
X
~
4
(
k
)
∥
2
⩽
4
∥
X
~
(
k
-
1
)
∥
2
-
μ
[
2
-
μ
(
∥
A
1
∥
2
∥
B
1
∥
2
+
∥
A
2
∥
2
∥
B
2
∥
2
)
]
∥
ξ
~
(
k
-
1
)
∥
2
-
μ
[
2
-
μ
(
∥
C
1
∥
2
∥
D
1
∥
2
+
∥
C
2
∥
2
∥
D
2
∥
2
)
]
∥
η
~
(
k
-
1
)
∥
2
.
Taking the norm of both sides of (23) gives
(32)
∥
X
~
(
k
)
∥
2
=
∥
X
~
1
(
k
)
+
X
~
2
(
k
)
+
X
~
3
(
k
)
+
X
~
4
(
k
)
∥
2
16
⩽
∥
X
~
1
(
k
)
∥
2
+
∥
X
~
2
(
k
)
∥
2
+
∥
X
~
3
(
k
)
∥
2
+
∥
X
~
4
(
k
)
∥
2
4
.
Substituting (31) into (32) yields
(33)
∥
X
~
(
k
)
∥
2
⩽
∥
X
~
(
k
-
1
)
∥
2
-
μ
4
[
2
-
μ
(
∥
A
1
∥
2
∥
B
1
∥
2
+
∥
A
2
∥
2
∥
B
2
∥
2
)
]
×
∥
ξ
~
(
k
-
1
)
∥
2
-
μ
4
[
2
-
μ
(
∥
C
1
∥
2
∥
D
1
∥
2
+
∥
C
2
∥
2
∥
D
2
∥
2
)
]
×
∥
η
~
(
k
-
1
)
∥
2
⩽
∥
X
~
(
0
)
∥
2
-
μ
4
[
2
-
μ
(
∥
A
1
∥
2
∥
B
1
∥
2
+
∥
A
2
∥
2
∥
B
2
∥
2
)
]
×
∑
i
=
0
k
-
1
∥
ξ
~
(
i
)
∥
2
-
μ
4
[
2
-
μ
(
∥
C
1
∥
2
∥
D
1
∥
2
+
∥
C
2
∥
2
∥
D
2
∥
2
)
]
×
∑
i
=
0
k
-
1
∥
η
~
(
i
)
∥
2
.
Thus, we have
(34)
μ
4
(
2
-
μ
∑
i
=
1
2
‍
∥
A
i
∥
2
∥
B
i
∥
2
)
∑
i
=
0
k
-
1
‍
∥
ξ
~
(
i
)
∥
2
+
μ
4
(
2
-
μ
∑
i
=
1
2
‍
∥
C
i
∥
2
∥
D
i
∥
2
)
∑
i
=
0
k
-
1
‍
∥
η
~
(
i
)
∥
2
⩽
∥
X
~
(
0
)
∥
2
<
∞
.
Since
(35)
μ
0
=
max
⁡
{
∑
i
=
1
2
‍
∥
A
i
∥
2
∥
B
i
∥
2
,
∑
i
=
1
2
‍
∥
C
i
∥
2
∥
D
i
∥
2
}
,
that is,
(36)
μ
4
(
2
-
μ
μ
0
)
[
∑
i
=
0
k
-
1
‍
(
∥
ξ
~
(
i
)
∥
2
+
∥
η
~
(
i
)
∥
2
)
]
⩽
∥
X
~
(
0
)
∥
2
,
if the convergence factor μ is chosen to satisfy 0
<
μ
<
2
/
μ
0, then we have
(37)
∑
k
=
1
∞
‍
(
∥
ξ
~
(
k
)
∥
2
+
∥
η
~
(
k
)
∥
2
)
⩽
∥
X
~
(
0
)
∥
2
.
It follows that ξ
~
(
k
)
→
0 and η
~
(
k
)
→
0, as i
→
∞ or as k
→
∞, we obtain
(38)
A
1
X
~
(
k
)
B
1
+
A
2
X
~
(
k
)
B
2
⟶
0
,
C
1
X
~
(
k
)
D
1
+
C
2
X
~
(
k
)
D
2
⟶
0
.
From Theorem 3, we have X
~
(
k
)
→
0, as k
→
∞. This completes the proof of Theorem 4.


## 3.1. The Exact Solution

According to Lemma1, we rewrite (6) as
(7)
[
B
1
T
⊗
A
1
+
B
2
T
⊗
A
2
D
1
T
⊗
C
1
+
D
2
T
⊗
C
2
]
col
[
X
]
=
col
[
F
1
,
F
2
]
,
where [
F
1
,
F
2
] is a block matrix. The exact solution of (6) can be given by the following theorem.Theorem 3.
Equation (6) has a unique solution if and only if S
1 is a full-column rank matrix; in this case, the unique solution is
(8)
col
[
X
]
=
(
S
1
T
S
1
)
-
1
S
1
T
col
[
F
1
,
F
2
]
,
where
(9)
S
1
:
=
[
B
1
T
⊗
A
1
+
B
2
T
⊗
A
2
D
1
T
⊗
C
1
+
D
2
T
⊗
C
2
]
∈
ℝ
(
2
m
n
)
×
(
m
n
)
.
If F
1
=
F
2
=
0, then (6) has a unique solution X
=
0.


## 3.2. The Gradient Iterative Algorithm

The hierarchical identification principle implies that the related system can be decomposed into several subsystems, so we introduce four intermediate matrices:(10)
M
1
≔
F
1
-
A
2
X
B
2
,
M
2
≔
F
1
-
A
1
X
B
1
,
N
1
≔
F
2
-
C
2
X
D
2
,
N
2
≔
F
2
-
C
1
X
D
1
.
Using these four intermediate matrices, we decompose (6) into four subsystems
(11)
S
1
:
A
1
X
B
1
=
M
1
,
S
2
:
A
2
X
B
2
=
M
2
,
S
3
:
C
1
X
D
1
=
N
1
,
S
4
:
C
2
X
D
2
=
N
2
.
According to Lemma 2, it is not hard to get the iterative solutions X
1
(
k
), X
2
(
k
), X
3
(
k
), and X
4
(
k
) of these four subsystems S
1, S
2, S
3, and S
4 as follows:
(12)
X
1
(
k
)
=
X
1
(
k
-
1
)
+
μ
A
1
T
[
M
1
-
A
1
X
1
(
k
-
1
)
B
1
]
B
1
T
,
X
2
(
k
)
=
X
2
(
k
-
1
)
+
μ
A
2
T
[
M
2
-
A
2
X
2
(
k
-
1
)
B
2
]
B
2
T
,
X
3
(
k
)
=
X
3
(
k
-
1
)
+
μ
C
1
T
[
N
1
-
C
1
X
3
(
k
-
1
)
D
1
]
D
1
T
,
X
4
(
k
)
=
X
4
(
k
-
1
)
+
μ
C
2
T
[
N
2
-
C
2
X
4
(
k
-
1
)
D
2
]
D
2
T
.
We will determine the convergence factor μ later. Using (10), we have
(13)
X
1
(
k
)
=
X
1
(
k
-
1
)
+
μ
A
1
T
[
F
1
-
A
2
X
B
2
-
A
1
X
1
(
k
-
1
)
B
1
]
B
1
T
,
X
2
(
k
)
=
X
2
(
k
-
1
)
+
μ
A
2
T
[
F
1
-
A
1
X
B
1
-
A
2
X
2
(
k
-
1
)
B
2
]
B
2
T
,
X
3
(
k
)
=
X
3
(
k
-
1
)
+
μ
C
1
T
[
F
2
-
C
2
X
D
2
-
C
1
X
3
(
k
-
1
)
D
1
]
D
1
T
,
X
4
(
k
)
=
X
4
(
k
-
1
)
+
μ
C
2
T
[
F
2
-
C
1
X
D
1
-
C
2
X
4
(
k
-
1
)
D
2
]
D
2
T
.
Because the unknown matrix X appears in the right-hand side, this algorithm is impossible to realize. By using the hierarchical identification principle, we replace the unknown matrix X with its estimate at iteration k
-
1. Hence, we obtain
(14)
X
1
(
k
)
=
X
1
(
k
-
1
)
+
μ
A
1
T
×
[
F
1
-
A
2
X
1
(
k
-
1
)
B
2
-
A
1
X
1
(
k
-
1
)
B
1
]
B
1
T
,
X
2
(
k
)
=
X
2
(
k
-
1
)
+
μ
A
2
T
×
[
F
1
-
A
1
X
2
(
k
-
1
)
B
1
-
A
2
X
2
(
k
-
1
)
B
2
]
B
2
T
,
X
3
(
k
)
=
X
3
(
k
-
1
)
+
μ
C
1
T
×
[
F
2
-
C
2
X
3
(
k
-
1
)
D
2
-
C
1
X
3
(
k
-
1
)
D
1
]
D
1
T
,
X
4
(
k
)
=
X
4
(
k
-
1
)
+
μ
C
2
T
×
[
F
2
-
C
1
X
4
(
k
-
1
)
D
1
-
C
2
X
4
(
k
-
1
)
D
2
]
D
2
T
.
In fact, only an iterative solution X
(
k
) is needed in this algorithm; we take the average of X
1
(
k
), X
2
(
k
), X
3
(
k
), and X
4
(
k
) as X
(
k
) and obtain the gradient-based iterative algorithm:
(15)
X
(
k
)
=
X
1
(
k
)
+
X
2
(
k
)
+
X
3
(
k
)
+
X
4
(
k
)
4
,
(16)
X
1
(
k
)
=
X
(
k
-
1
)
+
μ
A
1
T
×
[
F
1
-
A
2
X
(
k
-
1
)
B
2
-
A
1
X
(
k
-
1
)
B
1
]
B
1
T
,
(17)
X
2
(
k
)
=
X
(
k
-
1
)
+
μ
A
2
T
×
[
F
1
-
A
1
X
(
k
-
1
)
B
1
-
A
2
X
(
k
-
1
)
B
2
]
B
2
T
,
(18)
X
3
(
k
)
=
X
(
k
-
1
)
+
μ
C
1
T
×
[
F
2
-
C
2
X
(
k
-
1
)
D
2
-
C
1
X
(
k
-
1
)
D
1
]
D
1
T
,
(19)
X
4
(
k
)
=
X
(
k
-
1
)
+
μ
C
2
T
×
[
F
2
-
C
1
X
(
k
-
1
)
D
1
-
C
2
X
(
k
-
1
)
D
2
]
D
2
T
,
(20)
0
<
μ
<
2
μ
0
,
μ
0
:
=
max
⁡
{
∑
i
=
1
2
‍
∥
A
i
∥
2
∥
B
i
∥
2
,
∑
i
=
1
2
‍
∥
C
i
∥
2
∥
D
i
∥
2
}
.Theorem 4.
If (6) has a unique solution X, then the iterative solution X
(
k
) given by (15)–(20) converges to X, that is, lim
⁡
k
→
∞
X
(
k
)
=
X, or the error matrix X
(
k
)
-
X converges to zero for any initial value X
(
0
).Proof.
We define estimation error matrices:(21)
X
~
i
(
k
)
:
=
X
i
(
k
)
-
X
,
i
=
1,2
,
3,4
,
X
~
(
k
)
:
=
X
(
k
)
-
X
.
From this definition and (15)–(19), we have the following error formulas:
(22)
X
~
1
(
k
)
=
X
~
(
k
-
1
)
-
μ
A
1
T
[
A
1
X
~
(
k
-
1
)
B
1
+
A
2
X
~
(
k
-
1
)
B
2
]
B
1
T
,
X
~
2
(
k
)
=
X
~
(
k
-
1
)
-
μ
A
2
T
[
A
1
X
~
(
k
-
1
)
B
1
+
A
2
X
~
(
k
-
1
)
B
2
]
B
2
T
,
X
~
3
(
k
)
=
X
~
(
k
-
1
)
-
μ
C
1
T
[
C
1
X
~
(
k
-
1
)
D
1
+
C
2
X
~
(
k
-
1
)
D
2
]
D
1
T
,
X
~
4
(
k
)
=
X
~
(
k
-
1
)
-
μ
C
2
T
[
C
1
X
~
(
k
-
1
)
D
1
+
C
2
X
~
(
k
-
1
)
D
2
]
D
2
T
,
(23)
X
~
(
k
)
=
[
X
~
1
(
k
)
+
X
~
2
(
k
)
+
X
~
3
(
k
)
+
X
~
4
(
k
)
]
4
.
Set
(24)
ξ
~
(
k
-
1
)
:
=
A
1
X
~
(
k
-
1
)
B
1
+
A
2
X
~
(
k
-
1
)
B
2
,
η
~
(
k
-
1
)
:
=
C
1
X
~
(
k
-
1
)
D
1
+
C
2
X
~
(
k
-
1
)
D
2
.
By the trace formula ∥
A
∥
2
=
tr
⁡
[
A
T
A
] and from (22), we obtain
(25)
∥
X
~
1
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
A
1
T
ξ
~
(
k
-
1
)
B
1
T
]
+
μ
2
∥
A
1
T
ξ
~
(
k
-
1
)
B
1
T
∥
2
,
(26)
∥
X
~
2
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
A
2
T
ξ
~
(
k
-
1
)
B
2
T
]
+
μ
2
∥
A
2
T
ξ
~
(
k
-
1
)
B
2
T
∥
2
,
(27)
∥
X
~
3
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
C
1
T
η
~
(
k
-
1
)
D
1
T
]
+
μ
2
∥
C
1
T
ξ
~
(
k
-
1
)
D
1
T
∥
2
,
(28)
∥
X
~
4
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
C
2
T
η
~
(
k
-
1
)
D
2
T
]
+
μ
2
∥
C
2
T
ξ
~
(
k
-
1
)
D
2
T
∥
2
.
Adding (25) to (26) gives
(29)
∥
X
~
1
(
k
)
∥
2
+
∥
X
~
2
(
k
)
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
A
1
T
ξ
~
(
k
-
1
)
B
1
T
]
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
A
2
T
ξ
~
(
k
-
1
)
B
2
T
]
+
μ
2
∥
A
1
T
ξ
~
(
k
-
1
)
B
1
T
∥
2
+
μ
2
∥
A
2
T
ξ
~
(
k
-
1
)
B
2
T
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
{
[
B
1
T
X
~
T
(
k
-
1
)
A
1
T
+
B
2
T
X
~
T
(
k
-
1
)
A
2
T
]
×
ξ
~
(
k
-
1
)
[
B
1
T
X
~
T
(
k
-
1
)
A
1
T
}
+
μ
2
∥
A
1
T
ξ
~
(
k
-
1
)
B
1
T
∥
2
+
μ
2
∥
A
2
T
ξ
~
(
k
-
1
)
B
2
T
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
{
[
A
1
X
~
(
k
-
1
)
B
1
+
A
2
X
~
(
k
-
1
)
B
2
]
T
×
ξ
~
(
k
-
1
)
}
+
μ
2
∥
A
1
T
ξ
~
(
k
-
1
)
B
1
T
∥
2
+
μ
2
∥
A
2
T
ξ
~
(
k
-
1
)
B
2
T
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
∥
ξ
~
(
k
-
1
)
∥
2
+
μ
2
∥
A
1
T
ξ
~
(
k
-
1
)
B
1
T
∥
2
+
μ
2
∥
A
2
T
ξ
~
(
k
-
1
)
B
2
T
∥
2
⩽
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
∥
ξ
~
(
k
-
1
)
∥
2
+
μ
2
(
∥
A
1
∥
2
∥
B
1
∥
2
+
∥
A
2
∥
2
∥
B
2
∥
2
)
×
∥
ξ
~
(
k
-
1
)
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
μ
[
2
-
μ
(
∥
A
1
∥
2
∥
B
1
∥
2
+
∥
A
2
∥
2
∥
B
2
∥
2
)
]
×
∥
ξ
~
(
k
-
1
)
∥
2
.
Similarly, adding (27) to (28) yields
(30)
∥
X
~
3
(
k
)
∥
2
+
∥
X
~
4
(
k
)
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
C
1
T
η
~
(
k
-
1
)
D
1
T
]
-
2
μ
tr
⁡
[
X
~
T
(
k
-
1
)
C
2
T
η
~
(
k
-
1
)
D
2
T
]
+
μ
2
∥
C
1
T
ξ
~
(
k
-
1
)
D
1
T
∥
2
+
μ
2
∥
C
2
T
ξ
~
(
k
-
1
)
D
2
T
∥
2
=
2
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
∥
η
~
(
k
-
1
)
∥
2
+
μ
2
∥
C
1
T
ξ
~
(
k
-
1
)
D
1
T
∥
2
+
μ
2
∥
C
2
T
ξ
~
(
k
-
1
)
D
2
T
∥
2
⩽
2
∥
X
~
(
k
-
1
)
∥
2
-
μ
[
2
-
μ
(
∥
C
1
∥
2
∥
D
1
∥
2
+
∥
C
2
∥
2
∥
D
2
∥
2
)
]
×
∥
η
~
(
k
-
1
)
∥
2
.
Using (29) and (30), we have
(31)
∥
X
~
1
(
k
)
∥
2
+
∥
X
~
2
(
k
)
∥
2
+
∥
X
~
3
(
k
)
∥
2
+
∥
X
~
4
(
k
)
∥
2
⩽
4
∥
X
~
(
k
-
1
)
∥
2
-
μ
[
2
-
μ
(
∥
A
1
∥
2
∥
B
1
∥
2
+
∥
A
2
∥
2
∥
B
2
∥
2
)
]
∥
ξ
~
(
k
-
1
)
∥
2
-
μ
[
2
-
μ
(
∥
C
1
∥
2
∥
D
1
∥
2
+
∥
C
2
∥
2
∥
D
2
∥
2
)
]
∥
η
~
(
k
-
1
)
∥
2
.
Taking the norm of both sides of (23) gives
(32)
∥
X
~
(
k
)
∥
2
=
∥
X
~
1
(
k
)
+
X
~
2
(
k
)
+
X
~
3
(
k
)
+
X
~
4
(
k
)
∥
2
16
⩽
∥
X
~
1
(
k
)
∥
2
+
∥
X
~
2
(
k
)
∥
2
+
∥
X
~
3
(
k
)
∥
2
+
∥
X
~
4
(
k
)
∥
2
4
.
Substituting (31) into (32) yields
(33)
∥
X
~
(
k
)
∥
2
⩽
∥
X
~
(
k
-
1
)
∥
2
-
μ
4
[
2
-
μ
(
∥
A
1
∥
2
∥
B
1
∥
2
+
∥
A
2
∥
2
∥
B
2
∥
2
)
]
×
∥
ξ
~
(
k
-
1
)
∥
2
-
μ
4
[
2
-
μ
(
∥
C
1
∥
2
∥
D
1
∥
2
+
∥
C
2
∥
2
∥
D
2
∥
2
)
]
×
∥
η
~
(
k
-
1
)
∥
2
⩽
∥
X
~
(
0
)
∥
2
-
μ
4
[
2
-
μ
(
∥
A
1
∥
2
∥
B
1
∥
2
+
∥
A
2
∥
2
∥
B
2
∥
2
)
]
×
∑
i
=
0
k
-
1
∥
ξ
~
(
i
)
∥
2
-
μ
4
[
2
-
μ
(
∥
C
1
∥
2
∥
D
1
∥
2
+
∥
C
2
∥
2
∥
D
2
∥
2
)
]
×
∑
i
=
0
k
-
1
∥
η
~
(
i
)
∥
2
.
Thus, we have
(34)
μ
4
(
2
-
μ
∑
i
=
1
2
‍
∥
A
i
∥
2
∥
B
i
∥
2
)
∑
i
=
0
k
-
1
‍
∥
ξ
~
(
i
)
∥
2
+
μ
4
(
2
-
μ
∑
i
=
1
2
‍
∥
C
i
∥
2
∥
D
i
∥
2
)
∑
i
=
0
k
-
1
‍
∥
η
~
(
i
)
∥
2
⩽
∥
X
~
(
0
)
∥
2
<
∞
.
Since
(35)
μ
0
=
max
⁡
{
∑
i
=
1
2
‍
∥
A
i
∥
2
∥
B
i
∥
2
,
∑
i
=
1
2
‍
∥
C
i
∥
2
∥
D
i
∥
2
}
,
that is,
(36)
μ
4
(
2
-
μ
μ
0
)
[
∑
i
=
0
k
-
1
‍
(
∥
ξ
~
(
i
)
∥
2
+
∥
η
~
(
i
)
∥
2
)
]
⩽
∥
X
~
(
0
)
∥
2
,
if the convergence factor μ is chosen to satisfy 0
<
μ
<
2
/
μ
0, then we have
(37)
∑
k
=
1
∞
‍
(
∥
ξ
~
(
k
)
∥
2
+
∥
η
~
(
k
)
∥
2
)
⩽
∥
X
~
(
0
)
∥
2
.
It follows that ξ
~
(
k
)
→
0 and η
~
(
k
)
→
0, as i
→
∞ or as k
→
∞, we obtain
(38)
A
1
X
~
(
k
)
B
1
+
A
2
X
~
(
k
)
B
2
⟶
0
,
C
1
X
~
(
k
)
D
1
+
C
2
X
~
(
k
)
D
2
⟶
0
.
From Theorem 3, we have X
~
(
k
)
→
0, as k
→
∞. This completes the proof of Theorem 4.


## 4. General Real Coupled Matrix Equations

In this section, we study the general real coupled matrix equations of form(39)
A
1
X
B
1
+
A
2
X
B
2
+
⋯
+
A
p
X
B
p
=
F
1
,
C
1
X
D
1
+
C
2
X
D
2
+
⋯
+
C
q
X
D
q
=
F
2
,
where A
i
,
C
j
∈
ℝ
m
×
m, B
i
,
D
j
∈
ℝ
n
×
n, i
=
1,2
,
…
,
p, j
=
1,2
,
…
,
q, and F
1
,
F
2
∈
ℝ
m
×
n are the given matrices, and X
∈
ℝ
m
×
n is the unknown matrix to be determined.


### 4.1. The Exact Solution

According to Lemma1, (39) can be rewritten as
(40)
[
∑
i
=
1
p
‍
B
i
T
⊗
A
i
∑
i
=
1
q
‍
D
i
T
⊗
C
i
]
col
[
X
]
=
col
[
F
1
,
F
2
]
,
where [
F
1
,
F
2
] is a block matrix. Set
(41)
S
2
:
=
[
∑
i
=
1
p
‍
B
i
T
⊗
A
i
∑
i
=
1
q
‍
D
i
T
⊗
C
i
]
∈
ℝ
(
2
m
n
)
×
(
m
n
)
.
The exact solution of (39) can be given by the following theorem.Theorem 5.
Equation (39) has a unique solution if and only if S
2 is a full-column rank matrix; in this case, the unique solution is
(42)
col
[
X
]
=
(
S
2
T
S
2
)
-
1
S
2
T
col
[
F
1
,
F
2
]
.
IfF
1
=
F
2
=
0, then (39) has a unique solution X
=
0.


### 4.2. The Gradient-Based Iterative Algorithm

We define the intermediate matrices:(43)
G
i
:
=
F
1
-
∑
r
=
1
,
r
≠
i
p
A
r
X
B
r
,
i
=
1,2
,
…
,
p
,
H
j
:
=
F
2
-
∑
s
=
1
,
s
≠
j
q
C
s
X
D
s
,
j
=
1,2
,
…
,
q
.
By using the hierarchical identification principle, we decompose (39) into
(44)
G
i
=
A
i
X
B
i
,
i
=
1,2
,
…
,
p
,
H
j
=
C
j
X
D
j
,
j
=
1,2
,
…
,
q
.
According to Lemma 2, we have
(45)
Y
i
(
k
)
=
Y
i
(
k
-
1
)
+
μ
A
i
T
[
G
i
-
A
i
Y
i
(
k
-
1
)
B
i
]
B
i
T
,
(46)
Z
j
(
k
)
=
Z
j
(
k
-
1
)
+
μ
C
j
T
[
H
j
-
C
j
Z
j
(
k
-
1
)
D
j
]
D
j
T
,
where Y
i
(
k
) and Z
j
(
k
) denote the iterative solutions of equations G
i
=
A
i
X
B
i and H
j
=
C
j
X
D
j at iteration k, respectively. Substituting G
i into (45) and H
j into (46) gives
(47)
Y
i
(
k
)
=
Y
i
(
k
-
1
)
+
μ
A
i
T
[
F
1
-
∑
r
=
1
,
r
≠
i
p
A
r
X
B
r
-
A
i
Y
i
(
k
-
1
)
B
i
]
B
i
T
,
Z
j
(
k
)
=
Z
j
(
k
-
1
)
+
μ
C
j
T
[
F
2
-
∑
s
=
1
,
s
≠
j
q
C
s
X
D
s
-
C
j
Z
j
(
k
-
1
)
D
j
]
D
j
T
.
To realize the above algorithm, we replace the unknown matrix X in the right-hand side of the above two equations with Y
i
(
k
-
1
) and Z
j
(
k
-
1
), respectively, and obtain
(48)
Y
i
(
k
)
=
Y
i
(
k
-
1
)
+
μ
A
i
T
[
F
1
-
∑
r
=
1
p
‍
A
r
Y
i
(
k
-
1
)
B
r
]
B
i
T
,
0000000000000000000000000000
i
=
1,2
,
…
,
p
,
Z
j
(
k
)
=
Z
j
(
k
-
1
)
+
μ
C
j
T
[
F
2
-
∑
s
=
1
q
‍
C
s
Z
j
(
k
-
1
)
D
s
]
D
j
T
,
000000000000000000000000000000
j
=
1,2
,
…
,
q
.
In fact, only an iterative solution X
(
k
) is needed; taking the average of Y
i
(
k
) and Z
j
(
k
) as X
(
k
), we obtain the gradient iterative algorithm for solving (39) as follows:
(49)
X
(
k
)
=
∑
i
=
1
p
Y
i
(
k
)
+
∑
j
=
1
q
Z
j
(
k
)
p
+
q
,
(50)
Y
i
(
k
)
=
X
(
k
-
1
)
+
μ
A
i
T
[
F
1
-
∑
r
=
1
p
‍
A
r
X
(
k
-
1
)
B
r
]
B
i
T
,
(51)
Z
j
(
k
)
=
X
(
k
-
1
)
+
μ
C
j
T
[
F
2
-
∑
s
=
1
q
‍
C
s
X
(
k
-
1
)
D
s
]
D
j
T
,
(52)
0
<
μ
<
2
μ
0
,
μ
0
:
=
max
⁡
⁡
{
∑
i
=
1
p
‍
∥
A
i
∥
2
∥
B
i
∥
2
,
∑
j
=
1
q
‍
∥
C
j
∥
2
∥
D
j
∥
2
}
.Theorem 6.
If (39) has a unique solution X, then the iterative solution X
(
k
) given by (49)–(52) converges to X, that is, lim
⁡
k
→
∞
X
(
k
)
=
X, or the error matrix X
(
k
)
-
X converges to zero for any initial value X
(
0
).Proof.
The estimation error matrices are defined as(53)
Y
~
i
(
k
)
≔
Y
i
(
k
)
-
X
,
i
=
1,2
,
…
,
p
,
Z
~
j
(
k
)
≔
Z
j
(
k
)
-
X
,
j
=
1,2
,
…
,
q
,
X
~
(
k
)
≔
X
(
k
)
-
X
.
By using (49)–(51), we have
(54)
Y
~
i
(
k
)
=
X
~
(
k
-
1
)
-
μ
A
i
T
[
∑
r
=
1
p
‍
A
r
X
~
(
k
-
1
)
B
r
]
B
i
T
,
000000000000000000000000000000
i
=
1,2
,
…
,
p
,
Z
~
j
(
k
)
=
X
~
(
k
-
1
)
-
μ
C
j
T
[
∑
s
=
1
q
‍
C
s
X
~
(
k
-
1
)
D
s
]
D
j
T
,
000000000000000000000000000000
j
=
1,2
,
…
,
q
,
(55)
X
~
(
k
)
=
∑
i
=
1
p
Y
~
i
(
k
)
+
∑
j
=
1
q
Z
~
i
(
k
)
p
+
q
.
Set
(56)
ξ
~
(
k
-
1
)
:
=
∑
r
=
1
p
‍
A
r
X
~
(
k
-
1
)
B
r
,
η
~
(
k
-
1
)
:
=
∑
s
=
1
q
‍
C
s
X
~
(
k
-
1
)
D
s
.
According to the trace formula ∥
X
∥
2
=
tr
⁡
[
X
T
X
] and using (54), for i
=
1,2
,
…
,
p and j
=
1,2
,
…
,
q, we obtain
(57)
∥
Y
~
i
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
(
k
-
1
)
T
A
i
T
ξ
~
(
k
-
1
)
B
i
T
]
+
μ
2
∥
A
i
T
ξ
~
(
k
-
1
)
B
i
T
∥
2
,
(58)
∥
Z
~
j
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
(
k
-
1
)
T
C
j
T
η
~
(
k
-
1
)
D
j
T
]
+
μ
2
∥
C
j
T
η
~
(
k
-
1
)
D
j
T
∥
2
.
For i
=
1,2
,
…
,
p and j
=
1,2
,
…
,
q, adding (57) to (58) gives
(59)
∑
i
=
1
p
‍
∥
Y
~
i
(
k
)
∥
2
+
∑
j
=
1
q
∥
Z
~
i
(
k
)
∥
2
=
(
p
+
q
)
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
(
∥
ξ
~
(
k
-
1
)
∥
2
+
∥
η
~
(
k
-
1
)
∥
2
)
+
μ
2
∑
i
=
1
p
‍
∥
A
i
T
ξ
~
(
k
-
1
)
B
i
T
∥
2
+
μ
2
∑
j
=
1
q
‍
∥
C
j
T
η
~
(
k
-
1
)
D
j
T
∥
2
⩽
(
p
+
q
)
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
(
∥
ξ
~
(
k
-
1
)
∥
2
+
∥
η
~
(
k
-
1
)
∥
2
)
+
μ
2
∥
ξ
~
(
k
-
1
)
∥
2
∑
i
=
1
p
‍
∥
A
i
∥
2
∥
B
i
∥
2
+
μ
2
∥
η
~
(
k
-
1
)
∥
2
∑
j
=
1
q
‍
∥
C
j
∥
2
∥
D
j
∥
2
=
(
p
+
q
)
∥
X
~
(
k
-
1
)
∥
2
-
μ
(
2
-
μ
∑
i
=
1
p
‍
∥
A
i
∥
2
∥
B
i
∥
2
)
∥
ξ
~
(
k
-
1
)
∥
2
-
μ
(
2
-
μ
∑
j
=
1
q
‍
∥
C
j
∥
2
∥
D
j
∥
2
)
∥
η
~
(
k
-
1
)
∥
2
.
Taking the norm of both sides of (55) and using the norm inequality ∥
∑
i
=
1
m
‍
A
i
∥
2
⩽
m
∑
i
=
1
m
‍
∥
A
i
∥
2 gives
(60)
∥
X
~
(
k
)
∥
2
=
∥
∑
i
=
1
p
Y
~
i
(
k
)
+
∑
j
=
1
q
Z
~
j
(
k
)
∥
2
(
p
+
q
)
2
⩽
∑
i
=
1
p
∥
Y
~
i
(
k
)
∥
2
+
∑
j
=
1
q
∥
Z
~
j
(
k
)
∥
2
p
+
q
.
Dividing (59) by p
+
q and using the inequality (60) gives
(61)
∥
X
~
(
k
)
∥
2
⩽
∥
X
~
(
k
-
1
)
∥
2
-
μ
p
+
q
(
2
-
μ
∑
i
=
1
p
‍
∥
A
i
∥
2
∥
B
i
∥
2
)
∥
ξ
~
(
k
-
1
)
∥
2
-
μ
p
+
q
(
2
-
μ
∑
j
=
1
q
‍
∥
C
j
∥
2
∥
D
j
∥
2
)
∥
η
~
(
k
-
1
)
∥
2
⩽
∥
X
~
(
k
-
1
)
∥
2
-
μ
p
+
q
(
2
-
μ
μ
0
)
×
(
∥
ξ
~
(
k
-
1
)
∥
2
+
∥
η
~
(
k
-
1
)
∥
2
)
⩽
∥
X
~
(
0
)
∥
2
-
μ
p
+
q
(
2
-
μ
μ
0
)
∑
t
=
0
k
-
1
(
∥
ξ
~
(
t
)
∥
2
+
∥
η
~
(
t
)
∥
2
)
.
In the second “⩽,” the assumption that μ
0
=
max
⁡
{
∑
i
=
1
p
‍
∥
A
i
∥
2
∥
B
i
∥
2
,
∑
j
=
1
q
‍
∥
C
j
∥
2
∥
D
j
∥
2
} is used. Since ∥
X
~
(
k
)
∥
2
⩾
0, we have
(62)
μ
p
+
q
(
2
-
μ
μ
0
)
∑
t
=
0
k
-
1
‍
(
∥
ξ
~
(
t
)
∥
2
+
∥
η
~
(
t
)
∥
2
)
⩽
∥
X
~
(
0
)
∥
2
<
∞
.
If the convergence factor μ is chosen to satisfy 0
<
μ
<
2
/
μ
0, then we have
(63)
∑
k
=
0
∞
‍
(
∥
ξ
~
(
k
)
∥
2
+
∥
η
~
(
k
)
∥
2
)
⩽
∥
X
~
(
0
)
∥
2
.
It follows that ∥
η
~
(
k
)
∥
2
+
∥
η
~
(
k
)
∥
2
→
0, as k
→
∞ or as k
→
∞, we obtain
(64)
A
1
X
~
(
k
-
1
)
B
1
+
A
2
X
~
(
k
-
1
)
B
2
+
⋯
+
A
p
X
~
(
k
-
1
)
B
p
⟶
0
,
C
1
X
~
(
k
-
1
)
D
1
+
C
2
X
~
(
k
-
1
)
D
2
+
⋯
+
C
q
X
~
(
k
-
1
)
D
q
⟶
0
.
According to Theorem 5, we have X
~
(
k
)
→
0, as k
→
∞. This completes the proof of Theorem 6.


## 4.1. The Exact Solution

According to Lemma1, (39) can be rewritten as
(40)
[
∑
i
=
1
p
‍
B
i
T
⊗
A
i
∑
i
=
1
q
‍
D
i
T
⊗
C
i
]
col
[
X
]
=
col
[
F
1
,
F
2
]
,
where [
F
1
,
F
2
] is a block matrix. Set
(41)
S
2
:
=
[
∑
i
=
1
p
‍
B
i
T
⊗
A
i
∑
i
=
1
q
‍
D
i
T
⊗
C
i
]
∈
ℝ
(
2
m
n
)
×
(
m
n
)
.
The exact solution of (39) can be given by the following theorem.Theorem 5.
Equation (39) has a unique solution if and only if S
2 is a full-column rank matrix; in this case, the unique solution is
(42)
col
[
X
]
=
(
S
2
T
S
2
)
-
1
S
2
T
col
[
F
1
,
F
2
]
.
IfF
1
=
F
2
=
0, then (39) has a unique solution X
=
0.


## 4.2. The Gradient-Based Iterative Algorithm

We define the intermediate matrices:(43)
G
i
:
=
F
1
-
∑
r
=
1
,
r
≠
i
p
A
r
X
B
r
,
i
=
1,2
,
…
,
p
,
H
j
:
=
F
2
-
∑
s
=
1
,
s
≠
j
q
C
s
X
D
s
,
j
=
1,2
,
…
,
q
.
By using the hierarchical identification principle, we decompose (39) into
(44)
G
i
=
A
i
X
B
i
,
i
=
1,2
,
…
,
p
,
H
j
=
C
j
X
D
j
,
j
=
1,2
,
…
,
q
.
According to Lemma 2, we have
(45)
Y
i
(
k
)
=
Y
i
(
k
-
1
)
+
μ
A
i
T
[
G
i
-
A
i
Y
i
(
k
-
1
)
B
i
]
B
i
T
,
(46)
Z
j
(
k
)
=
Z
j
(
k
-
1
)
+
μ
C
j
T
[
H
j
-
C
j
Z
j
(
k
-
1
)
D
j
]
D
j
T
,
where Y
i
(
k
) and Z
j
(
k
) denote the iterative solutions of equations G
i
=
A
i
X
B
i and H
j
=
C
j
X
D
j at iteration k, respectively. Substituting G
i into (45) and H
j into (46) gives
(47)
Y
i
(
k
)
=
Y
i
(
k
-
1
)
+
μ
A
i
T
[
F
1
-
∑
r
=
1
,
r
≠
i
p
A
r
X
B
r
-
A
i
Y
i
(
k
-
1
)
B
i
]
B
i
T
,
Z
j
(
k
)
=
Z
j
(
k
-
1
)
+
μ
C
j
T
[
F
2
-
∑
s
=
1
,
s
≠
j
q
C
s
X
D
s
-
C
j
Z
j
(
k
-
1
)
D
j
]
D
j
T
.
To realize the above algorithm, we replace the unknown matrix X in the right-hand side of the above two equations with Y
i
(
k
-
1
) and Z
j
(
k
-
1
), respectively, and obtain
(48)
Y
i
(
k
)
=
Y
i
(
k
-
1
)
+
μ
A
i
T
[
F
1
-
∑
r
=
1
p
‍
A
r
Y
i
(
k
-
1
)
B
r
]
B
i
T
,
0000000000000000000000000000
i
=
1,2
,
…
,
p
,
Z
j
(
k
)
=
Z
j
(
k
-
1
)
+
μ
C
j
T
[
F
2
-
∑
s
=
1
q
‍
C
s
Z
j
(
k
-
1
)
D
s
]
D
j
T
,
000000000000000000000000000000
j
=
1,2
,
…
,
q
.
In fact, only an iterative solution X
(
k
) is needed; taking the average of Y
i
(
k
) and Z
j
(
k
) as X
(
k
), we obtain the gradient iterative algorithm for solving (39) as follows:
(49)
X
(
k
)
=
∑
i
=
1
p
Y
i
(
k
)
+
∑
j
=
1
q
Z
j
(
k
)
p
+
q
,
(50)
Y
i
(
k
)
=
X
(
k
-
1
)
+
μ
A
i
T
[
F
1
-
∑
r
=
1
p
‍
A
r
X
(
k
-
1
)
B
r
]
B
i
T
,
(51)
Z
j
(
k
)
=
X
(
k
-
1
)
+
μ
C
j
T
[
F
2
-
∑
s
=
1
q
‍
C
s
X
(
k
-
1
)
D
s
]
D
j
T
,
(52)
0
<
μ
<
2
μ
0
,
μ
0
:
=
max
⁡
⁡
{
∑
i
=
1
p
‍
∥
A
i
∥
2
∥
B
i
∥
2
,
∑
j
=
1
q
‍
∥
C
j
∥
2
∥
D
j
∥
2
}
.Theorem 6.
If (39) has a unique solution X, then the iterative solution X
(
k
) given by (49)–(52) converges to X, that is, lim
⁡
k
→
∞
X
(
k
)
=
X, or the error matrix X
(
k
)
-
X converges to zero for any initial value X
(
0
).Proof.
The estimation error matrices are defined as(53)
Y
~
i
(
k
)
≔
Y
i
(
k
)
-
X
,
i
=
1,2
,
…
,
p
,
Z
~
j
(
k
)
≔
Z
j
(
k
)
-
X
,
j
=
1,2
,
…
,
q
,
X
~
(
k
)
≔
X
(
k
)
-
X
.
By using (49)–(51), we have
(54)
Y
~
i
(
k
)
=
X
~
(
k
-
1
)
-
μ
A
i
T
[
∑
r
=
1
p
‍
A
r
X
~
(
k
-
1
)
B
r
]
B
i
T
,
000000000000000000000000000000
i
=
1,2
,
…
,
p
,
Z
~
j
(
k
)
=
X
~
(
k
-
1
)
-
μ
C
j
T
[
∑
s
=
1
q
‍
C
s
X
~
(
k
-
1
)
D
s
]
D
j
T
,
000000000000000000000000000000
j
=
1,2
,
…
,
q
,
(55)
X
~
(
k
)
=
∑
i
=
1
p
Y
~
i
(
k
)
+
∑
j
=
1
q
Z
~
i
(
k
)
p
+
q
.
Set
(56)
ξ
~
(
k
-
1
)
:
=
∑
r
=
1
p
‍
A
r
X
~
(
k
-
1
)
B
r
,
η
~
(
k
-
1
)
:
=
∑
s
=
1
q
‍
C
s
X
~
(
k
-
1
)
D
s
.
According to the trace formula ∥
X
∥
2
=
tr
⁡
[
X
T
X
] and using (54), for i
=
1,2
,
…
,
p and j
=
1,2
,
…
,
q, we obtain
(57)
∥
Y
~
i
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
(
k
-
1
)
T
A
i
T
ξ
~
(
k
-
1
)
B
i
T
]
+
μ
2
∥
A
i
T
ξ
~
(
k
-
1
)
B
i
T
∥
2
,
(58)
∥
Z
~
j
(
k
)
∥
2
=
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
tr
⁡
[
X
~
(
k
-
1
)
T
C
j
T
η
~
(
k
-
1
)
D
j
T
]
+
μ
2
∥
C
j
T
η
~
(
k
-
1
)
D
j
T
∥
2
.
For i
=
1,2
,
…
,
p and j
=
1,2
,
…
,
q, adding (57) to (58) gives
(59)
∑
i
=
1
p
‍
∥
Y
~
i
(
k
)
∥
2
+
∑
j
=
1
q
∥
Z
~
i
(
k
)
∥
2
=
(
p
+
q
)
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
(
∥
ξ
~
(
k
-
1
)
∥
2
+
∥
η
~
(
k
-
1
)
∥
2
)
+
μ
2
∑
i
=
1
p
‍
∥
A
i
T
ξ
~
(
k
-
1
)
B
i
T
∥
2
+
μ
2
∑
j
=
1
q
‍
∥
C
j
T
η
~
(
k
-
1
)
D
j
T
∥
2
⩽
(
p
+
q
)
∥
X
~
(
k
-
1
)
∥
2
-
2
μ
(
∥
ξ
~
(
k
-
1
)
∥
2
+
∥
η
~
(
k
-
1
)
∥
2
)
+
μ
2
∥
ξ
~
(
k
-
1
)
∥
2
∑
i
=
1
p
‍
∥
A
i
∥
2
∥
B
i
∥
2
+
μ
2
∥
η
~
(
k
-
1
)
∥
2
∑
j
=
1
q
‍
∥
C
j
∥
2
∥
D
j
∥
2
=
(
p
+
q
)
∥
X
~
(
k
-
1
)
∥
2
-
μ
(
2
-
μ
∑
i
=
1
p
‍
∥
A
i
∥
2
∥
B
i
∥
2
)
∥
ξ
~
(
k
-
1
)
∥
2
-
μ
(
2
-
μ
∑
j
=
1
q
‍
∥
C
j
∥
2
∥
D
j
∥
2
)
∥
η
~
(
k
-
1
)
∥
2
.
Taking the norm of both sides of (55) and using the norm inequality ∥
∑
i
=
1
m
‍
A
i
∥
2
⩽
m
∑
i
=
1
m
‍
∥
A
i
∥
2 gives
(60)
∥
X
~
(
k
)
∥
2
=
∥
∑
i
=
1
p
Y
~
i
(
k
)
+
∑
j
=
1
q
Z
~
j
(
k
)
∥
2
(
p
+
q
)
2
⩽
∑
i
=
1
p
∥
Y
~
i
(
k
)
∥
2
+
∑
j
=
1
q
∥
Z
~
j
(
k
)
∥
2
p
+
q
.
Dividing (59) by p
+
q and using the inequality (60) gives
(61)
∥
X
~
(
k
)
∥
2
⩽
∥
X
~
(
k
-
1
)
∥
2
-
μ
p
+
q
(
2
-
μ
∑
i
=
1
p
‍
∥
A
i
∥
2
∥
B
i
∥
2
)
∥
ξ
~
(
k
-
1
)
∥
2
-
μ
p
+
q
(
2
-
μ
∑
j
=
1
q
‍
∥
C
j
∥
2
∥
D
j
∥
2
)
∥
η
~
(
k
-
1
)
∥
2
⩽
∥
X
~
(
k
-
1
)
∥
2
-
μ
p
+
q
(
2
-
μ
μ
0
)
×
(
∥
ξ
~
(
k
-
1
)
∥
2
+
∥
η
~
(
k
-
1
)
∥
2
)
⩽
∥
X
~
(
0
)
∥
2
-
μ
p
+
q
(
2
-
μ
μ
0
)
∑
t
=
0
k
-
1
(
∥
ξ
~
(
t
)
∥
2
+
∥
η
~
(
t
)
∥
2
)
.
In the second “⩽,” the assumption that μ
0
=
max
⁡
{
∑
i
=
1
p
‍
∥
A
i
∥
2
∥
B
i
∥
2
,
∑
j
=
1
q
‍
∥
C
j
∥
2
∥
D
j
∥
2
} is used. Since ∥
X
~
(
k
)
∥
2
⩾
0, we have
(62)
μ
p
+
q
(
2
-
μ
μ
0
)
∑
t
=
0
k
-
1
‍
(
∥
ξ
~
(
t
)
∥
2
+
∥
η
~
(
t
)
∥
2
)
⩽
∥
X
~
(
0
)
∥
2
<
∞
.
If the convergence factor μ is chosen to satisfy 0
<
μ
<
2
/
μ
0, then we have
(63)
∑
k
=
0
∞
‍
(
∥
ξ
~
(
k
)
∥
2
+
∥
η
~
(
k
)
∥
2
)
⩽
∥
X
~
(
0
)
∥
2
.
It follows that ∥
η
~
(
k
)
∥
2
+
∥
η
~
(
k
)
∥
2
→
0, as k
→
∞ or as k
→
∞, we obtain
(64)
A
1
X
~
(
k
-
1
)
B
1
+
A
2
X
~
(
k
-
1
)
B
2
+
⋯
+
A
p
X
~
(
k
-
1
)
B
p
⟶
0
,
C
1
X
~
(
k
-
1
)
D
1
+
C
2
X
~
(
k
-
1
)
D
2
+
⋯
+
C
q
X
~
(
k
-
1
)
D
q
⟶
0
.
According to Theorem 5, we have X
~
(
k
)
→
0, as k
→
∞. This completes the proof of Theorem 6.


## 5. A Numerical Example

This section offers a numerical example to illustrate the performance of the proposed algorithm. Consider (6) with
(65)
A
1
=
[
-
0.8
-
0.3
-
0.2
-
2.1
]
,
B
1
=
[
-
1.7
0.4
-
0.8
-
1.5
]
,
A
2
=
[
0.9
-
1.5
0.6
-
0.2
]
,
B
2
=
[
0.3
1.4
1.1
0.9
]
,
C
1
=
[
-
2.4
1.3
0.4
-
0.2
]
,
D
1
=
[
0.6
-
0.3
0.6
1
]
,
C
2
=
[
0.9
0.6
-
0.4
-
0.1
]
,
D
2
=
[
0
1.2
0.1
-
1.5
]
,
F
1
=
[
-
24.34
-
38.9
25.73
1.74
]
,
F
2
=
[
24.18
31.11
-
3.7
-
7.52
]
.
From Theorem 3, the exact solution is
(66)
X
=
[
x
11
x
12
x
21
x
22
]
=
[
-
2
-
8
7
6
]
.TakingX
(
0
)
=
1
0
-
6
1
2
×
2 as initial iterative value, we apply the gradient-based iterative algorithm in (15)–(19) to compute X
(
k
). The iterative solutions X
(
k
) are shown in Table 1, where δ
=
∥
X
(
k
)
-
X
∥
/
∥
X
∥ is the relative error and the convergence factor μ is 25
/
139. The relation of the relative error δ with different convergence factors is shown in Figure 1. From Table 1 and Figure 1, we can find that δ becomes smaller and smaller and tends to zero as the iterative times increase. This demonstrates that the algorithm proposed in this paper is effective.Table 1
The gradient iterative solution (μ
=
25
/
139).
k
x
11
x
12
x
21
x
22
δ (%)
2
−1.60864
−7.04633
6.47331
5.28946
10.98102
4
−1.89505
−7.89255
6.97696
5.89765
1.48121
6
−1.96855
−7.98772
7.00523
5.98269
0.30968
8
−1.99026
−7.99846
7.00259
5.99669
0.08667
10
−1.99696
−7.99976
7.00084
5.99931
0.02617
12
−1.99905
−7.99995
7.00025
5.99985
0.00805
14
−1.99970
−7.99999
7.00007
5.99996
0.00249
16
−1.99991
−8.00000
7.00002
5.99999
0.00077
18
−1.99997
−8.00000
7.00001
6.00000
0.00024
20
−1.99999
−8.00000
7.00000
6.00000
0.00007
Solution
−2.00000
−8.00000
7.00000
6.00000Figure 1
The relative errorδ versus k  (dots: μ
0
=
max
⁡
⁡
{
∑
i
=
1
2
‍
∥
A
i
∥
2
∥
B
i
∥
2
,
∑
j
=
1
2
‍
∥
C
j
∥
2
∥
D
j
∥
2
}
=
11.12).A simple calculation indicates that the range of the convergence factor is conservative. Based on the deep analysis of Figure1, we find that the rate of convergence increases when μ enlarges from 25
/
556 to 25
/
278. However, if we keep enlarging μ from 25
/
139 to 125
/
556, the rate of convergence will drop. This shows that the best convergence factor is uncovered in this algorithm. How to find the best convergence factor subtly is our work in the future.


## 6. Conclusions

This paper has proposed a gradient-based iterative algorithm for solving a class of the real coupled matrix equations. By using the hierarchical identification principle, we prove that the iterative solution is convergent if the unique solution exists. An example demonstrates that the algorithm is effective and indicates that the best convergence factor is existent. We will find the best convergence factor of this algorithm in the future.



---

*Source: 649524-2014-05-04.xml*