# Riemannian Gradient Algorithm for the Numerical Solution of Linear Matrix Equations

Authors: Xiaomin Duan; Huafei Sun; Xinyu Zhao
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

A Riemannian gradient algorithm based on geometric structures of a manifold consisting of all positive definite matrices is proposed to calculate the numerical solution of the linear matrix equationQ
=
X
+
∑
i
=
1
m
A
i
T
X
A
i. In this algorithm, the geodesic distance on the curved Riemannian manifold is taken as an objective function and the geodesic curve is treated as the convergence path. Also the optimal variable step sizes corresponding to the minimum value of the objective function are provided in order to improve the convergence speed. Furthermore, the convergence speed of the Riemannian gradient algorithm is compared with that of the traditional conjugate gradient method in two simulation examples. It is found that the convergence speed of the provided algorithm is faster than that of the conjugate gradient method.


---

## Body

## 1. Introduction

The linear matrix equation(1)
Q
=
X
+
∑
i
=
1
m
‍
A
i
T
X
A
i
,
where A
1
,
A
2
,
…
,
A
m are arbitrary n
×
n real matrices, m is a nonnegative integer, and A
i
T denotes the transpose of the matrix A
i, arises from many fields, such as the control theory, the dynamic programming, and the stochastic filtering [1–4]. In the past decades, there has been increasing interest in the solution problems of this equation. In the case of m
=
1, some numerical methods, such as Bartels-Stewart method, Hessenberg-Schur method, and Schur and QR decompositions method, were proposed in [5, 6]. Based on the Kronecker product and the fixed point theorem in partially ordered sets, some sufficient conditions for the existence of a unique symmetric positive definite solution are given in [7, 8]. Ran and Reurings ([7, Theorem 3.3] and [9, Theorem 3.1]) pointed out that if Q
-
∑
i
=
1
m
‍
A
i
T
Q
A
i is a positive definite matrix, then there exists a unique solution and it is symmetric positive definite. Recently, under the condition that (1) is consistent, Su and Chen presented an efficient numerical iterative method based on the conjugate gradient method (CGM) [10].In addition, based on geometric structures on a Riemannian manifold, Duan et al. proposed a natural gradient descent algorithm to solve algebraic Lyapunov equations [11, 12]. Following them, we investigate the solution problem of (1) in the view of Riemannian manifolds. Note that this solution of (1) is a symmetric positive definite matrix and the set of all the symmetric positive definite matrices can be considered as a manifold. Thus, it is more convenient to investigate the solution problem with the help of these geometric structures on this manifold. To address such a need, a new framework is presented in this paper to calculate the numerical solution, which is based on the geometric structures on the Riemannian manifold of positive definite symmetric matrices. We will first describe briefly some fundamental knowledge on manifold. Then a Riemannian gradient algorithm is proposed, where the geodesic distance on the curved Riemannian manifold is taken as an objective function and the geodesic curve is treated as the convergence path. Also the optimal variable step sizes corresponding to the minimum value of the objective function are provided in order to improve the convergence speed. Finally, the behavior of the provided algorithm and the traditional CGM is compared and demonstrated in two simulation examples.


## 2. A Survey of Some Geometrical Concepts

In the present section, we briefly survey the geometry of smooth manifolds by recalling concepts as tangent spaces, the Riemannian gradient, geodesic curves, and geodesic distance. More details can be found in [13, 14]. Specially, these concepts of manifold consisting of all symmetric positive definite matrices have also been introduced in [15–20].


### 2.1. Riemannian Manifolds

