# Generalized Preconditioned MHSS Method for a Class of Complex Symmetric Linear Systems

Authors: Cui-Xia Li; Yan-Jun Liang; Shi-Liang Wu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Based on the modified Hermitian and skew-Hermitian splitting (MHSS) and preconditioned MHSS (PMHSS) methods, a generalized preconditioned MHSS (GPMHSS) method for a class of complex symmetric linear systems is presented. Theoretical analysis gives an upper bound for the spectral radius of the iteration matrix. From a practical point of view, we have analyzed and implemented inexact GPMHSS (IGPMHSS) iteration, which employs Krylov subspace methods as its inner processes. Numerical experiments are reported to confirm the efficiency of the proposed methods.


---

## Body

## 1. Introduction

Consider an iterative solution of the system of linear equations as follows:(1)
A
x
=
b
,
A
∈
C
n
×
n
,
x
,
b
∈
C
n
,
where A
∈
C
n
×
n is a complex symmetric matrix of the following form:
(2)
A
=
W
+
i
T
,
with W
∈
R
n
×
n being symmetric positive definite and T
∈
R
n
×
n being symmetric positive semidefinite. Here and in the sequel, we use i
=
-
1 as the imaginary unit. One can readily verify that A is non-Hermitian, that is to say, the linear system (1) is a non-Hermitian linear system. System such as (1) is important and arises in a variety of scientific and engineering applications, including structural dynamics [1–4], diffuse optical tomography [5, 6], FFT-based solution of certain time-dependent PDEs [7], lattice quantum chromodynamics [8], molecular dynamics and fluid dynamics [9], quantum chemistry, and eddy current problem [10, 11]. One can see [12, 13] for more examples and additional references. In order to solve (1) more effectively, many efficient numerical algorithms have been proposed in the literature, see [14–20].Based on the specific structure of the coefficient matrixA, one can verify that the Hermitian and skew-Hermitian parts of the coefficient matrix A, respectively, are
(3)
H
=
1
2
(
A
+
A
*
)
=
W
,
S
=
1
2
(
A
-
A
*
)
=
i
T
.
Obviously, the above Hermitian and skew-Hermitian splitting of the coefficient matrix A is in line with the real and imaginary parts splitting of the coefficient matrix A. Based on the HSS method [21], Bai et al. [2] skillfully designed a modified HSS (MHSS) method to solve the complex symmetric linear system (1) which is described below.The MHSS Method. Let x
(
0
)
∈
C
n be an arbitrary initial guess. For k
=
0,1
,
2
,
…, until the sequence of iterates {
x
(
k
)
}
k
=
0
∞ converges, compute the next iterate x
(
k
+
1
) according to the following procedure:
(4)
(
α
I
+
W
)
x
(
k
+
1
/
2
)
=
(
α
I
-
i
T
)
x
(
k
)
+
b
,
(
α
I
+
T
)
x
(
k
+
1
)
=
(
α
I
+
i
W
)
x
(
k
+
1
/
2
)
-
i
b
,
where α is a given positive constant and I is an identity matrix.The potential advantage of the MHSS method over the HSS method [21] for solving the complex symmetric linear system (1) is that only two linear subsystems with coefficient matrices α
I
+
W and α
I
+
T, both being real and symmetric positive definite, need to be solved at each step. Therefore, in this case, these two linear subsystems can be solved efficiently using mostly real arithmetic either exactly by a sparse Cholesky factorization or inexactly by conjugated gradient scheme. That is to say, the MHSS method successfully avoids solving a shifted skew-Hermitian linear subsystem with coefficient matrix α
I
+
i
T.Theoretical analysis in [2] shows that the MHSS method converges unconditionally to the unique solution of the complex symmetric linear system (1) when W
∈
R
n
×
n is real symmetric positive definite and T
∈
R
n
×
n is real symmetric positive semidefinite. The corresponding optimum parameter α
=
λ
min
⁡
(
W
)
λ
max
⁡
(
W
) is obtained to minimize an upper bound on the spectral radius of the iteration matrix associated with (4).The MHSS method immediately attracted considerable attention and resulted in many papers devoted to various aspects of the new algorithms. For instance, preconditioned modified Hermitian and skew-Hermitian splitting (PMHSS) iteration in [3], lopsided preconditioned modified Hermitian and skew-Hermitian splitting (LPMHSS) iteration in [22], amongst others. On the other hand, the MHSS method was successfully extended to the solution of control problems in [23].In this paper, based on the splitting (3), we generalize the MHSS iterative scheme into a new approach, called generalized preconditioned MHSS (GPMHSS) iteration. By introducing two symmetric positive definite matrices P
1
∈
R
n
×
n and P
2
∈
R
n
×
n, the GPMHSS iterative scheme works as follows.The GPMHSS Method. Let x
(
0
)
∈
C
n be an arbitrary initial guess. For k
=
0,1
,
2
,
…, until the sequence of iterates {
x
(
k
)
}
k
=
0
∞ converges, compute the next iterate x
(
k
+
1
) according to the following procedure:
(5)
(
α
P
1
+
W
)
x
(
k
+
1
/
2
)
=
(
α
P
1
-
i
T
)
x
(
k
)
+
b
,
(
β
P
2
+
T
)
x
(
k
+
1
)
=
(
β
P
2
+
i
W
)
x
(
k
+
1
/
2
)
-
i
b
,
where α is a nonnegative constant and β is a positive constant.Note that the GPMHSS iteration method can cover many existing variants of the standard MHSS iteration. For instance, whenα
=
β and P
1
=
P
2
=
I, the GPMHSS iteration method is equivalent to the standard MHSS iteration method in [2]; when α
=
β and P
1
=
P
2, the GPMHSS iteration method is equivalent to the standard PMHSS iteration method in [3]; when α
=
0 and P
2
=
I, it leads to the LPMHSS iteration method in [22].Theoretical analysis gives an upper bound about the contraction factor of GPMHSS iteration method, which shows the relations among GPMHSS, MHSS, and other existing variants. From a practical point of view, we also discuss the inexact variants of the GPMHSS iteration method and their implementation. A number of numerical experiments are presented to illustrate the advantages of the GPMHSS and IGPMHSS methods.This paper is organized as follows. In Section2, we study the convergence properties of the GPMHSS iteration method. In Section 3, we discuss the implementation of GPMHSS iteration method and the corresponding inexact GPMHSS (IPGMHSS) iteration method. Numerical experiments are reported to confirm the efficiency of the proposed methods in Section 4. Finally, we end the paper with concluding remarks in Section 5.