Let us denote a Riemannian manifold by𝕄. Its tangent space at point Y
∈
𝕄 is denoted by T
Y
𝕄. Given any pair of points V
,
W
∈
T
Y
𝕄, an inner product 〈
V
,
W
〉
Y
∈
ℝ is defined. The specification of the inner product for a Riemannian manifold turns it into a metric space. In fact, the length of a curve ϕ
Y
,
V
:
[
0,1
]
→
𝕄, such that ϕ
Y
,
V
(
0
)
=
Y
∈
𝕄 and ϕ
˙
Y
,
V
(
0
)
=
V
∈
T
Y
𝕄, is given by
(2)
ℓ
(
ϕ
Y
,
V
)
:
=
∫
0
1
‍
〈
ϕ
˙
Y
,
V
(
t
)
,
ϕ
˙
Y
,
V
(
t
)
〉
ϕ
Y
,
V
(
t
)
d
t
.
Given arbitrary Y
∈
𝕄 and V
∈
T
Y
𝕄, the curve γ
Y
,
V
:
[
0,1
]
→
𝕄 of shortest length is called geodesic curve. Such minimal length ℓ
(
γ
Y
,
V
) is called geodesic distance between two points, namely, γ
Y
,
V
(
0
)
=
Y and γ
Y
,
V
(
1
). The Riemannian distance between two points is denoted by
(3)
d
(
γ
Y
,
V
(
0
)
,
γ
Y
,
V
(
1
)
)
:
=
ℓ
(
γ
Y
,
V
)
.
It is obvious that if any pair of points on manifold 𝕄 can be connected by a geodesic curve, then it is possible to measure the distance between any given pair of points on it. Given a regular function f
:
𝕄
→
ℝ, its Riemannian gradient ∇
Y
f in the direction of the vector V
∈
T
Y
𝕄 measures the rate of change of the function f in the direction V. Namely, given any smooth curve ϕ
Y
,
V
:
[
0,1
]
→
𝕄, such that ϕ
Y
,
V
(
0
)
=
Y
∈
𝕄 and ϕ
˙
Y
,
V
(
0
)
=
V, the Riemannian gradient ∇
Y
f is the unique vector in T
Y
𝕄 such that
(4)
〈
V
,
∇
Y
f
〉
Y
=
d
d
t
f
(
ϕ
Y
,
V
(
t
)
)
|
t
=
0
.Instead of the Euclidean space, by using the Riemannian gradient the optimization method may be readily extended to smooth manifolds [21, Chapter 7]. For this purpose, let us consider the differential equation on manifold 𝕄
(5)
y
˙
=
-
∇
y
(
t
)
f
,
y
(
0
)
=
y
¯
∈
𝕄
.If the functionf
:
𝕄
→
ℝ is bounded and the smooth manifold 𝕄 is compact, then the solution of such differential equation tends to a local minimum of the function f in 𝕄, depending on the initial point y
¯. In fact, it is achieved by definition of the Riemannian gradient
(6)
d
d
t
f
(
y
(
t
)
)
=
〈
y
˙
(
t
)
,
∇
y
(
t
)
f
〉
y
(
t
)
=
-
〈
∇
y
(
t
)
f
,
∇
y
(
t
)
f
〉
y
(
t
)
≤
0
.


### 2.2. Manifold of Symmetric Positive Definite Matrices

Let us denote the general linear group by𝔾
𝕃
(
n
), which consists of all nonsingular real matrices. And X
>
0 means that X is a symmetric positive definite matrix. For a different notation X
-
Y
>
0, we will use X
>
Y. And the manifold consisting of all symmetric positive definite matrices is defined by 𝕊
+
(
n
)
=
{
X
∈
ℝ
n
×
n
∣
X
T
=
X
,
X
>
0
}. The tangent space at a point X
∈
𝕊
+
(
n
) is given by T
X
𝕊
+
(
n
)
=
{
V
∈
ℝ
n
×
n
∣
V
T
=
V
}. On manifold 𝕊
+
(
n
), the Riemannian metric at the point X is defined by
(7)
〈
V
1
,
V
2
〉
X
=
tr
[
X
-
1
V
1
X
-
1
V
2
]
,
where V
1
,
V
2
∈
T
X
𝕊
+
(
n
) and tr denotes the trace of the matrix. Then, manifold 𝕊
+
(
n
) with the Riemannian metric (7) becomes a Riemannian manifold. Moreover, it is a curved manifold since this metric is invariant under the group action of 𝔾
𝕃
(
n
). The geodesic curve at the point X
∈
𝕊
+
(
n
) in the direction V
∈
T
X
𝕊
+
(
n
) can be expressed as
(8)
γ
X
,
V
(
t
)
=
X
1
/
2
exp
(
t
X
-
1
/
2
V
X
-
1
/
2
)
X
1
/
2
.
Obviously, this geodesic curve is entirely contained in manifold. The Hopf-Rinow theorem implies that manifold 𝕊
+
(
n
) with the Riemannian metric (7) is geodesically complete ([22, Theorem 2.8]). This means, for any given pair X
,
Y
∈
𝕊
+
(
n
), we can find a geodesic curve γ
X
,
V
(
t
), such that γ
X
,
V
(
0
)
=
X and γ
X
,
V
(
1
)
=
Y, by taking the initial velocity γ
˙
X
,
V
(
0
)
=
X
1
/
2
log
(
X
-
1
/
2
Y
X
-
1
/
2
)
X
1
/
2
∈
T
X
𝕊
+
(
n
). According to (2), the geodesic distance d
(
X
,
Y
) can be computed explicitly by
(9)
d
2
(
X
,
Y
)
=
tr
[
log
2
(
X
-
1
/
2
Y
X
-
1
/
2
)
]
.


## 2.1. Riemannian Manifolds

Let us denote a Riemannian manifold by𝕄. Its tangent space at point Y
∈
𝕄 is denoted by T
Y
𝕄. Given any pair of points V
,
W
∈
T
Y
𝕄, an inner product 〈
V
,
W
〉
Y
∈
ℝ is defined. The specification of the inner product for a Riemannian manifold turns it into a metric space. In fact, the length of a curve ϕ
Y
,
V
:
[
0,1
]
→
𝕄, such that ϕ
Y
,
V
(
0
)
=
Y
∈
𝕄 and ϕ
˙
Y
,
V
(
0
)
=
V
∈
T
Y
𝕄, is given by
(2)
ℓ
(
ϕ
Y
,
V
)
:
=
∫
0
1
‍
〈
ϕ
˙
Y
,
V
(
t
)
,
ϕ
˙
Y
,
V
(
t
)
〉
ϕ
Y
,
V
(
t
)
d
t
.
Given arbitrary Y
∈
𝕄 and V
∈
T
Y
𝕄, the curve γ
Y
,
V
:
[
0,1
]
→
𝕄 of shortest length is called geodesic curve. Such minimal length ℓ
(
γ
Y
,
V
) is called geodesic distance between two points, namely, γ
Y
,
V
(
0
)
=
Y and γ
Y
,
V
(
1
). The Riemannian distance between two points is denoted by
(3)
d
(
γ
Y
,
V
(
0
)
,
γ
Y
,
V
(
1
)
)
:
=
ℓ
(
γ
Y
,
V
)
.
It is obvious that if any pair of points on manifold 𝕄 can be connected by a geodesic curve, then it is possible to measure the distance between any given pair of points on it. Given a regular function f
:
𝕄
→
ℝ, its Riemannian gradient ∇
Y
f in the direction of the vector V
∈
T
Y
𝕄 measures the rate of change of the function f in the direction V. Namely, given any smooth curve ϕ
Y
,
V
:
[
0,1
]
→
𝕄, such that ϕ
Y
,
V
(
0
)
=
Y
∈
𝕄 and ϕ
˙
Y
,
V
(
0
)
=
V, the Riemannian gradient ∇
Y
f is the unique vector in T
Y
𝕄 such that
(4)
〈
V
,
∇
Y
f
〉
Y
=
d
d
t
f
(
ϕ
Y
,
V
(
t
)
)
|
t
=
0
.Instead of the Euclidean space, by using the Riemannian gradient the optimization method may be readily extended to smooth manifolds [21, Chapter 7]. For this purpose, let us consider the differential equation on manifold 𝕄
(5)
y
˙
=
-
∇
y
(
t
)
f
,
y
(
0
)
=
y
¯
∈
𝕄
.If the functionf
:
𝕄
→
ℝ is bounded and the smooth manifold 𝕄 is compact, then the solution of such differential equation tends to a local minimum of the function f in 𝕄, depending on the initial point y
¯. In fact, it is achieved by definition of the Riemannian gradient
(6)
d
d
t
f
(
y
(
t
)
)
=
〈
y
˙
(
t
)
,
∇
y
(
t
)
f
〉
y
(
t
)
=
-
〈
∇
y
(
t
)
f
,
∇
y
(
t
)
f
〉
y
(
t
)
≤
0
.


## 2.2. Manifold of Symmetric Positive Definite Matrices

Let us denote the general linear group by𝔾
𝕃
(
n
), which consists of all nonsingular real matrices. And X
>
0 means that X is a symmetric positive definite matrix. For a different notation X
-
Y
>
0, we will use X
>
Y. And the manifold consisting of all symmetric positive definite matrices is defined by 𝕊
+
(
n
)
=
{
X
∈
ℝ
n
×
n
∣
X
T
=
X
,
X
>
0
}. The tangent space at a point X
∈
𝕊
+
(
n
) is given by T
X
𝕊
+
(
n
)
=
{
V
∈
ℝ
n
×
n
∣
V
T
=
V
}. On manifold 𝕊
+
(
n
), the Riemannian metric at the point X is defined by
(7)
〈
V
1
,
V
2
〉
X
=
tr
[
X
-
1
V
1
X
-
1
V
2
]
,
where V
1
,
V
2
∈
T
X
𝕊
+
(
n
) and tr denotes the trace of the matrix. Then, manifold 𝕊
+
(
n
) with the Riemannian metric (7) becomes a Riemannian manifold. Moreover, it is a curved manifold since this metric is invariant under the group action of 𝔾
𝕃
(
n
). The geodesic curve at the point X
∈
𝕊
+
(
n
) in the direction V
∈
T
X
𝕊
+
(
n
) can be expressed as
(8)
γ
X
,
V
(
t
)
=
X
1
/
2
exp
(
t
X
-
1
/
2
V
X
-
1
/
2
)
X
1
/
2
.
Obviously, this geodesic curve is entirely contained in manifold. The Hopf-Rinow theorem implies that manifold 𝕊
+
(
n
) with the Riemannian metric (7) is geodesically complete ([22, Theorem 2.8]). This means, for any given pair X
,
Y
∈
𝕊
+
(
n
), we can find a geodesic curve γ
X
,
V
(
t
), such that γ
X
,
V
(
0
)
=
X and γ
X
,
V
(
1
)
=
Y, by taking the initial velocity γ
˙
X
,
V
(
0
)
=
X
1
/
2
log
(
X
-
1
/
2
Y
X
-
1
/
2
)
X
1
/
2
∈
T
X
𝕊
+
(
n
). According to (2), the geodesic distance d
(
X
,
Y
) can be computed explicitly by
(9)
d
2
(
X
,
Y
)
=
tr
[
log
2
(
X
-
1
/
2
Y
X
-
1
/
2
)
]
.