## 2. Convergence Analysis for the GPMHSS Method

In this section, the convergence of the GPMHSS method is studied and an upper bound for the contraction factor of the GPMHSS method is derived. The GPMHSS iteration method can be generalized into a two-step splitting iteration framework. The following lemma is required to study the convergence properties of the GPMHSS method.The spectral radius of the matrixA is the nonnegative real number ρ
(
A
)
=
max
⁡
{
|
λ
|
:
λ
∈
σ
(
A
)
}, where σ
(
A
) denotes the spectrum of matrix A. In fact, we have the basic property on spectral radius of the product of two matrices that ρ
(
A
B
)
=
ρ
(
B
A
), which is used in the proof of the following theorem.Lemma 1 (see [21]).
LetA
∈
C
n
×
n, A
=
M
i
-
N
i
(
i
=
1,2
) be two splittings of A, and x
(
0
)
∈
C
n be a given initial vector. If {
x
(
k
)
} is a two-step iteration sequence defined by
(6)
M
1
x
(
k
+
1
/
2
)
=
N
1
x
(
k
)
+
b
,
M
2
x
(
k
+
1
)
=
N
2
x
(
k
+
1
/
2
)
+
b
,
k
=
0,1
,
…, then
(7)
x
(
k
+
1
)
=
M
2
-
1
N
2
M
1
-
1
N
1
x
(
k
)
+
M
2
-
1
(
I
+
N
2
M
1
-
1
)
b
,
0000000000000000000000000000000
k
=
0,1
,
…
.
Moreover, if the spectral radius ρ
(
M
2
-
1
N
2
M
1
-
1
N
1
)
<
1, then the iterative sequence {
x
(
k
)
} converges to the unique solution x
*
∈
C
n of the system (1) for all initial vectors x
(
0
)
∈
C
n.Applying this lemma to the GPMHSS method, we get the convergence property in the following theorem.Theorem 2.
LetP
1
∈
R
n
×
n and P
2
∈
R
n
×
n be two symmetric positive definite matrices. Let A
=
W
+
i
T
∈
C
n
×
n, with W
∈
R
n
×
n and T
∈
R
n
×
n symmetric positive definite and symmetric positive semidefinite, respectively, and let α be a nonnegative constant and let β be a positive constant. Then the iteration matrix M
α
,
β of GPMHSS method is
(8)
M
α
,
β
=
(
β
P
2
+
T
)
-
1
(
β
P
2
+
i
W
)
(
α
P
1
+
W
)
-
1
(
α
P
1
-
i
T
)
.
Denote γ
(
P
^
P
^
T
)
,
λ
(
W
^
), and μ
(
T
^
) to be the spectral sets of the matrices P
^
P
^
T, W
^, and T
^, respectively, where P
^
=
P
1
-
1
/
2
P
2
1
/
2, W
^
=
P
1
-
1
/
2
W
P
1
-
1
/
2, and T
^
=
P
2
-
1
/
2
T
P
2
-
1
/
2. Then,
(9)
ρ
(
M
α
,
β
)
≤
δ
=
κ
(
P
^
)
δ
1
δ
2
,
where κ
(
P
^
) denotes the spectral condition number of the matrix P
^, and
(10)
δ
1
=
(
β
max
⁡
γ
i
∈
γ
(
P
^
P
^
T
)
|
γ
i
-
1
|
+
|
β
-
α
|
)
max
⁡
λ
i
∈
λ
(
W
^
)
1
α
+
λ
i
+
max
⁡
λ
i
∈
λ
(
W
^
)
α
2
+
λ
i
2
α
+
λ
i
,
δ
2
=
(
α
max
⁡
γ
i
∈
λ
(
P
^
P
^
T
)
|
1
γ
i
-
1
|
+
|
β
-
α
|
)
max
⁡
μ
i
∈
μ
(
T
^
)
1
β
+
μ
i
+
1
.Proof.
LetM
1
=
α
P
1
+
W, N
1
=
α
P
1
-
i
T, M
2
=
β
P
2
+
T, and N
2
=
β
P
2
+
i
W. Obviously, α
P
1
+
W and β
P
2
+
T are nonsingular for any nonnegative constants α and positive constants β. So formula (8) is valid.
LetW
^
=
P
1
-
1
/
2
W
P
1
-
1
/
2, T
^
=
P
2
-
1
/
2
T
P
2
-
1
/
2, and P
^
=
P
1
-
1
/
2
P
2
1
/
2. Then,
(11)
α
P
1
+
W
=
P
1
1
/
2
(
α
I
+
W
^
)
P
1
1
/
2
,
β
P
2
+
i
W
=
P
1
1
/
2
(
β
P
^
P
^
T
+
i
W
^
)
P
1
1
/
2
,
β
P
2
+
T
=
P
2
1
/
2
(
β
I
+
T
^
)
P
2
1
/
2
,
α
P
1
-
i
T
=
P
2
1
/
2
(
α
P
^
-
1
P
^
-
T
-
i
T
^
)
P
2
1
/
2
.
Hence,
(12)
(
β
P
2
+
i
W
)
(
α
P
1
+
W
)
-
1
=
P
1
1
/
2
(
β
P
^
P
^
T
+
i
W
^
)
(
α
I
+
W
^
)
-
1
P
1
-
1
/
2
=
P
1
1
/
2
[
(
β
P
^
P
^
T
-
α
I
)
(
α
I
+
W
^
)
-
1
00000000000
+
(
α
I
+
i
W
^
)
(
α
I
+
W
^
)
-
1
]
P
1
-
1
/
2
=
P
1
1
/
2
[
β
(
P
^
P
^
T
-
I
)
(
α
I
+
W
^
)
-
1
+
(
β
-
α
)
(
α
I
+
W
^
)
-
1
00000000000
+
(
α
I
+
i
W
^
)
(
α
I
+
W
^
)
-
1
]
P
1
-
1
/
2
,
(
α
P
1
-
i
T
)
(
β
P
2
+
T
)
-
1
=
P
2
1
/
2
(
α
P
^
-
1
P
^
-
T
-
i
T
^
)
(
β
I
+
T
^
)
-
1
P
2
-
1
/
2
=
P
2
1
/
2
[
(
α
P
^
-
1
P
^
-
T
-
β
I
)
(
β
I
+
T
^
)
-
1
00000000000
+
(
β
I
-
i
T
^
)
(
β
I
+
T
^
)
-
1
]
P
2
-
1
/
2
=
P
2
1
/
2
[
α
(
P
^
-
1
P
^
-
T
-
I
)
(
β
I
+
T
^
)
-
1
+
(
α
-
β
)
00000000000
×
(
β
I
+
T
^
)
-
1
+
(
β
I
-
i
T
^
)
(
β
I
+
T
^
)
-
1
]
P
2
-
1
/
2
.
Further, we have
(13)
ρ
(
M
α
,
β
)
=
ρ
(
(
β
P
2
+
T
)
-
1
(
β
P
2
+
i
W
)
(
α
P
1
+
W
)
-
1
(
α
P
1
-
i
T
)
)
=
ρ
(
(
β
P
2
+
i
W
)
(
α
P
1
+
W
)
-
1
(
α
P
1
-
i
T
)
(
β
P
2
+
T
)
-
1
)
=
ρ
(
P
1
1
/
2
(
β
P
^
P
^
T
+
i
W
^
)
(
α
I
+
W
^
)
-
1
×
P
^
(
α
P
^
-
1
P
^
-
T
-
i
T
^
)
(
β
I
+
T
^
)
-
1
P
2
-
1
/
2
)
=
ρ
(
(
β
P
^
P
^
T
+
i
W
^
)
(
α
I
+
W
^
)
-
1
P
^
(
α
P
^
-
1
P
^
-
T
-
i
T
^
)
×
(
β
I
+
T
^
)
-
1
P
^
-
1
)
≤
∥
(
β
P
^
P
^
T
+
i
W
^
)
(
α
I
+
W
^
)
-
1
P
^
(
α
P
^
-
1
P
^
-
T
-
i
T
^
)
×
(
β
I
+
T
^
)
-
1
P
^
-
1
∥
2
≤
κ
(
P
^
)
[
β
∥
P
^
P
^
T
-
I
∥
2
∥
(
α
I
+
W
^
)
-
1
∥
2
+
|
β
-
α
|
×
∥
(
α
I
+
W
^
)
-
1
∥
2
+
∥
(
α
I
+
i
W
^
)
(
α
I
+
W
^
)
-
1
∥
2
]
·
[
α
∥
P
^
-
1
P
^
-
T
-
I
∥
2
∥
(
β
I
+
T
^
)
-
1
∥
2
+
|
α
-
β
|
×
∥
(
β
I
+
T
^
)
-
1
∥
2
+
∥
(
β
I
-
i
T
^
)
(
β
I
+
T
^
)
-
1
∥
2
]
.
Through simple calculations, we can get that
(14)
∥
(
α
I
+
i
W
^
)
(
α
I
+
W
^
)
-
1
∥
2
≤
max
⁡
λ
i
∈
λ
(
W
^
)
α
2
+
λ
i
2
α
+
λ
i
,
∥
(
β
I
-
i
T
^
)
(
β
I
+
T
^
)
-
1
∥
2
≤
max
⁡
μ
i
∈
μ
(
T
^
)
β
2
+
μ
i
2
β
+
μ
i
≤
1
,
∥
P
^
P
^
T
-
I
∥
2
≤
max
⁡
γ
i
∈
γ
(
P
^
P
^
T
)
|
γ
i
-
1
|
,
∥
P
^
-
1
P
^
-
T
-
I
∥
2
≤
max
⁡
γ
i
∈
γ
(
P
^
P
^
T
)
|
1
γ
i
-
1
|
.
The above inequalities give the upper bound for ρ
(
M
α
,
β
) in (9).Some remarks on Theorem2 are given below.(i) WhenP
1
=
P
2, GPMHSS reduces to a two-parameter PMHSS method. In this case, P
^
=
P
1
-
1
/
2
P
2
1
/
2
=
I and κ
(
P
^
)
=
1. In the meantime, the upper bound in (9) results in
(15)
δ
=
(
max
⁡
λ
i
∈
λ
(
W
^
)
|
β
-
α
|
α
+
λ
i
+
max
⁡
λ
i
∈
λ
(
W
^
)
α
2
+
λ
i
2
α
+
λ
i
)
×
(
max
⁡
μ
i
∈
μ
(
T
^
)
|
β
-
α
|
β
+
μ
i
+
1
)
,
which not only gives the upper bound for the PMHSS method and complements the theoretical results of paper [24], but also includes the special case in [2] as follows:
(16)
δ
=
max
⁡
λ
i
∈
λ
(
W
^
)
α
2
+
λ
i
2
α
+
λ
i
,
when α
=
β.(ii) Based on (9), it is obvious that the convergence rate not only depends on the choices of two parameters α and β, but also depends on the choices of two auxiliary matrices P
1 and P
2. Obviously, the efficiency of GPMHSS is best when one can obtain two iteration parameters and auxiliary matrices to minimize an upper bound on the spectral radius of the iteration matrix. That is to say, it is necessary to discuss the choices of α, β, P
1, and P
2. In fact, under certain conditions, in theory, some optimal iteration parameters α and β can be derived. For example, if P
1
=
0, then the GPMHSS method reduces to the LPMHSS method and the corresponding optimal parameter is λ
min
⁡
2
/
μ
max
⁡, where λ
min
⁡ and μ
max
⁡, respectively, denote the smallest eigenvalue of P
2
-
1
W and the largest eigenvalue of P
2
-
1
T [22]. If P
1
=
P
2 and α
=
β, then the GPMHSS method reduces to the PMHSS method and the corresponding optimal parameter is λ
min
⁡
λ
max
⁡, where λ
min
⁡ and λ
max
⁡, respectively, denote the smallest and largest eigenvalues of P
2
-
1
W [3]. Whereas, if P
1
=
P
2 and α
≠
β, in theory, one cannot find the optimal values of the iteration parameters α and β in general. With respect to this point, one can see [24] for more details. Based on this fact, in theory, a conclusion may be obtained, that is, one cannot find the optimal values of the iteration parameters α and β for GPMHSS in general. But even so, the optimal values of the iteration parameters α and β for GPMHSS may be obtained experimentally (to see the fourth section). In practice, to further improve the efficiency of the GPMHSS method, it is desirable to determine or find a good estimate of the optimal parameter that minimizes the convergence factor. In fact, to find the actual optimal estimates α and β for the GPMHSS method is a hard task because its solution strongly depends on the particular structures and properties of the coefficient matrix A, as well as the splittings matrices W and T, and needs further in-depth study, both from theory and computation point of view. When the optimal parameter cannot be derived in theory, the value of parameter α
=
β with P
1
=
P
2 is selected by the statement on the choice of the iteration parameter in [25], that is to say, experience suggests that in most applications and for an appropriate scaling of the problem, a “small” value of parameter (usually between 0.01 and 0.5) may give good results. Whereas choosing a parameter so as to minimize the spectral radius of the iteration matrix is not necessarily the best choice when the algorithm is used as a preconditioner for a Krylov subspace method. Remarkably, it can be shown that, for certain problems, the alternating iteration results in an h-independent preconditioner for GMRES when parameter is chosen sufficiently small, correspond to a spectral radius very close to 1. In fact, if we define the optimal value of parameter as the one that minimizes the total amount of work needed to compute an approximate solution, this will not necessarily be the same as the parameter that minimizes the number of (outer) iterations. Overall, the analytic determination of such an optimal value for parameter appears to be daunting. With respect to the choices of two auxiliary matrices P
1 and P
2, P
1 and P
2 may have a better degree of freedom. For instance, if W
+
T is positive definite, then P
1
=
T and P
2
=
W with α
>
0 would be a natural choice for the preconditioners.(iii) WhenP
1
≠
P
2 and α
=
β, the upper bound in (9) results in
(17)
δ
=
(
max
⁡
γ
i
∈
γ
(
P
^
P
^
T
)
|
γ
i
-
1
|
max
⁡
λ
i
∈
λ
(
W
^
)
α
α
+
λ
i
+
max
⁡
λ
i
∈
λ
(
W
^
)
α
2
+
λ
i
2
α
+
λ
i
)
·
(
max
⁡
γ
i
∈
λ
(
P
^
P
^
T
)
|
1
γ
i
-
1
|
max
⁡
μ
i
∈
μ
(
T
^
)
α
α
+
μ
i
+
1
)
.The approach to minimize the upper bound is very important in theoretical viewpoint. However, it is not practical since the corresponding spectral radius of the iteration matrixM
α
,
β is not optimal. How to choose the suitable preconditioners and parameters for practical problem is still a great challenge.


## 3. The IGPMHSS Iteration

In the GPMHSS method, it is required to solve two systems of linear equations whose coefficient matrices areα
P
1
+
W and β
P
2
+
T, respectively. However, this may be very costly and impractical in actual implementations. To overcome this disadvantage and improve the computational efficiency of the GPMHSS iteration method, we propose to solve the two subproblems iteratively [21, 26], which leads to the inexact GPMHSS (IGPMHSS) iterative scheme. Its convergence can be shown in a similar way to that of the IHSS iteration method, using Theorem 3.1 of [21]. Since α
P
1
+
W and β
P
2
+
T are symmetric positive definite, there can employ some Krylov subspace methods (such as CG) to gain its solution easily. Of course, if good preconditioners for matrices α
P
1
+
W and β
P
2
+
T are available, we can use the preconditioned conjugate gradient (PCG) method instead of CG for the two inner systems, this yields a better performance of IGPMHSS method. If either α
P
1
+
W or β
P
2
+
T (or both) is Toeplitz, we can use fast algorithms for solution of the corresponding subsystems [27]. Here, just like the IHSS iteration, the IGPMHSS iterative scheme is presented in Algorithm 1.Algorithm 1
k
=
0;
while (not convergent)
r
(
k
)
=
b
-
A
x
(
k
);
approximately solve(
α
P
1
+
W
)
z
(
k
)
=
r
(
k
) by employing CG method,
such that the residualp
(
k
)
=
r
(
k
)
-
(
α
P
1
+
W
)
z
(
k
) of the iteration satisfies
∥
p
(
k
)
∥
2
≤
η
k
∥
r
(
k
)
∥
2;
x
(
k
+
1
/
2
)
=
x
(
k
)
+
z
(
k
);
r
(
k
+
1
/
2
)
=
b
+
i
A
x
(
k
+
1
/
2
);
approximately solve(
β
P
2
+
T
)
z
(
k
+
1
/
2
)
=
r
(
k
+
1
/
2
) by employing CG method,
such that the residualq
(
k
+
1
/
2
)
=
r
(
k
+
1
/
2
)
-
(
β
P
2
+
T
)
z
(
k
+
1
/
2
) of the
iteration satisfies∥
q
(
k
+
1
/
2
)
∥
2
≤
τ
k
∥
r
(
k
+
1
/
2
)
∥
2;
x
(
k
+
1
)
=
x
(
k
+
1
/
2
)
+
z
(
k
+
1
/
2
);
k
=
k
+
1;
endRemark 3.
It is not difficult to find the fact that if the inner systems can be solved exactly withη
k
=
0 and τ
k
=
0, then the IGPMHSS iteration reduces to the GPMHSS iteration. In fact, to guarantee the convergence of the IGPMHSS iteration, it is not necessary for η
k and τ
k to go to zero when k increases.
To derive the convergence properties for IGPMHSS iteration, the following lemma is required, which is presented by Bai et al. [21]. Here ∥
|
x
|
∥
M
=
∥
M
x
∥ for any x
∈
C
n, which immediately induces the matrix norm ∥
|
X
|
∥
M
=
∥
M
X
M
-
1
∥ for any X
∈
C
n
×
n.Lemma 4.
LetA
∈
C
n
×
n, A
=
M
i
-
N
i
(
i
=
1,2
) be two splittings of A. If {
x
(
k
)
} is an iteration sequence defined as follows:
(18)
x
(
k
+
1
/
2
)
=
x
(
k
)
+
z
(
k
)
w
i
t
h
M
1
z
(
k
)
=
r
(
k
)
+
p
(
k
)
,
satisfying ∥
p
(
k
)
∥
2
≤
η
k
∥
r
(
k
)
∥
2, where r
(
k
)
=
b
-
A
x
(
k
), and
(19)
00000000
x
(
k
+
1
)
=
x
(
k
+
1
/
2
)
+
z
(
k
+
1
/
2
)
w
i
t
h
M
2
z
(
k
+
1
/
2
)
=
r
(
k
+
1
/
2
)
+
q
(
k
+
1
/
2
)
,
satisfying ∥
q
(
k
+
1
/
2
)
∥
2
≤
τ
k
∥
r
(
k
+
1
/
2
)
∥
2, where r
(
k
+
1
/
2
)
=
b
-
A
x
(
k
+
1
/
2
), then {
x
(
k
)
} is of the form
(20)
x
(
k
+
1
)
=
M
2
-
1
N
2
M
1
-
1
N
1
x
(
k
)
+
M
2
-
1
(
I
+
N
2
M
1
-
1
)
b
+
M
2
-
1
(
N
2
M
1
-
1
p
(
k
)
+
q
(
k
+
1
/
2
)
)
.
Moreover, if x
*
∈
C
n is the exact solution of the system of linear equations (1), then we have
(21)
∥
|
x
(
k
+
1
)
-
x
*
|
∥
M
2
≤
(
ζ
+
μ
θ
η
k
+
θ
(
ρ
+
θ
ν
η
k
)
τ
k
)
∥
|
x
(
k
)
-
x
*
|
∥
M
2
,
0000000000000000000000000000000
k
=
1,2
,
…
,
where
(22)
ζ
=
∥
N
2
M
1
-
1
N
1
M
2
-
1
∥
2
,
ρ
=
∥
M
2
M
1
-
1
N
1
M
2
-
1
∥
2
,
μ
=
∥
N
2
M
1
-
1
∥
2
,
θ
=
∥
A
M
2
-
1
∥
2
,
ν
=
∥
M
2
M
1
-
1
∥
2
.
In particular, if
(23)
ζ
+
μ
θ
η
max
⁡
+
θ
(
ρ
+
θ
ν
η
max
⁡
)
τ
max
⁡
<
1
,
then the iteration sequence {
x
(
k
)
} converges to x
*
∈
C
n, where η
max
⁡
=
max
⁡
k
{
η
k
} and τ
max
⁡
=
max
⁡
k
{
τ
k
}.Applying Lemma4, the following theorem gives a convergence analysis of the IGPMHSS iteration.Let(24)
δ
α
,
β
=
∥
(
α
P
1
-
i
T
)
(
α
P
1
+
W
)
-
1
(
β
P
2
+
i
W
)
(
β
P
2
+
T
)
-
1
∥
2
.
Then we have the following result.Theorem 5.
LetP
1
∈
R
n
×
n and P
2
∈
R
n
×
n be two symmetric positive definite matrices. Let A
=
W
+
i
T
∈
C
n
×
n, with W
∈
R
n
×
n and T
∈
R
n
×
n symmetric positive definite and symmetric positive semidefinite, respectively, and let α be a nonnegative constant and let β be a positive constant. If {
x
(
k
)
} is an iteration sequence generated by the IGPMHSS iteration method (Algorithm 1) and if x
*
∈
C
n is the exact solution of the system of linear equations (1), then it holds that
(25)
∥
|
x
(
k
+
1
)
-
x
*
|
∥
≤
(
δ
α
,
β
+
μ
θ
η
k
+
θ
(
ρ
+
θ
υ
η
k
)
τ
k
)
∥
|
x
(
k
)
-
x
*
|
∥
,
00000000000000000000000000000000
k
=
1,2
,
…
,
where
(26)
ρ
=
∥
(
β
P
2
+
T
)
(
α
P
1
+
W
)
-
1
(
α
P
1
-
i
T
)
(
β
P
2
+
T
)
-
1
∥
2
,
θ
=
∥
A
(
β
P
2
+
T
)
-
1
∥
2
,
μ
=
∥
(
β
P
2
+
i
W
)
(
α
P
1
+
W
)
-
1
∥
2
,
υ
=
∥
(
β
P
2
+
T
)
(
α
P
1
+
W
)
-
1
∥
2
.
In particular, if
(27)
δ
α
,
β
+
μ
θ
η
max
⁡
+
θ
(
ρ
+
θ
ν
η
max
⁡
)
τ
max
⁡
<
1
,
then the iteration sequence {
x
(
k
)
} converges to x
*, where η
max
⁡
=
max
⁡
k
{
η
k
} and τ
max
⁡
=
max
⁡
k
{
τ
k
}.In fact, replacingM
i and N
i
(
i
=
1,2
) in Lemma 4 with M
1
=
α
P
1
+
W, N
1
=
α
P
1
-
i
T, M
2
=
β
P
2
+
T, and N
2
=
β
P
2
+
i
W, we straightforwardly obtain the proof of Theorem 5.Theorem5 shows that the choices of the tolerances η
k and τ
k are to make the IGPMHSS convergent. In fact, as mentioned previously, to get the convergence of the IGPMHSS iteration, the tolerances η
k and τ
k are not required to approach zero as k increases. However, the optimal tolerances η
k and τ
k are not easy to analyze.