## 3. Riemannian Gradient Algorithm

Some researches show that the solution equation (1) is unique positive definite if A
i and Q satisfy ([7, Theorem 3.3] and [9, Theorem 3.1])
(10)
Q
>
∑
i
=
1
m
‍
A
i
T
Q
A
i
.
For convenience, we set ℒ
(
X
)
:
=
X
+
∑
i
=
1
m
‍
A
i
T
X
A
i. And it is obvious that ℒ
(
X
) is positive definite for any X
∈
𝕊
+
(
n
). Our purpose is to seek a matrix X
∈
𝕊
+
(
n
) so that ℒ
(
X
) is as close as possible to the given matrix Q. We figure out some key points in designing such an algorithm as follows.  (1)
For someX, the difference between the given positive definite matrix Q and the matrix ℒ
(
X
) should be determined and calculated firstly.
(2)
Then the Riemannian gradient descent algorithm can be used to adjustX so that the difference will be as small as possible.Note that the Riemannian manifold𝕊
+
(
n
) is geodesically complete; hence we take the geodesic distance between Q and ℒ
(
X
) as the objective function J
(
X
); namely,
(11)
J
(
X
)
=
d
2
(
Q
,
ℒ
(
X
)
)
.Then the exact solution of (1) may be defined as
(12)
X
*
=
arg
min
X
∈
𝕊
+
(
n
)
J
(
X
)
.
In order to find a minimizer for the criterion d
2
(
Q
,
ℒ
(
X
)
), we may make use of the differential equation (5). The minimum is apparently achieved for a value X
* such that
(13)
∇
X
J
(
X
)
|
X
=
X
*
=
0
.
It is necessary to solve the optimization problem (12) by a numerical optimization algorithm that takes the geometry of the Riemannian manifold 𝕊
+
(
n
) into account. In particular, starting from an initial guess X
0, it is possible to provide an iterative learning algorithm that generates a pattern X
t
∈
𝕊
+
(
n
) at any learning step t
∈
ℕ. Following [21, 23], such sequence may be generated by moving from each point X
t
∈
𝕊
+
(
n
) to the next point X
t
+
1 along a short geodesic curve in the opposite direction of the Riemannian gradient of the objective function J
(
X
t
). Namely, if we set
(14)
V
t
=
∇
X
t
J
(
X
t
)
,
then the Riemannian gradient algorithm may be expressed as
(15)
X
t
+
1
=
γ
X
t
,
-
V
t
(
η
)
,
where η denotes any suitable step-size schedule that drives the iterative algorithm (15) to be convergent. Moreover, with the aim to employ the learning algorithm (14) and (15), we first need to compute the Riemannian gradient of the objective function J
(
X
t
) to be optimized. In the definition of the Riemannian gradient (4), the generic smooth curve ϕ
Y
,
V may be replaced with a geodesic curve. We obtain the following theorem.Theorem 1.
The Riemannian gradient of the objective functionJ
(
X
) at point X is given by
(16)
∇
X
J
(
X
)
=
2
X
(
∑
j
=
1
m
‍
log
(
Q
-
1
ℒ
(
X
)
)
ℒ
(
X
)
-
1
h
h
h
h
h
+
∑
j
=
1
m
‍
A
j
log
(
Q
-
1
ℒ
(
X
)
)
ℒ
(
X
)
-
1
A
j
T
)
X
.Proof.
First, we compute the derivative of the real-valued function(17)
J
(
s
(
t
)
)
=
tr
[
log
2
(
Q
-
1
/
2
ℒ
(
s
(
t
)
)
Q
-
1
/
2
)
]
with respect to t, where s
(
t
)
=
X
1
/
2
exp
(
t
X
-
1
/
2
V
X
-
1
/
2
)
X
1
/
2 is the geodesic curve emanating from X in the direction V
=
s
˙
(
0
). Using Proposition 2.1 in [16] and some properties of the trace, we have
(18)
d
d
t
J
(
s
(
t
)
)
|
t
=
0
=
2
tr
[
(
∑
j
=
1
m
‍
log
(
Q
-
1
ℒ
(
X
)
)
ℒ
(
X
)
-
1
h
h
h
h
h
h
h
h
+
∑
j
=
1
m
‍
A
j
log
(
Q
-
1
ℒ
(
X
)
)
ℒ
(
X
)
-
1
A
j
T
)
V
]
=
〈
V
,
2
X
(
∑
j
=
1
m
‍
log
(
Q
-
1
ℒ
(
X
)
)
ℒ
(
X
)
-
1
h
h
h
h
h
h
h
h
h
h
+
∑
j
=
1
m
‍
A
j
log
(
Q
-
1
ℒ
(
X
)
)
ℒ
(
X
)
-
1
A
j
T
)
X
〉
X
.
From (4), it is shown that formula (16) is valid.Remark 2.
In fact, note thatlog
(
A
-
1
B
A
)
=
A
-
1
(
log
B
)
A, for all A
∈
𝔾
𝕃
(
n
) and for all B
∈
𝕊
+
(
n
); then ∇
X
J
(
X
) can be rewritten as
(19)
∇
X
J
(
X
)
=
2
X
(
ℒ
(
X
)
-
1
/
2
∑
j
=
1
m
‍
log
(
ℒ
(
X
)
1
/
2
Q
-
1
ℒ
(
X
)
1
/
2
)
ℒ
(
X
)
-
1
/
2
h
h
h
h
h
h
h
h
+
∑
j
=
1
m
‍
A
j
ℒ
(
X
)
-
1
/
2
h
h
h
h
h
h
h
h
∑
j
=
1
m
‍
×
log
(
ℒ
(
X
)
1
/
2
Q
-
1
ℒ
(
X
)
1
/
2
)
ℒ
(
X
)
-
1
/
2
A
j
T
)
X
.
Therefore it can be shown that ∇
X
J
(
X
) belongs to T
X
𝕊
+
(
n
).Corollary 3.
The Riemannian gradient∇
X
J
(
X
) has the unique zero point X
*. Furthermore, for the suitable step-size η, the sequence X
t generated by (15) converges to the solution X
* of (1).Proof.
Obviously,X
* is a zero point of the Riemannian gradient ∇
X
J
(
X
). If there exists X
*
′
≠
X
* such that the iterative process (15) stops, then it means ∇
X
*
′
J
(
X
*
′
)
=
0, which is equivalent to
(20)
log
(
Q
-
1
ℒ
(
X
*
′
)
)
ℒ
(
X
*
′
)
-
1
+
∑
j
=
1
m
‍
A
j
log
(
Q
-
1
ℒ
(
X
*
′
)
)
ℒ
(
X
*
′
)
-
1
A
j
T
=
0
.
On (20), the matrices are multiplied by X
*
′ and the traces of the matrices are taken; then we have
(21)
tr
[
log
(
Q
-
1
ℒ
(
X
*
′
)
)
]
=
0
.
Also, if we replace X
*
′ with X
* and follow the similar procedure as the above equation, then the following equalities are valid:
(22)
tr
[
log
(
Q
-
1
ℒ
(
X
*
′
)
)
ℒ
(
X
*
′
)
-
1
(
X
*
+
∑
i
=
1
m
‍
A
i
T
X
*
A
i
)
]
=
tr
[
log
(
Q
-
1
ℒ
(
X
*
′
)
)
ℒ
(
X
*
′
)
-
1
Q
]
=
0
.
Combining (21) and (22), it can be obtained that
(23)
tr
[
log
(
Q
-
1
ℒ
(
X
*
′
)
)
]
-
tr
[
log
(
Q
-
1
ℒ
(
X
*
′
)
)
Q
]
=
tr
[
log
(
Q
-
1
ℒ
(
X
*
′
)
)
(
I
-
ℒ
(
X
*
′
)
-
1
Q
)
]
=
tr
[
log
(
Q
-
1
/
2
ℒ
(
X
*
′
)
Q
-
1
/
2
)
×
(
I
-
Q
1
/
2
ℒ
(
X
*
′
)
-
1
Q
1
/
2
)
]
=
0
.
Note that the matrix Q
1
/
2
ℒ
(
X
*
′
)
-
1
Q
1
/
2 is positive definite; hence its orthogonal similarity diagonalization can be expressed as U
Λ
U
T, where U is an orthogonal matrix, Λ
=
diag
(
λ
1
,
…
,
λ
n
), and λ
i
>
0
(
1
≤
i
≤
n
). Thus, we have
(24)
tr
[
(
Λ
-
I
)
log
Λ
]
=
∑
i
=
1
n
‍
(
λ
i
-
1
)
log
λ
i
=
0
.
Note that the second equality of (24) is equivalent to
(25)
∏
i
=
1
n
‍
λ
i
λ
i
-
1
=
1
.
Thus we have λ
1
=
λ
2
=
⋯
=
λ
n
=
1. That means that X
*
′ is also the solution of (1), which is contrary to the uniqueness of the solution. Therefore, the Riemannian gradient ∇
X
J
(
X
) has the unique zero point X
*. This completes the proof of Corollary 3.Now, the Riemannian gradient algorithm with a constant step-size (RGACS) becomes(26)
X
t
+
1
=
X
t
1
/
2
exp
(
-
η
X
t
-
1
/
2
∇
X
t
J
(
X
t
)
X
t
-
1
/
2
)
X
t
1
/
2
,
where the initial value X
0
∈
𝕊
+
(
n
). The step-size η corresponding to the best convergence speed can be obtained experimentally and the details will be described in our following simulations.


### 3.1. Riemannian Gradient Algorithm with Variable Step Sizes

In order to improve the convergence speed of the RGACS, an optimal variable step-size scheduleη
t will be considered and defined as follows. The step-size η
t should be evaluated in such a way that the objective function J
(
X
t
+
1
) at step t
+
1 is reduced as much as possible with respect to J
(
X
t
). We may regard J
(
X
t
+
1
) as a function of the step-size η
t that
(27)
J
(
X
t
+
1
)
=
d
2
(
Q
,
ℒ
(
γ
X
t
,
-
V
t
(
η
t
)
)
)
.
Thus, we can optimize the value of η
t to ensure that the difference J
(
X
t
)
-
J
(
X
t
+
1
) could be as large as possible. Since it is difficult to solve this nonlinear problem exactly, we use a suboptimal approximation in practice. Under the hypothesis that the step-size value η
t is small enough, we may invoke the expansion of the function J
(
X
t
+
1
) around the point η
t
=
0 as follows:
(28)
J
(
X
t
+
1
)
≈
J
(
X
t
)
+
C
1
,
t
η
t
+
1
2
C
2
,
t
η
t
2
.
Then the optimal step-size η
t
*, corresponding to the maximum of the difference J
(
X
t
)
-
J
(
X
t
+
1
), can be expressed as
(29)
η
t
*
≈
-
C
1
,
t
C
2
,
t
.
The coefficients C
1
,
t, and C
2
,
t can be calculated using the first-order and the second-order derivatives of the function J
(
X
t
+
1
) with respect to the parameter η
t at the point η
t
=
0. On the basis of the above discussion, we may obtain the optimal step-size
(30)
η
t
*
≈
tr
[
(
X
ℒ
(
log
(
Q
-
1
ℒ
(
X
t
)
)
ℒ
(
X
t
)
-
1
)
)
2
]
tr
[
(
ℒ
(
X
t
)
-
1
ℒ
(
∇
X
t
J
(
X
t
)
)
)
2
]
.Then, the Riemannian gradient algorithm with variable step-sizes (RGAVS) can be written explicitly as(31)
X
t
+
1
=
X
t
1
/
2
exp
(
-
η
t
*
X
t
-
1
/
2
∇
X
t
J
(
X
t
)
X
t
-
1
/
2
)
X
t
1
/
2
,
where X
0
∈
𝕊
+
(
n
) is the initial value.