## 4. Numerical Experiments

In this section, we give some numerical experiments to demonstrate the performance of the GPMHSS and IGPMHSS methods for solving the linear system (1)-(2). Numerical comparisons with the MHSS and PMHSS methods are also presented to show the advantage of the GPMHSS method.In our implementations, the initial guess is chosen to bex
(
0
)
=
0 and the stopping criteria for outer iterations (when MHSS, PMHSS, and GPMHSS methods are used as solvers) are
(28)
∥
b
-
A
x
(
k
)
∥
2
∥
b
∥
2
≤
1
0
-
6
.The preconditionerP used in PMHSS method is chosen to be P
=
W. For the sake of comparing, the corresponding preconditioners P
1 and P
2 used in GPMHSS method are chosen to be P
1
=
P
2
=
W. Similarly, if the preconditioner P used in PMHSS method is chosen to be P
=
T, then corresponding preconditioners P
1 and P
2 used in GPMHSS method are chosen to be P
1
=
P
2
=
T. Since the numerical results in [2, 3] show that the PMHSS method outperforms the MHSS and HSS methods when they are employed as preconditioners for the GMRES method or its restarted variants [28], we just examine the efficiency of the GPMHSS method as a solver for solving complex symmetric linear system (1)-(2) by comparing the iteration numbers (denoted as IT) and CPU times (in seconds, denoted as CPU) of the GPMHSS method with those of the MHSS and PMHSS methods.Example 6 (see [2, 3]).
The complex symmetric linear system (1)-(2) is of the following form:
(29)
A
x
≡
(
W
+
i
T
)
x
=
b
,
with
(30)
T
=
I
⊗
V
+
V
⊗
I
,
(31)
W
=
10
(
I
⊗
V
c
+
V
c
⊗
I
)
+
9
(
e
1
e
m
T
+
e
m
e
1
T
)
⊗
I
,
where V
= tridag (
-
1,2
,
-
1
)
∈
R
m
×
m, V
c
=
V
-
e
1
e
m
T
-
e
m
e
1
T
∈
R
m
×
m, and e
1 and e
m are the first and the last unit vectors in R
m, respectively. We take the right-hand side vector b to be b
=
(
1
+
i
)
A
1, with1 being the vector of all entries equal to 1. Here T and W correspond to the five-point centered difference matrices approximating the negative Laplacian operator with homogeneous Dirichlet boundary conditions and periodic boundary conditions, respectively, on a uniform mesh in the unit square [
0,1
]
×
[
0,1
] with the mesh-size h
=
1
/
(
m
+
1
).As is known, the spectral radius of the iteration matrix may be decisive for the convergence of the iteration method. The spectral radius corresponding to the iteration method is necessary to consider. The comparisons of the spectral radius of the three different iteration matrices derived by MHSS, PMHSS, and GPMHSS methods with different mesh-sizes are performed in Tables1 and 2. In Tables 1 and 2, we use the optimal values of the parameters α and β, denoted by α
M
* for MHSS method, α
P
* for PMHSS method, and α
G
*, β
G
* for GPMHSS method. These parameters are obtained experimentally with the least spectral radius for the iteration matrices of the three methods.Table 1
The experimentally optimal parameters and the spectral radii for the iteration matrices of MHSS, PMHSS, and GPMHSS withP
=
P
1
=
P
2
=
W.
m
8
×
8
16
×
16
24
×
24
32
×
32
MHSS
α
M
*
3.7
2.1
1.5
1.2
ρ
(
M
α
M
*
)
0.7203
0.7989
0.8407
0.8662
PMHSS
α
P
*
0.8
0.8
0.8
0.8
ρ
(
M
α
P
*
)
0.6638
0.6672
0.6684
0.6690
GPMHSS
α
G
*
0.8
0.8
0.8
0.8
β
G
*
3
2
1.6
1.4
ρ
(
M
α
G
*
,
β
G
*
)
0.4618
0.4851
0.5096
0.5321Table 2
The experimentally optimal parameters and the spectral radii for the iteration matrices of MHSS, PMHSS, and GPMHSS withP
=
P
1
=
P
2
=
T.
m
8
×
8
16
×
16
24
×
24
32
×
32
MHSS
α
M
*
3.7
2.1
1.5
1.2
ρ
(
M
α
M
*
)
0.7203
0.7989
0.8407
0.8662
PMHSS
α
P
*
1.2
1.2
1.2
1.2
ρ
(
M
α
P
*
)
0.6636
0.6671
0.6684
0.6690
GPMHSS
α
G
*
1.2
1.2
1.2
1.2
β
G
*
6.4
2.2
1.9
1.8
ρ
(
M
α
G
*
,
β
G
*
)
0.3389
0.4814
0.5184
0.5461From Tables1 and 2, one can see that with the mesh-size m increasing, the trend of the experimentally optimal parameter of the MHSS method decreases. Whereas, with the mesh-size m increasing, the experimentally optimal parameter of the PMHSS method may be changeless. Fixing the parameter α in GPMHSS method, the trend of the experimentally optimal parameter β of the MHSS method decreases. No matter what P
=
P
1
=
P
2
=
W or P
=
P
1
=
P
2
=
T, we observe that the optimal spectral radius for the iteration matrices of the three methods grows with problem size, and the optimal spectral radius of GPMHSS method is still smaller than those of MHSS and PMHSS methods. In this case, the GPMHSS method may outperform the MHSS and PMHSS methods. To this end, we need to examine the efficiencies of the MHSS, PHSS, and GPMHSS methods for solving the systems of linear equations A
x
=
b, where A is described above.In Tables3 and 4, we list the iteration numbers and computational times for MHSS, PMHSS, and GPMHSS iteration methods by using the optimal parameters in Tables 1 and 2. In Tables 3 and 4, “RES” denotes the relative residual error.Table 3
RES, CPUs, and IT for MHSS, PMHSS, and GPMHSS withP
=
P
1
=
P
2
=
W.
m
8
×
8
16
×
16
24
×
24
32
×
32
MHSS
RES
9.733
e
-
7
9.388
e
-
7
9.337
e
-
7
9.369
e
-
7
CPU(s)
0.015
0.828
7.422
42.328
IT
46
75
99
120
PMHSS
RES
6.585
e
-
7
6.530
e
-
7
6.473
e
-
7
6.429
e
-
7
CPU(s)
0.015
0.062
0.141
0.281
IT
31
31
31
31
Speed-up
1
13.35
52.64
150.63
GPMHSS
RES
7.884
e
-
7
8.197
e
-
7
8.210
e
-
7
8.032
e
-
7
CPU(s)
0.015
0.032
0.094
0.219
IT
18
19
20
21
Speed-up
1
25.88
78.96
193.28Table 4
RES, CPUs, and IT for MHSS, PMHSS, and GPMHSS withP
=
P
1
=
P
2
=
T.
m
8
×
8
16
×
16
24
×
24
32
×
32
MHSS
RES
9.733
e
-
7
9.388
e
-
7
9.337
e
-
7
9.369
e
-
7
CPU(s)
0.015
0.828
7.422
42.328
IT
46
75
99
120
PMHSS
RES
7.003
e
-
7
6.959
e
-
7
6.901
e
-
7
6.857
e
-
7
CPU(s)
0.015
0.047
0.14
0.25
IT
31
31
31
31
Speed-up
1
17.62
53.01
169.31
GPMHSS
RES
7.340
e
-
7
5.868
e
-
7
6.583
e
-
7
6.623
e
-
7
CPU(s)
0.015
0.031
0.078
0.172
IT
14
18
20
21
Speed-up
1
26.71
95.15
246.09From Tables3 and 4, we see that GPMHSS iteration method is the best among three methods in terms of the iteration numbers and computational time, and the PMHSS scheme requires fewer iteration numbers than the MHSS scheme. For the MHSS and GPMHSS methods, the number of iterations grows with the problem size. For the PMHSS method, its iteration numbers are relatively stable. That is to say, the PMHSS method does not have any growth in iterations numbers by increasing grid dimension. The presented results in Tables 3 and 4 show that in all cases GPMHSS is superior to another two methods in terms of the CPU time. That is to say, compared with the MHSS and PMHSS methods, the GPMHSS method may be given priority under certain conditions. In fact, the speed-up of GDMHSS/PMHSS with respect to MHSS is quite noticeable, where we define it by
(32)
speed-up
=
CPU
of
HSS
method
CPU
of
GPMHSS
(or
PMHSS)
method
.
Obviously, the efficiency of GPMHSS is superior to that of MHSS and PMHSS.In particular, here we test the efficiency of the GPMHSS method whenP
1
=
T and P
2
=
W. In this case, some numerical results are obtained in Table 5. In Table 5, it is not difficult to find that the spectral radius grows with problem size. Simultaneously, the number of iterations grows with the problem size. Numerical results in Table 5 show that, under certain conditions, the GPMHSS method is feasible and efficient, compared with the MHSS and PMHSS methods.Table 5
Numerical results for GPMHSS withP
1
=
T and P
2
=
W.
m
8
×
8
16
×
16
24
×
24
32
×
32
α
G
*
1.3
1.3
1.3
1.3
β
G
*
1.7
1.4
1.4
2.2
ρ
(
M
α
G
*
,
β
G
*
)
0.2644
0.2988
0.3044
0.3191
RES
4.7850
e
-
7
5.349
e
-
7
9.151
e
-
7
8.141
e
-
7
CPU(s)
0.015
0.031
0.078
0.172
IT
15
18
18
20
Speed-up
1
26.71
95.15
246.09As already noted, in the two half steps of the GPMHSS iteration, it is necessary to solve two systems of linear equations, whose coefficient matrices areα
P
1
+
W and β
P
2
+
T, respectively. This can be very costly and impractical in actual implementations. We use the IGPMHSS method to solve the systems of linear equations (1)-(2) in the actual implementations. That is, it is necessary to solve two systems of linear equations with α
P
1
+
W and β
P
2
+
T by using the IGPMHSS iteration. It is easy to know that α
P
1
+
W and β
P
2
+
T are symmetric positive definite. So, solving the above two subsystems, there can employ the CG method.In our computations, the inner CG iteration is terminated if the current residual of the inner iterations satisfies(33)
∥
p
(
j
)
∥
2
∥
r
(
k
)
∥
2
-
≤
0.1
τ
(
k
)
,
∥
q
(
j
)
∥
2
∥
r
(
k
)
∥
2
-
≤
0.1
τ
(
k
)
,
(cf. Algorithm 1) where p
(
j
) and q
(
j
) are, respectively, the residuals of the jth inner CG for α
P
1
+
W and β
P
2
+
T. r
(
k
) is the kth outer IGPMHSS iteration, and τ is a tolerance.Some results are listed in Tables6 and 7, which are the numbers of outer IGPMHSS iteration (it.s), the average numbers (avg1) of inner CG iteration for α
P
1
+
W, and the average numbers (avg2) of CG iteration for β
P
2
+
T.Table 6
Convergence results for IGPMHSS withP
1
=
P
2
=
W.
(
α
G
,
β
G
)
m
τ
=
0.95
τ
=
0.9
τ
=
0.85
it.s
avg1
avg2
it.s
avg1
avg2
it.s
avg1
avg2
(0.8, 2)
16
55
21.5
22.5
69
26.3
27.1
82
30.4
30.8
(0.8, 1.6)
24
82
32.1
33.1
105
39.4
40.1
125
45.8
46.1
(0.8, 1.4)
32
108
43.4
42.7
139
52.8
52.4
167
61.3
60.8
(0.8, 1.2)
64
160
95.4
54.2
190
96.4
70.2
198
96.9
80.8
(0.8, 1.2)
128
185
95.1
91.7
178
95.2
92.3
184
95.3
91.9Table 7
Convergence results for IGPMHSS withP
1
=
P
2
=
T.
(
α
G
,
β
G
)
m
τ
=
0.95
τ
=
0.9
τ
=
0.85
it.s
avg1
avg2
it.s
avg1
avg2
it.s
avg1
avg2
(1.2, 2.2)
16
49
21.7
15.6
60
26.7
19.5
71
30.8
22.2
(1.2, 1.9)
24
70
29.7
21.8
91
38.1
28.1
109
45.0
33.0
(1.2, 1.8)
32
90
37.1
27.6
119
48.7
36.8
144
58.4
43.4
(1.2, 1.6)
64
157
92.5
54.2
200
93.2
75.1
200
93.1
88.6
(1.2, 1.6)
128
192
92.3
89.1
192
92.3
90.4
192
92.4
91.6In our numerical computations, it is easy to find the fact that the choice ofτ is important to the convergence rate of the IGPMHSS method. According to Tables 6 and 7, the iteration numbers of the IGPMHSS method generally increases when τ decreases. This increase can probably be eliminated using a suitable preconditioner. It is noted that no preconditioning is used for these inner iterations in our numerical computations.


## 5. Conclusion

In this paper, we have generalized the MHSS method into the GPMHSS method for a class of complex symmetric linear systems. Theoretical analysis shows that for any initial guess the GPMHSS method converges to the unique solution of the linear system for a wide range of the parameters. Then, an inexact version of GPMHSS (IGPMHSS) has been presented and implemented for saving the computational cost. Numerical experiments show that the GPMHSS and IGPMHSS methods are efficient and competitive.



---

*Source: 206821-2014-05-26.xml*