### 3.2. Volume Feature of Riemannian Gradient Algorithm

The RGACS shows an interesting volume feature that the ratiodet
(
X
t
+
1
)
/
det
(
X
t
) always keeps a certain relationship. In order to prove the volume feature, let us recall two properties of determinant operator det; namely, for each A
,
B
∈
𝔾
𝕃
(
n
), there hold
(32)
det
(
A
B
)
=
det
(
A
)
det
(
B
)
,
det
(
exp
(
A
)
)
=
exp
(
tr
(
A
)
)
.
By applying the determinant operator to both sides of (31), we obtain
(33)
det
(
X
t
+
1
)
=
det
(
X
t
)
exp
(
tr
[
-
η
X
t
-
1
/
2
∇
X
t
J
(
X
t
)
X
t
-
1
/
2
]
)
=
det
(
X
t
)
exp
(
tr
[
η
(
∑
j
=
1
m
‍
log
(
ℒ
(
X
t
)
-
1
Q
)
ℒ
(
X
t
)
-
1
X
t
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
log
(
ℒ
(
X
t
)
-
1
Q
)
ℒ
(
X
t
)
-
1
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
×
∑
j
=
1
m
‍
A
j
T
X
t
A
j
)
]
)
=
det
(
X
t
)
exp
(
tr
[
η
log
(
ℒ
(
X
t
)
-
1
Q
)
]
)
=
det
(
X
t
)
det
(
(
ℒ
(
X
t
)
-
1
Q
)
η
)
.
Therefore, the following property is valid:
(34)
det
(
X
t
+
1
)
det
(
X
t
)
=
(
det
(
Q
)
det
(
ℒ
(
X
t
)
)
η
.
This equation shows that the determinant of (1) is closely linked with that of the sequence X
t. Moreover, the ratio det
(
X
t
+
1
)
/
det
(
X
t
) is directly proportional to det
(
Q
)
/
det
(
ℒ
(
X
t
)
) for the fixed step-size η. In the case of the RGAVS, the similar conclusion can also be obtained.


## 3.1. Riemannian Gradient Algorithm with Variable Step Sizes

In order to improve the convergence speed of the RGACS, an optimal variable step-size scheduleη
t will be considered and defined as follows. The step-size η
t should be evaluated in such a way that the objective function J
(
X
t
+
1
) at step t
+
1 is reduced as much as possible with respect to J
(
X
t
). We may regard J
(
X
t
+
1
) as a function of the step-size η
t that
(27)
J
(
X
t
+
1
)
=
d
2
(
Q
,
ℒ
(
γ
X
t
,
-
V
t
(
η
t
)
)
)
.
Thus, we can optimize the value of η
t to ensure that the difference J
(
X
t
)
-
J
(
X
t
+
1
) could be as large as possible. Since it is difficult to solve this nonlinear problem exactly, we use a suboptimal approximation in practice. Under the hypothesis that the step-size value η
t is small enough, we may invoke the expansion of the function J
(
X
t
+
1
) around the point η
t
=
0 as follows:
(28)
J
(
X
t
+
1
)
≈
J
(
X
t
)
+
C
1
,
t
η
t
+
1
2
C
2
,
t
η
t
2
.
Then the optimal step-size η
t
*, corresponding to the maximum of the difference J
(
X
t
)
-
J
(
X
t
+
1
), can be expressed as
(29)
η
t
*
≈
-
C
1
,
t
C
2
,
t
.
The coefficients C
1
,
t, and C
2
,
t can be calculated using the first-order and the second-order derivatives of the function J
(
X
t
+
1
) with respect to the parameter η
t at the point η
t
=
0. On the basis of the above discussion, we may obtain the optimal step-size
(30)
η
t
*
≈
tr
[
(
X
ℒ
(
log
(
Q
-
1
ℒ
(
X
t
)
)
ℒ
(
X
t
)
-
1
)
)
2
]
tr
[
(
ℒ
(
X
t
)
-
1
ℒ
(
∇
X
t
J
(
X
t
)
)
)
2
]
.Then, the Riemannian gradient algorithm with variable step-sizes (RGAVS) can be written explicitly as(31)
X
t
+
1
=
X
t
1
/
2
exp
(
-
η
t
*
X
t
-
1
/
2
∇
X
t
J
(
X
t
)
X
t
-
1
/
2
)
X
t
1
/
2
,
where X
0
∈
𝕊
+
(
n
) is the initial value.


## 3.2. Volume Feature of Riemannian Gradient Algorithm

The RGACS shows an interesting volume feature that the ratiodet
(
X
t
+
1
)
/
det
(
X
t
) always keeps a certain relationship. In order to prove the volume feature, let us recall two properties of determinant operator det; namely, for each A
,
B
∈
𝔾
𝕃
(
n
), there hold
(32)
det
(
A
B
)
=
det
(
A
)
det
(
B
)
,
det
(
exp
(
A
)
)
=
exp
(
tr
(
A
)
)
.
By applying the determinant operator to both sides of (31), we obtain
(33)
det
(
X
t
+
1
)
=
det
(
X
t
)
exp
(
tr
[
-
η
X
t
-
1
/
2
∇
X
t
J
(
X
t
)
X
t
-
1
/
2
]
)
=
det
(
X
t
)
exp
(
tr
[
η
(
∑
j
=
1
m
‍
log
(
ℒ
(
X
t
)
-
1
Q
)
ℒ
(
X
t
)
-
1
X
t
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
log
(
ℒ
(
X
t
)
-
1
Q
)
ℒ
(
X
t
)
-
1
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
×
∑
j
=
1
m
‍
A
j
T
X
t
A
j
)
]
)
=
det
(
X
t
)
exp
(
tr
[
η
log
(
ℒ
(
X
t
)
-
1
Q
)
]
)
=
det
(
X
t
)
det
(
(
ℒ
(
X
t
)
-
1
Q
)
η
)
.
Therefore, the following property is valid:
(34)
det
(
X
t
+
1
)
det
(
X
t
)
=
(
det
(
Q
)
det
(
ℒ
(
X
t
)
)
η
.
This equation shows that the determinant of (1) is closely linked with that of the sequence X
t. Moreover, the ratio det
(
X
t
+
1
)
/
det
(
X
t
) is directly proportional to det
(
Q
)
/
det
(
ℒ
(
X
t
)
) for the fixed step-size η. In the case of the RGAVS, the similar conclusion can also be obtained.


## 4. Simulations

Two simulation examples are given to compare the convergence speed of the Riemannian gradient algorithm (the RGACS and the RGAVS) with those of the traditional CGM. The error criterion used in these simulations isρ
(
ℒ
(
X
)
-
Q
)
<
1
0
-
15, where ρ
(
·
) denotes the spectral radius of the matrix.Example 4.
First we consider the matrices equation(35)
Q
=
X
+
A
1
T
X
A
1
+
A
2
T
X
A
2
+
A
3
T
X
A
3
in the case of n
=
5. And Q, A
1, A
2, A
3 are as follows:
(36)
Q
=
(
4.6350
0.0698
0.3104
0.1990
0.2918
0.0698
4.0376
0.0742
0.0694
0.0690
0.3104
0.0742
4.2242
0.1554
0.1826
0.1990
0.0694
0.1554
4.1472
0.1572
0.2918
0.0690
0.1826
0.1572
4.2134
)
,
A
1
=
(
0.0834
0.0007
0.0293
0.0142
0.0139
0.0141
0.2975
0.0154
0.0289
0.0339
0.0244
0.0152
0.0317
0.0270
0.0117
0.0251
0.0097
0.0140
0.0135
0.0240
0.0158
0.1499
0.0275
0.0077
0.0157
)
,
A
2
=
(
0.2624
0.0166
0.0454
0.0172
0.0160
0.0893
0.0608
0.0242
0.0029
0.0061
0.0054
0.0190
0.0589
0.0036
0.0160
0.0278
0.0050
0.0255
0.0689
0.0048
0.0319
0.0064
0.0103
0.0219
0.0018
)
,
A
3
=
(
8.1622
0.6020
0.4505
0.8258
0.1067
2.7943
7.0630
0.0838
0.5383
0.9619
0.3112
0.6541
6.2290
0.9961
0.0046
0.5285
0.6892
0.9133
1.0782
0.7749
0.1656
0.7482
0.1524
0.4427
9.0073
)
,
which satisfy the condition (10). Both the RGACS and RGAVS are used and compared with the traditional CGM, and Q is taken as the initial value X
0 in the following simulation. In the RGACS case, the step sizes versus the iterations are shown in Figure 1. It can be found that the iterations will reduce gradually as the step-size changes from zero to a critical value; then the iterations will increase gradually above this critical value. Therefore, the critical value about 0.64 can be selected experimentally as an optimal step-size and used in the RGACS.Figure 1
Step sizes versus the iterations in the RGACS.Figure2 shows the convergence comparison between the RGACS and the RGAVS and the traditional CGM. It can be found that the RGAVS possesses the most steady convergence and the fastest convergence speed, and it needs 23 iterations to obtain the numerical solution of (35) as follows:
(37)
(
0.0856
-
0.0265
0.0151
-
0.1065
0.0148
-
0.0265
0.0970
0.0166
-
0.1581
-
0.0023
0.0151
0.0166
0.1591
-
0.3608
0.0283
-
0.1065
-
0.1581
-
0.3608
2.3761
-
0.1723
0.0148
-
0.0023
0.0283
-
0.1723
0.0654
)
.
The case of the CGM realizes the given error through 33 iterations and the slowest one among them is the RGACS with 31 steps. By comparison, it is shown that the convergence speed of the RGAVS is faster than both those of the RGACS and the traditional CGM in solving (35).Figure 2
Comparison of convergence speeds (log scale).Example 5.
A 20-order linear matrix equation withm
=
2 is also considered here following the similar procedure as Example 4, where Q
∈
𝕊
+
(
20
), and A
1
,
A
2
∈
ℝ
20
×
20 satisfy the condition (10). The simulation indicates that the numerical solution needs 26, 45, and 37 iterations in the RGAVS, the RGACS (η
=
0.74), and the traditional CGM, respectively. Therefore, the convergence speed of the RGAVS is also better than those of the RGACS and the CGM.


## 5. Conclusion

Using geometric structures of manifold consisting of all symmetric positive definite matrices, the Riemannian gradient algorithm is developed to calculate the numerical solution for the linear matrix equationQ
=
X
+
∑
i
=
1
m
‍
A
i
T
X
A
i. In this algorithm, the geodesic distance on the curved Riemannian manifold is taken as an objective function and the geodesic curve is treated as the convergence path. Also, the variable step sizes algorithm corresponding to the minimum value of the objective function is provided in order to improve the convergence speed in the constant step-size case. Finally, the convergence speed of the Riemannian gradient algorithm has been compared with the traditional conjugate gradient method by two simulation examples in 5-order and 20-order linear matrix equations, respectively. It is shown that the convergence speed of the Riemannian gradient algorithm with variable step sizes is faster than that of the conjugate gradient method. Although the Riemannian gradient algorithm is used here for the linear matrix equation in the real number field, it may be generalized to the complex case as well.



---

*Source: 507175-2014-01-06.xml*