# Preconditioning Complex Symmetric Linear Systems

Authors: Enrico Bertolazzi; Marco Frego
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A new preconditioner for symmetric complex linear systems based on Hermitian and skew-Hermitian splitting (HSS) for complex symmetric linear systems is herein presented. It applies to conjugate orthogonal conjugate gradient (COCG) or conjugate orthogonal conjugate residual (COCR) iterative solvers and does not require any estimation of the spectrum of the coefficient matrix. An upper bound of the condition number of the preconditioned linear system is provided. To reduce the computational cost the preconditioner is approximated with an inexact variant based on incomplete Cholesky decomposition or on orthogonal polynomials. Numerical results show that the present preconditioner and its inexact variant are efficient and robust solvers for this class of linear systems. A stability analysis of the inexact polynomial version completes the description of the preconditioner.


---

## Body

## 1. Introduction

Focus of this paper is the solution of the complex linear system given byA
x
=
b, where the symmetric complex matrix A has the property that can be written as A
=
B
+
i
C with B, C two real symmetric positive semidefinite matrices (semi-SPD) and B
+
C a symmetric positive definite (SPD) matrix. This kind of linear system arises, for example, in the discretization of problems in computational electrodynamics [1] or time-dependent Schrödinger equations, or in conductivity problems [2, 3].IfA is Hermitian, a straightforward extension of the conjugate gradients (CG) algorithm can be used [4]. Unfortunately, the CG method cannot be directly employed when A is only complex symmetric; thus, some specialized iterative methods must be adopted. An effective one is the HSS with its variants (MHSS), which need, at each iteration, the solution of two real linear systems. Other standard procedures to solve this problem are given by numerical iterative methods based on Krylov spaces and designed for complex symmetric linear systems: COCG [1], COCR [5], CSYM [6], and CMRH [7]. Some iterative methods for non-SPD linear systems like BiCGSTAB [8], BiCGSTAB(l) [9, 10], GMRES [11], and QMR [12] can be adapted for complex symmetric matrices [4, 13, 14].In this work a preconditioned version of COCG and COCR is proposed: the aim is the solution of complex linear systems, where the preconditioner is a single preconditioned MHSS iteration (PMHSS). The PMHSS step may be too costly or impracticable, thus cheaper approximations must be used, in this work some approximations are considered, and a polynomial preconditioner as a possible approximation of the PMHSS step is presented.Methods based on Hermitian and skew-Hermitian splitting (HSS) [15–18] can be used as standalone solvers or combined (as preconditioner) together with CG like algorithms. The speed of convergence of CG like iterative schemes depends on the condition number of the matrix A; thus, preconditioning is a standard way to improve convergence [19, 20]. Incomplete LU is a standard and accepted way to precondition linear systems. Despite its popularity, incomplete LU is potentially unstable, is difficult to parallelize, and lacks algorithmic scalability. Nevertheless, when incomplete LU is feasible and the preconditioned linear system is well conditioned, the resulting algorithm is generally the best performing.In this work we focus on large problems, where incomplete LU preconditioning is too costly or not feasible. In this case, iterative methods like SSOR are used as preconditioners, but a better performance is obtained using HSS iterative methods, which allow reducing the condition number effectively. However, HSS iterative methods need the solution of two SPD real systems at each step. Standard preconditioned CG methods can be used at each iteration [18] which can again be preconditioned with an incomplete Cholesky factorization, although for very large problems, the incomplete Cholesky factorization may be not convenient or not feasible. As an alternative, in the present paper a polynomial preconditioner is proposed, which allows solving the linear system for large matrices. Polynomial preconditioners do not have a good reputation as preconditioners [19, 20] and research on this subject was dropped in the late 80s. In fact, polynomial preconditioners based on Chebyshev polynomials need accurate estimates of minimum and maximum eigenvalues, while least squares polynomials were not computed using stable recurrences, limiting the degree of available stable polynomials [21–25]. However, in the last years, polynomial preconditioner received attention again after the works of Lu et al. [26]. Notwithstanding anything contained above, here we propose as a preconditioner the use of a polynomial approximation of a modified HSS step. A specialization for Chebyshev and Jacobi orthogonal polynomials is discussed and the corresponding polynomial preconditioner is evaluated using a stable recurrence which permits using very high degree polynomials.The paper has this structure. Section1.1 describes the problem and gives a brief summary of existing methods for the resolution with the fundamental results and variants that lead to the present method; in particular, the HSS is considered. Section 2 shows how to use one step of the MHSS method as a preconditioner and gives a bound on the conditioning number of the MHSS iteration matrix. Section 3 explains the iterative solution method with the strategy to adopt when Cholesky factorization is possible or not. Section 4 presents a scale transformation of the system in order to move the eigenvalues to the range (
0,1
]. Section 5 describes the polynomial preconditioner based first on least squares and then in terms of orthogonal polynomials and furnishes a stable recurrence for the computation of polynomial preconditioners for high degrees. The specialization for Chebyshev and Jacobi orthogonal polynomials is presented. Section 6 studies the numerical stability of this process and Section 7 shows some numerical results; Section 8 concludes the paper.


### 1.1. The MHSS Iterative Solver

The complexN
×
N linear system A
x
=
b, where A is complex symmetric, is solved via an iterative method based on a splitting algorithm (HSS). The preconditioner requires to solve (each time it is applied) two real symmetric and positive definite (SPD) linear systems.The HSS scheme can be summarized in this manner: rewrite the vectorx of the unknowns as the sum of real and imaginary part, x
=
y
+
i
z, and accordingly, the right hand side b
=
c
+
i
d, then the system A
x
=
b is rewritten as(1)
B
+
i
C
y
+
i
z
=
c
+
i
d
,so that the two steps of the modified HSS method proposed in [15] result in(2)
V
+
B
x
k
+
1
/
2
=
V
-
i
C
x
k
+
b
,
W
+
C
x
(
k
+
1
)
=
W
+
i
B
x
(
k
+
1
/
2
)
-
i
b
,for suitable matrices V and W. The previous procedure can be rewritten as a single step of a splitting based scheme P
x
k
+
1
=
Q
x
k
+
b by posing(3)
P
=
V
+
B
W
-
i
V
-
1
W
+
C
,
Q
=
V
+
B
W
-
i
V
-
1
W
+
i
B
V
+
B
-
1
V
-
i
C
.This iterative method converges if the iteration matrix P
-
1
Q, for example,(4)
P
-
1
Q
=
W
+
C
-
1
W
+
i
B
V
+
B
-
1
V
-
i
C
,has spectral radius strictly less than one. It is well known that the choice V
=
W
=
α
I, for a given positive constant α, yields the standard HSS method [15–17] for which an estimate of the spectral radius is given by(5)
ϱ
P
-
1
Q
≤
max
⁡
j
=
1,2
,
…
,
N
α
2
+
λ
j
B
2
α
+
λ
j
B
,where λ
j
(
B
) are the eigenvalues of the SPD matrix B. The optimal value for α can also be computed [15–17] and is(6)
α
opt
=
arg
min
⁡
α
max
⁡
j
=
1,2
,
…
,
n
α
2
+
λ
j
B
2
α
+
λ
j
B
α
opt
=
λ
min
⁡
B
λ
max
⁡
B
.From the previous formulas, it is clear that those computations rely on the knowledge (or estimate) of the minimum and maximum eigenvalue of the matrix B, which is, in general, a hard problem.Another possible choice forV and W is V
=
α
B and W
=
β
B which yields, when α
=
β, a variant of the MHSS method by [15–17]. In the next lemma, an upper bound of the spectral radius of the iteration matrix is given.Lemma 1.
LetV
=
α
B and W
=
β
B in (4) with B a SPD matrix and C a semi-SPD matrix, α
,
β
>
0; then the spectral radius of P
-
1
Q satisfies the upper bound(7)
ϱ
P
-
1
Q
≤
U
α
,
β
,
U
α
,
β
=
1
+
β
2
1
+
α
max
⁡
1
,
α
βand the minimum value of the upper bound U
(
α
,
β
) is attained when α
=
β
=
1 where U
(
1,1
)
=
2
/
2
≈
0.707.Proof.
SubstitutingV
=
α
B and W
=
β
B in (4) yields(8)
P
-
1
Q
=
β
B
+
C
-
1
β
B
+
i
B
α
B
+
B
-
1
α
B
-
i
C
P
-
1
Q
=
β
+
i
1
+
α
β
B
+
C
-
1
α
B
-
i
C
.If λ is an eigenvalue of matrix P
-
1
Q, it satisfies(9)
0
=
det
⁡
P
-
1
Q
-
λ
I
0
⇓
0
=
det
⁡
β
+
i
1
+
α
α
B
-
i
C
-
λ
β
B
+
C
0
⇓
0
=
det
⁡
α
β
+
i
1
+
α
-
λ
β
B
-
i
β
+
i
1
+
α
+
λ
C
0
⇓
0
=
det
⁡
α
i
-
λ
β
+
β
α
-
λ
λ
α
+
1
-
1
+
β
i
B
-
C
.Thus, μ defined as(10)
μ
=
α
i
-
λ
β
+
β
α
-
λ
λ
α
+
1
-
1
+
β
i
,is a generalized eigenvalue; that is, it satisfies det
⁡
(
μ
B
-
C
)
=
0 and it is well known that μ must be nonnegative. Computing λ from (10), the function λ
(
μ
) is found to be(11)
λ
μ
∶
=
α
+
i
μ
β
+
i
1
+
α
β
+
μ
,which allows evaluating an upper bound U
(
α
,
β
) of the spectral radius of P
-
1
Q:(12)
ϱ
(
P
-
1
Q
)
≤
sup
⁡
μ
≥
0
λ
μ
=
sup
⁡
μ
≥
0
α
2
+
μ
2
1
+
β
2
1
+
α
β
+
μ
≤
U
α
,
β
.There are optimal values of α and β that minimize U
(
α
,
β
) for α
≥
0 and β
≥
0. To find them, set α
=
l
sin
θ and β
=
l
cos
⁡
θ with θ
∈
[
0
,
π
/
2
]; then(13)
U
α
,
β
=
1
+
l
2
cos
⁡
θ
2
1
+
l
sin
θ
max
⁡
1
,
tan
θ
.If θ is fixed, the minimum of this last expression is for l
=
sin
θ
/
(
cos
⁡
θ
)
2 corresponding to(14)
U
α
,
β
=
cos
⁡
θ
max
⁡
1
,
tan
θ
=
max
⁡
cos
⁡
θ
,
sin
θ
for
θ
∈
0
,
π
2
.The minimum of U
(
α
,
β
) is attained for θ
=
π
/
4, which corresponds to α
=
β
=
1. The computation of U
(
1,1
) is then straightforward.In this setting, sinceα
=
β
=
1 are optimal, from now on it is assumed α
=
β
=
1 and therefore V
=
W
=
B. Using these values, the two-step method (2) is recast as the one step method:(15)
x
(
k
+
1
)
=
P
-
1
Q
x
k
+
b
=
x
(
k
)
+
P
-
1
b
-
A
x
k
,where the simplified expressions for P and Q are(16)
P
=
1
+
i
B
+
C
,
Q
=
C
+
i
B
.Notice that P is well defined and nonsingular provided that B
+
C is not singular. Thus the assumptions of Lemma 1 are weakened when α
=
β
=
1, resulting in the next corollary.Corollary 2.
LetB and C be semi-SPD with B
+
C not singular and P and Q as defined in (16); then the spectral radius of P
-
1
Q satisfies the upper bound ϱ
(
P
-
1
Q
)
≤
2
/
2.Remark 3.
The spectral radius of the iteration matrix is bounded independently of its size; thus, once the tolerance is fixed, the maximum number of iterations is independent of the size of the problem.The iterative method (15) can be reorganized in Algorithm 1.Algorithm 1:MHSS iterative solver for A
x
=
(
B
+
i
C
)
x
=
b.
r
←
b;  x
←
0;
while  r
>
ε
b  do
Solve (B
+
C
)
h
=
(
1
-
i
/
2
)
r;
x
←
x
+
h;
r
←
b
-
A
x;
end  whileRemark 4.
The MHSS iterative solver of Algorithm1 needs at each iteration the resolution of two real linear systems, respectively, for the real and imaginary part, whose coefficient matrices are SPD, namely, B
+
C. For small matrices this can be efficiently done by using Cholesky decomposition. For large and sparse matrices a preconditioned conjugate gradient method is mandatory.Although Algorithm1 can be used to solve the linear system (1), better performance is obtained using one or more steps of Algorithm 1 not for solving the linear system (1) but as a preconditioner for a faster conjugate gradient like iterative solver such as COCG or COCR. The convergence rate estimation for these iterative schemes for a linear system M
x
=
(
B
+
C
)
x
=
b depends on the condition number κ
2
=
M
2
M
-
1
2, where M
2
=
ϱ
(
M
T
M
) is the classic spectral norm of a matrix. The energy norm ·
M induced by the (real) SPD matrix M is used to obtain the well known estimate(17)
x
(
k
)
-
x
⋆
M
x
(
0
)
-
x
⋆
M
≤
2
κ
2
-
1
κ
2
+
1
k
,
v
M
=
v
T
M
v
,where x
⋆ is the solution of the linear system. In general, conjugate gradient-like iterative schemes perform efficiently if a good preconditioner makes the system well conditioned.


## 1.1. The MHSS Iterative Solver

The complexN
×
N linear system A
x
=
b, where A is complex symmetric, is solved via an iterative method based on a splitting algorithm (HSS). The preconditioner requires to solve (each time it is applied) two real symmetric and positive definite (SPD) linear systems.The HSS scheme can be summarized in this manner: rewrite the vectorx of the unknowns as the sum of real and imaginary part, x
=
y
+
i
z, and accordingly, the right hand side b
=
c
+
i
d, then the system A
x
=
b is rewritten as(1)
B
+
i
C
y
+
i
z
=
c
+
i
d
,so that the two steps of the modified HSS method proposed in [15] result in(2)
V
+
B
x
k
+
1
/
2
=
V
-
i
C
x
k
+
b
,
W
+
C
x
(
k
+
1
)
=
W
+
i
B
x
(
k
+
1
/
2
)
-
i
b
,for suitable matrices V and W. The previous procedure can be rewritten as a single step of a splitting based scheme P
x
k
+
1
=
Q
x
k
+
b by posing(3)
P
=
V
+
B
W
-
i
V
-
1
W
+
C
,
Q
=
V
+
B
W
-
i
V
-
1
W
+
i
B
V
+
B
-
1
V
-
i
C
.This iterative method converges if the iteration matrix P
-
1
Q, for example,(4)
P
-
1
Q
=
W
+
C
-
1
W
+
i
B
V
+
B
-
1
V
-
i
C
,has spectral radius strictly less than one. It is well known that the choice V
=
W
=
α
I, for a given positive constant α, yields the standard HSS method [15–17] for which an estimate of the spectral radius is given by(5)
ϱ
P
-
1
Q
≤
max
⁡
j
=
1,2
,
…
,
N
α
2
+
λ
j
B
2
α
+
λ
j
B
,where λ
j
(
B
) are the eigenvalues of the SPD matrix B. The optimal value for α can also be computed [15–17] and is(6)
α
opt
=
arg
min
⁡
α
max
⁡
j
=
1,2
,
…
,
n
α
2
+
λ
j
B
2
α
+
λ
j
B
α
opt
=
λ
min
⁡
B
λ
max
⁡
B
.From the previous formulas, it is clear that those computations rely on the knowledge (or estimate) of the minimum and maximum eigenvalue of the matrix B, which is, in general, a hard problem.Another possible choice forV and W is V
=
α
B and W
=
β
B which yields, when α
=
β, a variant of the MHSS method by [15–17]. In the next lemma, an upper bound of the spectral radius of the iteration matrix is given.Lemma 1.
LetV
=
α
B and W
=
β
B in (4) with B a SPD matrix and C a semi-SPD matrix, α
,
β
>
0; then the spectral radius of P
-
1
Q satisfies the upper bound(7)
ϱ
P
-
1
Q
≤
U
α
,
β
,
U
α
,
β
=
1
+
β
2
1
+
α
max
⁡
1
,
α
βand the minimum value of the upper bound U
(
α
,
β
) is attained when α
=
β
=
1 where U
(
1,1
)
=
2
/
2
≈
0.707.Proof.
SubstitutingV
=
α
B and W
=
β
B in (4) yields(8)
P
-
1
Q
=
β
B
+
C
-
1
β
B
+
i
B
α
B
+
B
-
1
α
B
-
i
C
P
-
1
Q
=
β
+
i
1
+
α
β
B
+
C
-
1
α
B
-
i
C
.If λ is an eigenvalue of matrix P
-
1
Q, it satisfies(9)
0
=
det
⁡
P
-
1
Q
-
λ
I
0
⇓
0
=
det
⁡
β
+
i
1
+
α
α
B
-
i
C
-
λ
β
B
+
C
0
⇓
0
=
det
⁡
α
β
+
i
1
+
α
-
λ
β
B
-
i
β
+
i
1
+
α
+
λ
C
0
⇓
0
=
det
⁡
α
i
-
λ
β
+
β
α
-
λ
λ
α
+
1
-
1
+
β
i
B
-
C
.Thus, μ defined as(10)
μ
=
α
i
-
λ
β
+
β
α
-
λ
λ
α
+
1
-
1
+
β
i
,is a generalized eigenvalue; that is, it satisfies det
⁡
(
μ
B
-
C
)
=
0 and it is well known that μ must be nonnegative. Computing λ from (10), the function λ
(
μ
) is found to be(11)
λ
μ
∶
=
α
+
i
μ
β
+
i
1
+
α
β
+
μ
,which allows evaluating an upper bound U
(
α
,
β
) of the spectral radius of P
-
1
Q:(12)
ϱ
(
P
-
1
Q
)
≤
sup
⁡
μ
≥
0
λ
μ
=
sup
⁡
μ
≥
0
α
2
+
μ
2
1
+
β
2
1
+
α
β
+
μ
≤
U
α
,
β
.There are optimal values of α and β that minimize U
(
α
,
β
) for α
≥
0 and β
≥
0. To find them, set α
=
l
sin
θ and β
=
l
cos
⁡
θ with θ
∈
[
0
,
π
/
2
]; then(13)
U
α
,
β
=
1
+
l
2
cos
⁡
θ
2
1
+
l
sin
θ
max
⁡
1
,
tan
θ
.If θ is fixed, the minimum of this last expression is for l
=
sin
θ
/
(
cos
⁡
θ
)
2 corresponding to(14)
U
α
,
β
=
cos
⁡
θ
max
⁡
1
,
tan
θ
=
max
⁡
cos
⁡
θ
,
sin
θ
for
θ
∈
0
,
π
2
.The minimum of U
(
α
,
β
) is attained for θ
=
π
/
4, which corresponds to α
=
β
=
1. The computation of U
(
1,1
) is then straightforward.In this setting, sinceα
=
β
=
1 are optimal, from now on it is assumed α
=
β
=
1 and therefore V
=
W
=
B. Using these values, the two-step method (2) is recast as the one step method:(15)
x
(
k
+
1
)
=
P
-
1
Q
x
k
+
b
=
x
(
k
)
+
P
-
1
b
-
A
x
k
,where the simplified expressions for P and Q are(16)
P
=
1
+
i
B
+
C
,
Q
=
C
+
i
B
.Notice that P is well defined and nonsingular provided that B
+
C is not singular. Thus the assumptions of Lemma 1 are weakened when α
=
β
=
1, resulting in the next corollary.Corollary 2.
LetB and C be semi-SPD with B
+
C not singular and P and Q as defined in (16); then the spectral radius of P
-
1
Q satisfies the upper bound ϱ
(
P
-
1
Q
)
≤
2
/
2.Remark 3.
The spectral radius of the iteration matrix is bounded independently of its size; thus, once the tolerance is fixed, the maximum number of iterations is independent of the size of the problem.The iterative method (15) can be reorganized in Algorithm 1.Algorithm 1:MHSS iterative solver for A
x
=
(
B
+
i
C
)
x
=
b.
r
←
b;  x
←
0;
while  r
>
ε
b  do
Solve (B
+
C
)
h
=
(
1
-
i
/
2
)
r;
x
←
x
+
h;
r
←
b
-
A
x;
end  whileRemark 4.
The MHSS iterative solver of Algorithm1 needs at each iteration the resolution of two real linear systems, respectively, for the real and imaginary part, whose coefficient matrices are SPD, namely, B
+
C. For small matrices this can be efficiently done by using Cholesky decomposition. For large and sparse matrices a preconditioned conjugate gradient method is mandatory.Although Algorithm1 can be used to solve the linear system (1), better performance is obtained using one or more steps of Algorithm 1 not for solving the linear system (1) but as a preconditioner for a faster conjugate gradient like iterative solver such as COCG or COCR. The convergence rate estimation for these iterative schemes for a linear system M
x
=
(
B
+
C
)
x
=
b depends on the condition number κ
2
=
M
2
M
-
1
2, where M
2
=
ϱ
(
M
T
M
) is the classic spectral norm of a matrix. The energy norm ·
M induced by the (real) SPD matrix M is used to obtain the well known estimate(17)
x
(
k
)
-
x
⋆
M
x
(
0
)
-
x
⋆
M
≤
2
κ
2
-
1
κ
2
+
1
k
,
v
M
=
v
T
M
v
,where x
⋆ is the solution of the linear system. In general, conjugate gradient-like iterative schemes perform efficiently if a good preconditioner makes the system well conditioned.


## 2. Use of MHSS as Preconditioner

In this section the effect of a fixed number of steps of Algorithm1, used as a preconditioner for the linear system (1), is analyzed in terms of the reduction of the condition number. Performing n steps of Algorithm 1 with x
(
0
)
=
0 is equivalent to compute x
(
n
)
=
P
n
-
1
b, where(18)
P
n
-
1
=
I
+
P
-
1
Q
+
P
-
1
Q
2
+
⋯
+
P
-
1
Q
n
-
1
P
-
1
.Matrix P
n can be interpreted as an approximation of matrix A
=
B
+
i
C.Thus, it is interesting to obtain an estimate of the condition number of the preconditioned matrixP
n
-
1
A in order to check the effect of MHSS when used as preconditioner.For the estimation, we need to recall some classical results about spectral radii and norms. For any matrixM and any matrix norm, Gelfand’s Formula connects norm and spectral radius [27, 28]:(19)
ϱ
M
=
limsup
⁡
k
→
∞
M
k
1
/
k
.Notice that when ϱ
(
M
)
<
1, for k large enough, M
k
<
1.Lemma 5.
LetA
=
P
-
Q so that P
-
1
Q is such that ϱ
(
P
-
1
Q
)
<
1; then for any ɛ
>
0 satisfying ϱ
(
P
-
1
Q
)
+
ɛ
<
1 there is an integer n
ɛ
>
0 such that(20)
κ
P
n
-
1
A
≤
1
+
ϱ
P
-
1
Q
+
ɛ
n
1
-
ϱ
P
-
1
Q
+
ɛ
n
,
n
≥
n
ɛ
,where κ
(
M
)
=
M
M
-
1 is the condition number with respect to the norm · and P
n is defined by (18).Proof.
Observe thatP
-
1
A
=
P
-
1
(
P
-
Q
)
=
I
-
P
-
1
Q; hence using (18) the preconditioned matrix becomes P
n
-
1
A
=
I
-
(
P
-
1
Q
)
n. From Gelfand’s Formula (19) there exists n
ɛ such that(21)
P
-
1
Q
n
≤
ϱ
P
-
1
Q
+
ɛ
n
<
1
∀
n
≥
n
ɛand from (21), by setting M
=
(
P
-
1
Q
)
n, the convergent series (see [29, 30]) gives a bound for the norm of the inverse(22)
I
-
M
-
1
=
∑
k
=
0
∞
‍
M
k
,
I
-
M
-
1
≤
∑
k
=
0
∞
‍
M
k
=
1
1
-
M
,
κ
P
n
-
1
A
=
κ
I
-
M
=
I
-
M
I
-
M
-
1
≤
1
+
M
1
-
M
.The thesis follows trivially from (21).Corollary 6.
LetA
=
P
-
Q so that P
-
1
Q has the property that ϱ
(
P
-
1
Q
)
<
1; then there exists a matrix norm · such that the conditioning number of the matrix P
n
-
1
A with respect to this norm satisfies(23)
κ
P
n
-
1
A
≤
1
+
0
.
8
n
1
-
0
.
8
n
≤
9
,where κ
M
=
M
M
-
1.Proof.
Recall that for any matrixM and ϵ
>
0 there exists a matrix norm · such that M
≤
ϱ
(
M
)
+
ϵ. This is a classical result of linear algebra; see, for example, Section 2.3 of [29] or Section 6.9 of [30]. Thus, given P
-
1
Q and choosing 0
<
ɛ
≤
0.8
-
2
/
2 there exists a matrix norm · such that P
-
1
Q
≤
0.8. The proof follows from Lemma 5.From Corollary2 using the Euclidean norm and choosing ɛ such that ϱ
(
P
-
1
Q
)
+
ɛ
=
2
/
2
+
ɛ
=
0.8 for n
≥
n
ɛ, the condition number of the preconditioned matrix satisfies(24)
κ
2
P
n
-
1
A
≤
1
+
0
.
8
n
1
-
0
.
8
n
,
κ
2
M
=
M
2
M
-
1
2
.This estimate shows that usingn steps of MHSS, with n large enough, the condition number of the preconditioned system can be bounded independently of the size of the linear system. In practice, when n
=
1 the reduction of the condition number is enough; in fact, Corollary 6 shows that, using the appropriate norm, the condition number of the preconditioned linear system is less than 9, independently of its size.Remark 7.
From [31], when the condition number κ is large, the estimate of m conjugate gradient (CG) iterations satisfies m
∝
κ. The cost of computation is proportional to the number of iterations, whereas the cost of each iteration is proportional to 1
+
C
n, where C is the cost of an iteration of MHSS used as preconditioner, relative to the cost of a CG step. Thus, using κ
2 from (24), when n is large enough, a rough estimate of the computational cost is(25)
cost
∝
κ
2
P
n
-
1
A
1
+
C
n
∝
1
+
0
.
8
n
1
-
0
.
8
n
1
+
C
n
.Fixing C, the cost (25), as a function of n alone, is convex and has a minimum for n
=
n
min
⁡
(
C
) which satisfies(26)
C
=
-
a
n
ln
⁡
a
a
n
ln
⁡
a
-
a
2
n
+
1
,
a
=
0.8
.The inverse function C
(
n
) which satisfies n
min
⁡
(
C
(
n
)
)
=
n is the right hand side of (26); moreover, C
(
n
) is a monotone decreasing function so that also n
min
⁡
(
C
) is monotone decreasing.
Table in equation (27) shows the constant C as a function of n giving the critical values of C such that n steps of MHSS are better than n
-
1 steps,(27)
n
min
⁡
1
2
3
4
5
6
7
C
(
n
min
⁡
)
0.98
0.47
0.29
0.2
0.14
0.1
0.07
.This means that it is convenient to use n
>
1 steps in the preconditioner MHSS only if the cost of one step of MHSS is less than 0.47, that is, half of one step of the conjugate gradient method, a situation that never happens in practice.According to Remark7, only P
n with n
=
1 is considered; that is, P
1
=
P as preconditioner for the linear system (1).


## 3. Iterative Solution of Complex Linear System

From the previous section, it is clear that the use of one iteration step of MHSS is a good choice that lowers the condition number of the original complex linear system (1). The resulting preconditioner matrix is P
=
(
1
+
i
)
(
B
+
C
) as defined in (16). Preconditioner P will be used together with semi-iterative methods specialized in the solution of complex problems. Examples of those methods include COCG [4, 32] or COCR [5, 32–34]. They are briefly exposed in Algorithms 2 and 3.Algorithm 2:COCG.
r
←
b
-
A
x;
r
~
←
P
-
1
r;
p
←
r
~;
ρ
←
[
r
~
,
r
];
while
r
>
ε
b
do
q
←
A
p;
μ
←
[
q
,
p
];
α
←
ρ
/
μ;
x
←
x
+
α
p;
r
←
r
-
α
q;
r
~
←
P
-
1
r;
β
←
ρ;
ρ
←
[
r
~
,
r
];
β
←
ρ
/
β;
p
←
r
~
+
β
p;
end whileAlgorithm 3:COCR (as in [32]).
r
←
b
-
A
x;
r
~
←
P
-
1
r;
p
←
r
~;
q
←
A
p;
ρ
←
[
r
~
,
q
];
while
r
>
ε
b
do
q
~
←
P
-
1
q;
α
←
ρ
/
q
~
,
q;
x
←
x
+
α
p;
r
←
r
-
α
q;
r
~
←
r
~
-
α
q
~;
t
←
A
r
~;
β
←
ρ;
ρ
←
r
~
,
t;
β
←
ρ
/
β;
p
←
r
~
+
β
p;
q
←
t
+
β
q
end whileThere are also other methods available for performing this task, for example, CSYM [6] or QMR [35].The application of the preconditionerP for those methods is equivalent to the solution of two real SPD systems depending on B
+
C, as in Remark 4. Of course, one can use a direct method [20, 36–38], or a conjugate gradient method [18] with incomplete Cholesky factorizations, or approximate inverse as a preconditioner [20, 39–43].For very large linear systemsB
+
C may be expensive or impractical to be formed; moreover, well suited variants of the incomplete Cholesky factorization or of the sparse inverse must be considered [44–46]. In the next section, a polynomial preconditioner based onorthogonal polynomials is presented; it will allow solving very large complex linear systems with a simple algorithm that can be implemented with few lines of code.The suggested strategy to get the solution of the complex linear system of the formA
x
=
(
B
+
i
C
)
x
=
b is to use an iterative method like COCG (Algorithm 2) or COCR (Algorithm 3) preconditioned with P
=
(
1
+
i
)
(
B
+
C
). The fully populated factorization of P
-
1 for large or huge linear systems is not practically computable so that some approximation must be used. The strategy can be split in two.(1)
Compute(
B
+
C
)
z
=
v
/
(
1
+
i
) by using iterative methods. This can be organized as the solution of two SPD real system, one for the real and one for the imaginary part. This systems can be efficiently solved by using preconditioned conjugate gradient.
(2)
Approximate(
B
+
C
)
-
1
=
Q
+
E with E “small” and use 1
+
i
-
1
Q
-
1 asapproximate preconditioner.Approximations of strategy 2 can be used as preconditioners for strategy 1. In strategy 2 the approximation must not be too far from the inverse of (
1
+
i
)
(
B
+
C
); otherwise the convergence will fail. In general for strategy 2 the choice of the preconditioner can be done as follows.(i)
If the complete Cholesky ofB
+
C
=
L
L
T is available use P
=
(
1
+
i
)
L
L
T as preconditioner.
(ii)
If the incomplete Cholesky ofB
+
C
=
L
~
L
~
T
+
E is computationally not too expensive and E is “small” then use P
=
(
1
+
i
)
L
~
L
~
T as preconditioner.
(iii)
IfE is “too large” or the solution of (
B
+
C
)
z
=
v
/
(
1
+
i
) by means of PCG with Q
=
L
~
L
~
T as preconditioner is too costly, try some alternative. In particular try:
(a)
Approximate(
B
+
C
)
-
1 using sparse inverse [39, 47–50] or algebraic multigrid [51–53] or any other good approximation.
(b)
Use the polynomial preconditioner proposed in Section5.Incomplete LU decomposition is easy to set up and exhibits good performances, while better performance may be obtained by using algebraic multigrid method (AMG); however, comparisons with AMG are difficult because of the large number of parameters to set up, like the depth and the shape of the W-cycle and the smoother (or relaxation method) and the coarse grid correction step.If incomplete Cholesky fails or it gives poor preconditioning, we propose a simple but effective polynomial preconditioner for which we give details such as stability and convergence properties. Notice that the optimal preconditioner is problem dependent and a good preconditioner for a problem may be a bad preconditioner for another one. The polynomial preconditioner proposed in Section5, without being optimal, is extremely simple to implement and easily parallelizable and ensures convergence also for very large linear systems. Finally, polynomial preconditioners have the following benefits with respect to incomplete Cholesky and algebraic multigrid.(i)
If matrix-vector multiplication(
B
+
i
C
)
z can be computed but the matrices B and C cannot be explicitly formed, the polynomial preconditioner is computable while incomplete LU and algebraic multigrid cannot be built.
(ii)
Polynomial preconditioners are easily parallelizable.
(iii)
The polynomial preconditioner never breaks down as it can happen to incomplete Cholesky or algebraic multigrid.
(iv)
The polynomial preconditioner can be extended also for problems with a singular matrix.Strategy number 1, that is, computing(
B
+
C
)
z
=
v
/
(
1
+
i
) by using iterative methods (e.g., PCG), is in general not competitive. In fact, to be competitive, the system should be solved with very few iterations, but this means that the preconditioner for B
+
C is very good; therefore it can be used directly and successfully in strategy number 2.


## 4. Scaling the Complex Linear System

The polynomial preconditioner presented in the next section depends on the knowledge of an interval containing eigenvalues. Scaling is a cheap procedure to recast the problem into one with eigenvalues in the interval(
0,1
]. Using the diagonal matrix S, the linear system (1) becomes(28)
S
A
S
w
=
S
B
S
+
i
S
C
S
w
=
S
b
,
with
x
=
S
w
,where S is a real diagonal matrix with positive entries on the diagonal. The scaled system inherits the properties of the original and still has the matrices S
B
S and S
C
S semi-SPD with S
B
S
+
S
C
S SPD. The next lemma shows how to choose a good scaling factor S used forward.Lemma 8.
LetM be a SPD matrix and S a diagonal matrix with S
i
i
=
∑
j
=
1
n
M
i
j
-
1
/
2; then the scaled matrix S
M
S has the eigenvalues in the range (
0,1
].Proof.
Notice thatS
M
S is symmetric and positively defined and is similar to S
2
M. Moreover, the estimate λ
max
⁡
(
S
2
M
)
≤
S
2
M
∞
=
1 follows trivially.Assumption 9.
From Lemma8, the linear system (1) is scaled to satisfy the following:(i)
matricesB and C are semi-SPD;
(ii)
matrixB
+
C is SPD with eigenvalues in (
0,1
].


## 5. Preconditioning with Polynomials

On the basis of the results of the previous sections with Assumption9, the linear system to be preconditioned has the form A
x
=
b with A
=
B
+
i
C with B and C semi-SPD and M
=
B
+
C SPD with eigenvalues distributed in the interval (
0,1
].A good preconditioner for this linear system is one step of MHSS in Algorithm1, which results in a multiplication by P
-
1 where P
=
(
1
+
i
)
(
B
+
C
)
=
(
1
+
i
)
M. Here the following polynomial approximation of P
-
1 is proposed:(29)
P
-
1
=
1
-
i
2
M
-
1
≈
1
-
i
2
s
m
(
M
)
.The matrix polynomial s
m
(
M
) must be an approximation of the inverse of M, that is, s
m
(
M
)
M
≈
I, where s
m
(
x
) is a polynomial with degree m. A measure of the quality of the preconditioned matrix for a generic polynomial s
(
x
) is the distance from the identity matrix:(30)
Q
σ
s
=
s
M
M
-
I
2
=
max
⁡
λ
∈
σ
(
M
)
1
-
λ
s
λ
,where σ
(
M
)
=
{
λ
1
,
…
,
λ
n
} is the spectrum of M. If, in particular, the preconditioned matrix s
m
(
M
)
M is the identity matrix then Q
σ
(
s
m
)
=
0. Thus, the polynomial preconditioner s
m should concentrate the eigenvalues of s
m
(
M
)
M around 1 in order to be effective.A preconditioner polynomial can be constructed by minimizingQ
σ
(
s
) of (30) within the space Π
m of polynomials of degree at most m. This implies the knowledge of the spectrum of the matrix M which is in general not available making problem (30) unfeasible. The following approximation of quality measure (30) is feasible:(31)
Q
ϵ
,
1
s
=
max
⁡
λ
∈
[
ϵ
,
1
]
1
-
λ
s
λ
,
σ
(
M
)
⊂
[
ϵ
,
1
]and it needs the knowledge of [
ϵ
,
1
] and an interval for ϵ
>
0, containing the spectrum of M. The polynomial which minimizes Q
[
ϵ
,
1
]
(
s
) for s
∈
Π
m is well known and is connected to an appropriately scaled and shifted Chebyshev polynomial. The construction of such solution is described in Section 5.1 and was previously considered by Ashby et al. [21, 22], Johnson et al. [23], Freund [54], Saad [24], and Axelsson [19]. The computation of Q
[
ϵ
,
1
]
(
s
) needs the estimation of a positive lower bound of the minimum eigenvalue of M, which is, in general, expensive or infeasible. The estimate ϵ
=
0 cannot be used because Q
[
0,1
]
(
s
)
≥
1 for any polynomial s. A different way to choose ϵ is analysed later in this section. Saad observed that the use of Chebyshev polynomials with the conjugate gradient method, that is, the polynomial which minimizes the condition number of the preconditioned system, is in general far from being the best polynomial preconditioner, that is, the one that minimizes the CG iterations [24]. Practice shows that although nonoptimal, Chebyshev preconditioners perform well in many situations. The following integral average quality measure proposed in [22–24, 55] is a feasible alternative to (31):(32)
Q
s
=
∫
0
1
‍
1
-
λ
s
λ
2
d
λ
,
σ
M
⊂
0,1
.The preconditioner polynomial s
m proposed here is the solution of minimization of quality measure (31) or (32):(33)
s
m
=
argmin
s
∈
Π
m
Q
ϵ
,
1
s
,
or
s
m
=
argmin
s
∈
Π
m
Q
s
.Solution of problem (33) is detailed in the next sections. The proposed solution to the first problem is by means of the Chebyshev polynomials, while the solution of the second problem is done with the Jacobi weight.


### 5.1. Chebyshev Polynomial Preconditioner

The solution of minimization problem (33) with quality measure Q
[
ϵ
,
1
]
(
s
) is well known and can be written in terms of Chebyshev polynomials [21, 22]:(34)
1
-
λ
s
m
λ
=
T
m
+
1
ϵ
λ
T
m
+
1
ϵ
0
,where T
m
+
1
ϵ
(
λ
) is the (
m
+
1
)th Chebyshev polynomial scaled in the interval [
ϵ
,
1
]. Polynomials T
k
ϵ
(
λ
) satisfy the recurrence(35)
T
0
ϵ
x
=
1
,
T
1
ϵ
(
x
)
=
a
x
+
b
,
T
n
+
1
ϵ
x
=
2
a
x
+
b
T
n
ϵ
x
-
T
n
-
1
ϵ
x
,where(36)
a
=
2
1
-
ϵ
,
b
=
-
1
+
ϵ
1
-
ϵ
.From (34), the preconditioner polynomial s
m
(
λ
) becomes(37)
s
m
λ
=
1
λ
1
-
T
m
+
1
ϵ
λ
T
m
+
1
ϵ
0and from (35) it is possible to give a recursive definition for s
m
(
λ
) too.Lemma 10 (recurrence formula for preconditioner).
Given the polynomialsq
n defined by the recurrence(38)
q
0
x
=
1
,
q
1
x
=
a
0
x
+
b
0
,
q
n
+
1
x
=
a
n
x
+
b
n
q
n
x
+
c
n
q
n
-
1
x
,
n
=
1,2
,
3
,
…then the polynomials r
n
(
x
)
=
q
n
(
x
)
/
q
n
(
0
) and s
n
(
x
)
=
(
1
-
r
n
+
1
(
x
)
)
/
x satisfy the recurrences(39)
r
0
x
=
1
,
r
1
x
=
1
-
a
0
′
x
,
r
n
+
1
x
=
a
n
′
x
+
b
n
′
r
n
x
+
c
n
′
r
n
-
1
x
,
s
0
x
=
a
0
′
,
s
1
x
=
a
1
′
x
+
b
1
′
,
s
n
x
=
a
n
′
x
+
b
n
′
s
n
-
1
x
+
c
n
′
s
n
-
2
x
-
a
n
′
,where(40)
a
0
′
=
-
a
0
b
0
,
a
1
′
=
-
a
0
a
1
b
0
b
1
+
c
1
,
b
1
′
=
-
a
0
b
1
+
a
1
b
0
b
0
b
1
+
c
1
,
a
n
′
=
a
n
γ
n
,
b
n
′
=
b
n
γ
n
,
c
n
′
=
c
n
γ
n
-
1
γ
n
andγ
n
=
q
n
(
0
)
/
q
n
+
1
(
0
) satisfies the recurrence(41)
γ
1
=
b
0
b
0
b
1
+
c
1
,
γ
n
=
1
b
n
+
c
n
γ
n
-
1
.Proof.
Take the ratio(42)
q
n
+
1
x
q
n
+
1
0
=
a
n
x
+
b
n
q
n
x
q
n
+
1
0
+
c
n
q
n
-
1
x
q
n
+
1
0
,
r
n
+
1
x
=
a
n
x
+
b
n
r
n
x
q
n
0
q
n
+
1
0
r
n
+
1
x
=
+
c
n
r
n
-
1
x
q
n
-
1
0
q
n
0
q
n
0
q
n
+
1
0
,
r
n
+
1
x
=
a
n
x
+
b
n
r
n
x
γ
n
+
c
n
r
n
-
1
x
γ
n
γ
n
-
1
,and notice that r
n
(
0
)
=
1 for all n. Recurrence for γ
n is trivially deduced. From r
n
+
1
(
x
)
=
1
-
x
s
n
(
x
) and by using (41),(43)
r
n
+
1
x
=
a
n
x
+
b
n
r
n
x
γ
n
+
c
n
r
n
-
1
x
γ
n
γ
n
-
1
,
1
-
x
s
n
x
=
a
n
x
+
b
n
1
-
x
s
n
-
1
x
γ
n
1
-
x
s
n
x
=
+
c
n
1
-
x
s
n
-
2
x
γ
n
γ
n
-
1
,
1
-
x
s
n
x
=
a
n
x
1
-
x
s
n
-
1
x
γ
n
+
b
n
γ
n
1
-
x
s
n
x
=
-
b
n
x
s
n
-
1
x
γ
n
+
c
n
γ
n
γ
n
-
1
1
-
x
s
n
x
=
-
c
n
x
s
n
-
2
x
γ
n
γ
n
-
1
,
-
x
s
n
x
=
a
n
x
γ
n
-
a
n
x
2
s
n
-
1
x
γ
n
-
b
n
x
s
n
-
1
x
γ
n
-
x
s
n
x
=
-
c
n
x
s
n
-
2
x
γ
n
γ
n
-
1
,
-
x
s
n
x
=
-
x
γ
n
a
n
x
+
b
n
s
n
-
1
x
+
c
n
s
n
-
2
x
γ
n
-
1
-
a
n
,dividing by -
x recurrence (39) is retrieved. Polynomials s
0 and s
1 are trivially computed.Using Lemma10 the polynomial preconditioner (34) satisfies the recurrence (39) with(44)
a
0
′
=
2
1
+
ϵ
,
a
1
′
=
-
8
ϵ
2
+
6
ϵ
+
1
,
b
1
′
=
8
1
+
ϵ
ϵ
2
+
6
ϵ
+
1
,
a
n
′
=
4
γ
n
1
-
ϵ
,
b
n
′
=
-
2
γ
n
1
+
ϵ
1
-
ϵ
,
c
n
′
=
-
γ
n
γ
n
-
1
,where c
n
=
-
1 is used and γ
n
=
T
n
ϵ
(
0
)
/
T
n
+
1
ϵ
(
0
) is computed by solving recurrence (35) for x
=
0; that is,(45)
T
n
ϵ
0
=
1
2
c
n
+
c
-
n
,
c
=
ϵ
-
1
ϵ
+
1
,
⟹
γ
n
=
T
n
ϵ
0
T
n
+
1
ϵ
0
=
c
n
+
c
-
n
c
n
+
1
+
c
-
n
+
1
.Numerical stability of recurrence (44) is discussed in Section 6. The estimation of ϵ is the complex task and some authors perform it dynamically. As an alternative, the present approach is to move the eigenvalues of the coefficient matrix from the interval [
ϵ
,
1
] to a stripe [
1
-
δ
,
1
+
δ
], so that the condition number remains bounded. The value of ϵ is not determined from the estimate of the eigenvalues but from the degree of the preconditioner polynomial and from the amplitude of the stripe δ. Once δ is fixed, the higher the degree of the preconditioner, the lower the value of ϵ, which decreases to zero. Thus, if the degree of the preconditioner is high enough, the eigenvalues are moved in the interval [
1
-
δ
,
1
+
δ
]. The important fact is that even if the degree is not high enough to move the complete spectrum, the majority of the eigenvalues are moved in the desired stripe, improving the performance of the conjugate gradient method. An idea of this behaviour is showed in Figure 1 on the right. The end of this section is devoted to the explicit expression of the value of ϵ computed backwards from the value of δ
n: once the maximum condition number is fixed, it is possible to increase the degree of the polynomial preconditioner so that ϵ decreases until the whole (or at least the most) spectrum of the matrix is contained in the specified range.Figure 1
The polynomials described in Lemma10: in the first row the product λ
s
m
(
λ
) showing the approximation of the identity; in the second row the explicit graph of the polynomial preconditioner compared with the function 1
/
λ; in the third row the performance of the preconditioner in terms of the degree, condition number, and concentration of the eigenvalues of the coefficient matrix around 1. The left column represents the Jacobi weight and the right column the Chebyshev polynomials.In the interval[
ϵ
,
1
], Chebyshev polynomial T
n
ϵ
(
x
)
/
T
n
ϵ
(
0
) is bounded in the range [
-
δ
,
δ
] where δ
=
T
n
ϵ
(
0
)
-
1
=
2
/
(
c
n
+
1
+
c
-
(
n
+
1
)
) and solving for ϵ gives(46)
ϵ
=
c
-
1
c
+
1
2
,
c
=
1
+
1
-
δ
2
δ
1
/
n
+
1
.


### 5.2. Jacobi Polynomial Preconditioner

The solution of minimization problem (33) with quality measure Q
(
s
) is well known and can be written in terms of Jacobi orthogonal polynomials [24].Definition 11.
Given a nonnegative weight functionw
(
λ
)
:
[
0,1
]
↦
R
+ (see [23, 24, 56]) the scalar product ·
,
·
w and the relative induced norm ·
w are defined as(47)
p
,
q
w
=
∫
0
1
‍
p
λ
q
λ
w
λ
d
λ
,
p
w
=
p
,
p
w
=
∫
0
1
‍
p
λ
2
w
λ
d
λ
,where p and q are continuous functions. The orthogonal polynomials with respect to the scalar product ·
,
·
w are the polynomials p
k
(
λ
) which satisfy p
k
,
p
j
w
=
0 if k
≠
j.The orthogonal polynomials with respect to the weight(48)
w
α
,
β
λ
=
1
-
λ
α
λ
β
,
for
α
,
β
>
-
1
,defined in the interval [
0,1
], are the Jacobi polynomials and they satisfy the recurrence (see [57]):(49)
p
0
α
,
β
x
=
1
,
p
1
α
,
β
x
=
a
0
α
,
β
x
+
b
0
α
,
β
,
p
n
+
1
α
,
β
x
=
a
n
α
,
β
x
+
b
n
α
,
β
p
n
α
,
β
x
+
c
n
α
,
β
p
n
-
1
α
,
β
x
,for(50)
a
n
α
,
β
=
1
,
b
n
α
,
β
=
-
1
2
1
+
β
2
-
α
2
2
n
+
α
+
β
2
n
+
1
+
α
+
β
,
c
n
α
,
β
=
-
n
n
+
α
n
+
β
n
+
α
+
β
2
n
-
1
+
α
+
β
2
n
+
1
+
α
+
β
2
n
+
α
+
β
2
.The class of polynomials of the form 1
-
λ
s
(
λ
) with s
∈
Π
m
-
1 can be thought as polynomials r
∈
Π
m
+
1 with r
(
0
)
=
1; thus, the minimization problem (33) for Q
(
s
) can be recast to the following constrained minimization for w
(
λ
)
≡
1:(51)
r
m
+
1
=
argmin
r
∈
Π
m
+
1
,
r
(
0
)
=
1
r
w
.The preconditioner polynomial is s
m
(
λ
)
=
λ
-
1
(
1
-
r
m
+
1
(
λ
)
). Polynomial r
m
+
1
(
λ
) is expanded by means of Jacobi orthogonal polynomials with α
=
β
=
0:(52)
r
m
+
1
λ
=
∑
k
=
0
m
+
1
‍
α
k
p
k
0,0
λ
.Making use of the property of orthogonality with respect to the scalar product, one has(53)
r
m
+
1
w
2
=
∑
k
=
0
m
+
1
‍
α
k
2
p
k
0,0
w
2
,
1
=
r
m
+
1
(
0
)
=
∑
k
=
0
m
+
1
‍
α
k
p
k
0,0
(
0
)
,and thus the constrained minimum problem (51) is recast as(54)
minimize
∑
k
=
0
m
+
1
‍
α
k
2
p
k
0,0
w
2
,
subject
to
∑
k
=
0
m
+
1
‍
α
k
p
k
0,0
0
=
1
.Problem (54) is solved by using Lagrange multiplier with first order conditions resulting in(55)
α
k
=
p
k
0
/
p
k
0,0
w
2
∑
k
=
0
m
+
1
p
k
0,0
0
2
/
p
k
0,0
w
2
,
⟹
r
m
+
1
λ
=
∑
k
=
0
m
+
1
p
k
0,0
0
p
k
0,0
λ
/
p
k
0,0
w
2
∑
k
=
0
m
+
1
p
k
0,0
0
2
/
p
k
0,0
w
2
.Remark 12.
The solution (55) is only formal but barely useful from a computational point of view, because it requires computing explicitly the least squares polynomial. In fact, it is well known that the evaluation of a polynomial of high degree is a very unstable process. To make solution (55) practical, it is mandatory to obtain a stable recurrence formula that allows evaluating polynomials even of very high degree, for example, 1000 or more.To find a recurrence for (55), it must be rewritten as a ratio of orthogonal polynomials as for the Chebyshev preconditioner (34). To this scope, some classical theorems and definitions on orthogonal polynomials are here recalled for convenience. Christoffel-Darboux formulas and Kernel Polynomials, here recalled without proofs (see [55, 57]), are used to build the recurrence.Theorem 13 (Christoffel-Darboux formulas).
Orthogonal polynomials with respect to the scalar product·
,
·
w share the following identities:(56)
∑
j
=
0
k
p
j
x
p
j
y
p
j
w
2
=
1
p
k
w
2
p
k
+
1
(
x
)
p
k
(
y
)
-
p
k
+
1
y
p
k
x
x
-
y
.Theorem 14 (Kernel Polynomials).
Given orthogonal polynomialsp
k
(
x
) with respect to the scalar product ·
,
·
w, that is, with respect to weight function w
(
x
), then the polynomials(57)
q
k
x
=
p
k
+
1
x
p
k
0
-
p
k
+
1
0
p
k
x
xare orthogonal polynomials with respect to the scalar product ·
,
·
w defined as p
,
q
w
=
∫
0
1
‍
p
(
x
)
q
(
x
)
x
w
(
x
)
d
x, that is, with respect to the weight function x
w
(
x
). Moreover q
0
(
x
)
=
1.With the formulas of Christoffel-Darboux (56) and x
=
λ, y
=
0, it is possible to rewrite (55) as(58)
r
m
+
1
λ
=
1
C
p
m
+
2
0,0
λ
p
m
+
1
0,0
0
-
p
m
+
2
0,0
0
p
m
+
1
0,0
λ
λ
,
C
=
p
m
+
1
0,0
w
2
∑
k
=
0
m
+
1
‍
p
k
0,0
0
2
p
k
0,0
w
2
.Using the Kernel Polynomials of this last theorem, expression (58) becomes(59)
r
m
+
1
λ
=
p
m
+
1
0,1
λ
p
m
+
1
0,1
0
,where p
k
0,1
(
x
) are orthogonal polynomials with respect to the weight w
(
λ
)
=
λ. In fact, the Kernel Polynomials with respect to λ
w
α
,
β
(
λ
)
=
w
α
,
β
+
1
(
λ
) satisfy q
α
,
β
(
x
)
=
p
α
,
β
+
1
(
x
).Preconditioner polynomial can be computed recursively using Lemma10, where coefficients a
n
′, b
n
′, c
n
′, and γ
n are computed from a
n
0,1, b
n
0,1, and c
n
0,1. Given(60)
a
n
0,1
=
1
,
b
n
0,1
=
-
1
2
1
+
1
2
n
+
1
2
n
+
3
,
c
n
0,1
=
-
n
n
+
1
4
2
n
+
1
2
,the values of a
n
′, b
n
′, c
n
′, and γ
n from Lemma 10 become(61)
a
0
′
=
3
2
,
a
1
′
=
-
10
3
,
b
1
′
=
4
,
γ
n
=
a
n
′
=
-
4
+
2
3
n
+
5
n
+
2
2
,
b
n
′
=
2
-
Δ
,
c
n
′
=
-
1
+
Δ
,
Δ
=
2
3
n
2
+
6
n
+
2
2
n
+
1
n
+
2
2
.The only difficulty of the previous coefficients computation lies in the recursive solution of γ
n, which is here omitted for conciseness but that is a linear three-term recurrence with polynomial coefficients. Notice that Δ
→
0; thus the limit value of the above coefficients is evident.


### 5.3.  Recurrence Formula for the Preconditioner

Looking at Algorithms2 and 3, the polynomial preconditioner s
m
(
M
) is applied to a vector, that is, s
m
(
M
)
v. Thus, to avoid matrix-matrix multiplication, by defining v
(
k
)
=
s
k
(
M
)
v and using Lemma 10 with (61), the following recurrence is obtained for s
m
(
M
)
v
=
v
(
m
):(62)
v
(
0
)
=
a
0
′
v
,
v
(
1
)
=
a
1
′
M
v
+
b
1
′
v
,
v
(
n
)
=
a
n
′
M
v
n
-
1
-
v
+
b
n
′
v
(
n
-
1
)
+
c
n
′
v
(
n
-
2
)
,where n
=
2,3
,
…
,
m and recurrence (62) is the proposed preconditioner with coefficients given by (44) for Chebyshev and (61) for Jacobi the polynomials. Equation (62) yields Algorithm 4.Algorithm 4:Application of preconditioner s
m
(
M
) to a vector v.
t
1
←
a
0
′
v;
y
←
a
1
′
M
v
+
b
1
′
v;
for
n
=
2,3
,
…
,
m
do
t
0
←
t
1;
t
1
←
y;
y
←
a
n
′
M
t
1
-
v
+
b
n
′
t
1
+
c
n
′
t
0;
end for
return
y;


## 5.1. Chebyshev Polynomial Preconditioner

The solution of minimization problem (33) with quality measure Q
[
ϵ
,
1
]
(
s
) is well known and can be written in terms of Chebyshev polynomials [21, 22]:(34)
1
-
λ
s
m
λ
=
T
m
+
1
ϵ
λ
T
m
+
1
ϵ
0
,where T
m
+
1
ϵ
(
λ
) is the (
m
+
1
)th Chebyshev polynomial scaled in the interval [
ϵ
,
1
]. Polynomials T
k
ϵ
(
λ
) satisfy the recurrence(35)
T
0
ϵ
x
=
1
,
T
1
ϵ
(
x
)
=
a
x
+
b
,
T
n
+
1
ϵ
x
=
2
a
x
+
b
T
n
ϵ
x
-
T
n
-
1
ϵ
x
,where(36)
a
=
2
1
-
ϵ
,
b
=
-
1
+
ϵ
1
-
ϵ
.From (34), the preconditioner polynomial s
m
(
λ
) becomes(37)
s
m
λ
=
1
λ
1
-
T
m
+
1
ϵ
λ
T
m
+
1
ϵ
0and from (35) it is possible to give a recursive definition for s
m
(
λ
) too.Lemma 10 (recurrence formula for preconditioner).
Given the polynomialsq
n defined by the recurrence(38)
q
0
x
=
1
,
q
1
x
=
a
0
x
+
b
0
,
q
n
+
1
x
=
a
n
x
+
b
n
q
n
x
+
c
n
q
n
-
1
x
,
n
=
1,2
,
3
,
…then the polynomials r
n
(
x
)
=
q
n
(
x
)
/
q
n
(
0
) and s
n
(
x
)
=
(
1
-
r
n
+
1
(
x
)
)
/
x satisfy the recurrences(39)
r
0
x
=
1
,
r
1
x
=
1
-
a
0
′
x
,
r
n
+
1
x
=
a
n
′
x
+
b
n
′
r
n
x
+
c
n
′
r
n
-
1
x
,
s
0
x
=
a
0
′
,
s
1
x
=
a
1
′
x
+
b
1
′
,
s
n
x
=
a
n
′
x
+
b
n
′
s
n
-
1
x
+
c
n
′
s
n
-
2
x
-
a
n
′
,where(40)
a
0
′
=
-
a
0
b
0
,
a
1
′
=
-
a
0
a
1
b
0
b
1
+
c
1
,
b
1
′
=
-
a
0
b
1
+
a
1
b
0
b
0
b
1
+
c
1
,
a
n
′
=
a
n
γ
n
,
b
n
′
=
b
n
γ
n
,
c
n
′
=
c
n
γ
n
-
1
γ
n
andγ
n
=
q
n
(
0
)
/
q
n
+
1
(
0
) satisfies the recurrence(41)
γ
1
=
b
0
b
0
b
1
+
c
1
,
γ
n
=
1
b
n
+
c
n
γ
n
-
1
.Proof.
Take the ratio(42)
q
n
+
1
x
q
n
+
1
0
=
a
n
x
+
b
n
q
n
x
q
n
+
1
0
+
c
n
q
n
-
1
x
q
n
+
1
0
,
r
n
+
1
x
=
a
n
x
+
b
n
r
n
x
q
n
0
q
n
+
1
0
r
n
+
1
x
=
+
c
n
r
n
-
1
x
q
n
-
1
0
q
n
0
q
n
0
q
n
+
1
0
,
r
n
+
1
x
=
a
n
x
+
b
n
r
n
x
γ
n
+
c
n
r
n
-
1
x
γ
n
γ
n
-
1
,and notice that r
n
(
0
)
=
1 for all n. Recurrence for γ
n is trivially deduced. From r
n
+
1
(
x
)
=
1
-
x
s
n
(
x
) and by using (41),(43)
r
n
+
1
x
=
a
n
x
+
b
n
r
n
x
γ
n
+
c
n
r
n
-
1
x
γ
n
γ
n
-
1
,
1
-
x
s
n
x
=
a
n
x
+
b
n
1
-
x
s
n
-
1
x
γ
n
1
-
x
s
n
x
=
+
c
n
1
-
x
s
n
-
2
x
γ
n
γ
n
-
1
,
1
-
x
s
n
x
=
a
n
x
1
-
x
s
n
-
1
x
γ
n
+
b
n
γ
n
1
-
x
s
n
x
=
-
b
n
x
s
n
-
1
x
γ
n
+
c
n
γ
n
γ
n
-
1
1
-
x
s
n
x
=
-
c
n
x
s
n
-
2
x
γ
n
γ
n
-
1
,
-
x
s
n
x
=
a
n
x
γ
n
-
a
n
x
2
s
n
-
1
x
γ
n
-
b
n
x
s
n
-
1
x
γ
n
-
x
s
n
x
=
-
c
n
x
s
n
-
2
x
γ
n
γ
n
-
1
,
-
x
s
n
x
=
-
x
γ
n
a
n
x
+
b
n
s
n
-
1
x
+
c
n
s
n
-
2
x
γ
n
-
1
-
a
n
,dividing by -
x recurrence (39) is retrieved. Polynomials s
0 and s
1 are trivially computed.Using Lemma10 the polynomial preconditioner (34) satisfies the recurrence (39) with(44)
a
0
′
=
2
1
+
ϵ
,
a
1
′
=
-
8
ϵ
2
+
6
ϵ
+
1
,
b
1
′
=
8
1
+
ϵ
ϵ
2
+
6
ϵ
+
1
,
a
n
′
=
4
γ
n
1
-
ϵ
,
b
n
′
=
-
2
γ
n
1
+
ϵ
1
-
ϵ
,
c
n
′
=
-
γ
n
γ
n
-
1
,where c
n
=
-
1 is used and γ
n
=
T
n
ϵ
(
0
)
/
T
n
+
1
ϵ
(
0
) is computed by solving recurrence (35) for x
=
0; that is,(45)
T
n
ϵ
0
=
1
2
c
n
+
c
-
n
,
c
=
ϵ
-
1
ϵ
+
1
,
⟹
γ
n
=
T
n
ϵ
0
T
n
+
1
ϵ
0
=
c
n
+
c
-
n
c
n
+
1
+
c
-
n
+
1
.Numerical stability of recurrence (44) is discussed in Section 6. The estimation of ϵ is the complex task and some authors perform it dynamically. As an alternative, the present approach is to move the eigenvalues of the coefficient matrix from the interval [
ϵ
,
1
] to a stripe [
1
-
δ
,
1
+
δ
], so that the condition number remains bounded. The value of ϵ is not determined from the estimate of the eigenvalues but from the degree of the preconditioner polynomial and from the amplitude of the stripe δ. Once δ is fixed, the higher the degree of the preconditioner, the lower the value of ϵ, which decreases to zero. Thus, if the degree of the preconditioner is high enough, the eigenvalues are moved in the interval [
1
-
δ
,
1
+
δ
]. The important fact is that even if the degree is not high enough to move the complete spectrum, the majority of the eigenvalues are moved in the desired stripe, improving the performance of the conjugate gradient method. An idea of this behaviour is showed in Figure 1 on the right. The end of this section is devoted to the explicit expression of the value of ϵ computed backwards from the value of δ
n: once the maximum condition number is fixed, it is possible to increase the degree of the polynomial preconditioner so that ϵ decreases until the whole (or at least the most) spectrum of the matrix is contained in the specified range.Figure 1
The polynomials described in Lemma10: in the first row the product λ
s
m
(
λ
) showing the approximation of the identity; in the second row the explicit graph of the polynomial preconditioner compared with the function 1
/
λ; in the third row the performance of the preconditioner in terms of the degree, condition number, and concentration of the eigenvalues of the coefficient matrix around 1. The left column represents the Jacobi weight and the right column the Chebyshev polynomials.In the interval[
ϵ
,
1
], Chebyshev polynomial T
n
ϵ
(
x
)
/
T
n
ϵ
(
0
) is bounded in the range [
-
δ
,
δ
] where δ
=
T
n
ϵ
(
0
)
-
1
=
2
/
(
c
n
+
1
+
c
-
(
n
+
1
)
) and solving for ϵ gives(46)
ϵ
=
c
-
1
c
+
1
2
,
c
=
1
+
1
-
δ
2
δ
1
/
n
+
1
.


## 5.2. Jacobi Polynomial Preconditioner

The solution of minimization problem (33) with quality measure Q
(
s
) is well known and can be written in terms of Jacobi orthogonal polynomials [24].Definition 11.
Given a nonnegative weight functionw
(
λ
)
:
[
0,1
]
↦
R
+ (see [23, 24, 56]) the scalar product ·
,
·
w and the relative induced norm ·
w are defined as(47)
p
,
q
w
=
∫
0
1
‍
p
λ
q
λ
w
λ
d
λ
,
p
w
=
p
,
p
w
=
∫
0
1
‍
p
λ
2
w
λ
d
λ
,where p and q are continuous functions. The orthogonal polynomials with respect to the scalar product ·
,
·
w are the polynomials p
k
(
λ
) which satisfy p
k
,
p
j
w
=
0 if k
≠
j.The orthogonal polynomials with respect to the weight(48)
w
α
,
β
λ
=
1
-
λ
α
λ
β
,
for
α
,
β
>
-
1
,defined in the interval [
0,1
], are the Jacobi polynomials and they satisfy the recurrence (see [57]):(49)
p
0
α
,
β
x
=
1
,
p
1
α
,
β
x
=
a
0
α
,
β
x
+
b
0
α
,
β
,
p
n
+
1
α
,
β
x
=
a
n
α
,
β
x
+
b
n
α
,
β
p
n
α
,
β
x
+
c
n
α
,
β
p
n
-
1
α
,
β
x
,for(50)
a
n
α
,
β
=
1
,
b
n
α
,
β
=
-
1
2
1
+
β
2
-
α
2
2
n
+
α
+
β
2
n
+
1
+
α
+
β
,
c
n
α
,
β
=
-
n
n
+
α
n
+
β
n
+
α
+
β
2
n
-
1
+
α
+
β
2
n
+
1
+
α
+
β
2
n
+
α
+
β
2
.The class of polynomials of the form 1
-
λ
s
(
λ
) with s
∈
Π
m
-
1 can be thought as polynomials r
∈
Π
m
+
1 with r
(
0
)
=
1; thus, the minimization problem (33) for Q
(
s
) can be recast to the following constrained minimization for w
(
λ
)
≡
1:(51)
r
m
+
1
=
argmin
r
∈
Π
m
+
1
,
r
(
0
)
=
1
r
w
.The preconditioner polynomial is s
m
(
λ
)
=
λ
-
1
(
1
-
r
m
+
1
(
λ
)
). Polynomial r
m
+
1
(
λ
) is expanded by means of Jacobi orthogonal polynomials with α
=
β
=
0:(52)
r
m
+
1
λ
=
∑
k
=
0
m
+
1
‍
α
k
p
k
0,0
λ
.Making use of the property of orthogonality with respect to the scalar product, one has(53)
r
m
+
1
w
2
=
∑
k
=
0
m
+
1
‍
α
k
2
p
k
0,0
w
2
,
1
=
r
m
+
1
(
0
)
=
∑
k
=
0
m
+
1
‍
α
k
p
k
0,0
(
0
)
,and thus the constrained minimum problem (51) is recast as(54)
minimize
∑
k
=
0
m
+
1
‍
α
k
2
p
k
0,0
w
2
,
subject
to
∑
k
=
0
m
+
1
‍
α
k
p
k
0,0
0
=
1
.Problem (54) is solved by using Lagrange multiplier with first order conditions resulting in(55)
α
k
=
p
k
0
/
p
k
0,0
w
2
∑
k
=
0
m
+
1
p
k
0,0
0
2
/
p
k
0,0
w
2
,
⟹
r
m
+
1
λ
=
∑
k
=
0
m
+
1
p
k
0,0
0
p
k
0,0
λ
/
p
k
0,0
w
2
∑
k
=
0
m
+
1
p
k
0,0
0
2
/
p
k
0,0
w
2
.Remark 12.
The solution (55) is only formal but barely useful from a computational point of view, because it requires computing explicitly the least squares polynomial. In fact, it is well known that the evaluation of a polynomial of high degree is a very unstable process. To make solution (55) practical, it is mandatory to obtain a stable recurrence formula that allows evaluating polynomials even of very high degree, for example, 1000 or more.To find a recurrence for (55), it must be rewritten as a ratio of orthogonal polynomials as for the Chebyshev preconditioner (34). To this scope, some classical theorems and definitions on orthogonal polynomials are here recalled for convenience. Christoffel-Darboux formulas and Kernel Polynomials, here recalled without proofs (see [55, 57]), are used to build the recurrence.Theorem 13 (Christoffel-Darboux formulas).
Orthogonal polynomials with respect to the scalar product·
,
·
w share the following identities:(56)
∑
j
=
0
k
p
j
x
p
j
y
p
j
w
2
=
1
p
k
w
2
p
k
+
1
(
x
)
p
k
(
y
)
-
p
k
+
1
y
p
k
x
x
-
y
.Theorem 14 (Kernel Polynomials).
Given orthogonal polynomialsp
k
(
x
) with respect to the scalar product ·
,
·
w, that is, with respect to weight function w
(
x
), then the polynomials(57)
q
k
x
=
p
k
+
1
x
p
k
0
-
p
k
+
1
0
p
k
x
xare orthogonal polynomials with respect to the scalar product ·
,
·
w defined as p
,
q
w
=
∫
0
1
‍
p
(
x
)
q
(
x
)
x
w
(
x
)
d
x, that is, with respect to the weight function x
w
(
x
). Moreover q
0
(
x
)
=
1.With the formulas of Christoffel-Darboux (56) and x
=
λ, y
=
0, it is possible to rewrite (55) as(58)
r
m
+
1
λ
=
1
C
p
m
+
2
0,0
λ
p
m
+
1
0,0
0
-
p
m
+
2
0,0
0
p
m
+
1
0,0
λ
λ
,
C
=
p
m
+
1
0,0
w
2
∑
k
=
0
m
+
1
‍
p
k
0,0
0
2
p
k
0,0
w
2
.Using the Kernel Polynomials of this last theorem, expression (58) becomes(59)
r
m
+
1
λ
=
p
m
+
1
0,1
λ
p
m
+
1
0,1
0
,where p
k
0,1
(
x
) are orthogonal polynomials with respect to the weight w
(
λ
)
=
λ. In fact, the Kernel Polynomials with respect to λ
w
α
,
β
(
λ
)
=
w
α
,
β
+
1
(
λ
) satisfy q
α
,
β
(
x
)
=
p
α
,
β
+
1
(
x
).Preconditioner polynomial can be computed recursively using Lemma10, where coefficients a
n
′, b
n
′, c
n
′, and γ
n are computed from a
n
0,1, b
n
0,1, and c
n
0,1. Given(60)
a
n
0,1
=
1
,
b
n
0,1
=
-
1
2
1
+
1
2
n
+
1
2
n
+
3
,
c
n
0,1
=
-
n
n
+
1
4
2
n
+
1
2
,the values of a
n
′, b
n
′, c
n
′, and γ
n from Lemma 10 become(61)
a
0
′
=
3
2
,
a
1
′
=
-
10
3
,
b
1
′
=
4
,
γ
n
=
a
n
′
=
-
4
+
2
3
n
+
5
n
+
2
2
,
b
n
′
=
2
-
Δ
,
c
n
′
=
-
1
+
Δ
,
Δ
=
2
3
n
2
+
6
n
+
2
2
n
+
1
n
+
2
2
.The only difficulty of the previous coefficients computation lies in the recursive solution of γ
n, which is here omitted for conciseness but that is a linear three-term recurrence with polynomial coefficients. Notice that Δ
→
0; thus the limit value of the above coefficients is evident.


## 5.3.  Recurrence Formula for the Preconditioner

Looking at Algorithms2 and 3, the polynomial preconditioner s
m
(
M
) is applied to a vector, that is, s
m
(
M
)
v. Thus, to avoid matrix-matrix multiplication, by defining v
(
k
)
=
s
k
(
M
)
v and using Lemma 10 with (61), the following recurrence is obtained for s
m
(
M
)
v
=
v
(
m
):(62)
v
(
0
)
=
a
0
′
v
,
v
(
1
)
=
a
1
′
M
v
+
b
1
′
v
,
v
(
n
)
=
a
n
′
M
v
n
-
1
-
v
+
b
n
′
v
(
n
-
1
)
+
c
n
′
v
(
n
-
2
)
,where n
=
2,3
,
…
,
m and recurrence (62) is the proposed preconditioner with coefficients given by (44) for Chebyshev and (61) for Jacobi the polynomials. Equation (62) yields Algorithm 4.Algorithm 4:Application of preconditioner s
m
(
M
) to a vector v.
t
1
←
a
0
′
v;
y
←
a
1
′
M
v
+
b
1
′
v;
for
n
=
2,3
,
…
,
m
do
t
0
←
t
1;
t
1
←
y;
y
←
a
n
′
M
t
1
-
v
+
b
n
′
t
1
+
c
n
′
t
0;
end for
return
y;


## 6. Numerical Stability

Algorithm4, that is, the application of preconditioner s
m
(
M
) given by (62) to a vector v, also taking into account rounding errors, results in(63)
w
(
0
)
=
a
0
′
v
+
ϱ
(
0
)
,
w
(
1
)
=
a
1
′
M
v
+
b
1
′
v
+
ϱ
(
1
)
,
w
(
n
)
=
a
n
′
M
w
n
-
1
-
v
+
b
n
′
w
(
n
-
1
)
+
c
n
′
w
(
n
-
2
)
+
ϱ
(
n
)
,where ϱ
(
k
)
∞
≤
δ are the errors due to floating point operations with δ as an upper bound of such errors. The cumulative error e
(
k
)
=
w
(
k
)
-
v
(
k
) satisfies the linear recurrence(64)
e
(
0
)
=
ϱ
(
0
)
,
e
(
1
)
=
ϱ
(
1
)
,
e
(
n
)
=
a
n
′
M
e
(
n
-
1
)
+
b
n
′
e
(
n
-
1
)
+
c
n
′
e
(
n
-
2
)
+
ϱ
(
n
)
.The next definitions introduce the concept of generalized and joint spectral radius needed for the proof of the theorem of the matrix bound; they can be found in [58].Definition 15.
A matrix setΣ
=
{
A
k
∈
R
n
×
n
∣
k
∈
N
} is bounded if there is a constant C such that A
≤
C for all A
∈
Σ. An invariant subspace V for Σ is a vector space such that A
V
⊆
V for all A
∈
Σ. The set Σ is irreducible if the only invariant subspace is {
0
} or R
n.Definition 16.
The generalized spectral radiusϱ
(
Σ
) and the joint spectral radius ϱ
^
(
Σ
,
·
) of any set of matrices Σ are defined as(65)
ϱ
Σ
=
limsup
⁡
k
→
∞
ϱ
k
Σ
1
/
k
,
ϱ
k
Σ
=
sup
⁡
ϱ
∏
i
=
1
k
‍
A
i
∣
A
i
∈
Σ
,
ϱ
^
Σ
,
·
=
limsup
⁡
k
→
∞
ϱ
^
k
Σ
,
·
1
/
k
,
ϱ
^
k
Σ
,
·
=
sup
⁡
∏
i
=
1
k
‍
A
i
∣
A
i
∈
Σ
.Theorem 17.
LetΣ be a bounded and irreducible set of matrices with ϱ
(
Σ
)
>
0; then there is a constant C such that(66)
A
1
A
2
⋯
A
k
≤
C
ϱ
Σ
k
,
∀
A
j
∈
Σfor all k
>
0.Proof.
It is Theorem 2.1 by [58] with a slight modification to match the present case.Theorem 18.
Recurrence (64) satisfies(67)
e
(
n
)
∞
≤
C
δ
N
n
,where N is the dimension of the linear system, δ is the amplitude of the stripe for the eigenvalues, and C is an unknown constant coming from the norm inequalities which is found experimentally to be small.Proof.
The matrixM
=
B
+
C in (64), by Assumption 9, is SPD with eigenvalues in (
0,1
]. Thus M
=
T
T
Λ
T, with T orthogonal matrix, that is, T
T
T
=
I, and Λ diagonal. Multiplying on the left the recurrence (64) by T, the following error estimate is obtained:(68)
T
e
(
0
)
=
T
ϱ
(
0
)
,
T
e
(
1
)
=
T
ϱ
(
1
)
,
T
e
(
n
)
=
a
n
′
Λ
T
e
(
n
-
1
)
+
b
n
′
T
e
(
n
-
1
)
+
c
n
′
γ
n
-
1
T
e
(
n
-
2
)
+
T
ϱ
(
n
)
.Focusing on jth component of the transformed error, f
(
n
)
=
(
T
e
(
n
)
)
j and η
(
n
)
=
(
T
ϱ
(
n
)
)
j, a scalar recurrence is obtained:(69)
f
(
0
)
=
η
(
0
)
,
f
(
1
)
=
η
(
1
)
,
f
(
n
)
=
a
n
′
λ
j
+
b
n
′
f
(
n
-
1
)
+
c
n
′
f
(
n
-
2
)
+
η
(
n
)
.Recurrence (69) is restated in matrix form as(70)
f
n
=
A
n
f
n
-
1
+
b
n
,
A
n
=
a
n
′
λ
j
+
b
n
′
c
n
′
1
0
,
b
n
=
η
(
n
)
0
,
f
n
=
f
(
n
)
f
(
n
-
1
)
,with initial data f
1
T
=
(
η
(
1
)
,
η
(
0
)
). Notice that η
(
n
) is bounded by(71)
η
n
≤
T
ϱ
n
∞
≤
T
ϱ
n
2
η
n
=
ϱ
(
n
)
2
≤
N
ϱ
(
n
)
∞
≤
δ
Nand thus, f
1
∞
≤
δ
N and b
n
∞
≤
δ
N. From (70) it follows that(72)
f
n
=
A
n
A
n
-
1
⋯
A
2
f
1
+
∑
k
=
2
n
‍
A
n
A
n
-
1
⋯
A
n
-
k
+
1
b
k
.
The setΣ
=
{
A
i
∣
i
=
1
,
…
,
∞
} is bounded and irreducible; each matrix has spectral radius strictly less than 1 (see Lemma 19 for a proof); therefore the joint spectral radius is less than 1. From (72), with Theorem 17 using the infinity norm,(73)
f
n
≤
f
n
∞
≤
A
n
A
n
-
1
⋯
A
2
∞
f
1
∞
f
n
≤
f
n
∞
+
∑
k
=
2
n
‍
A
n
A
n
-
1
⋯
A
n
-
k
+
1
∞
b
k
∞
,
f
n
≤
f
n
∞
≤
C
ϱ
Σ
n
-
2
δ
N
+
C
∑
k
=
2
n
‍
ϱ
Σ
k
δ
N
f
n
≤
f
n
∞
≤
C
δ
N
n
,and, because of f
(
n
)
=
(
T
e
(
n
)
)
j, it follows that T
e
(
n
)
∞
≤
C
δ
N
n. A bound of the term e
(
n
) is done as(74)
e
(
n
)
∞
≤
e
(
n
)
2
=
T
e
(
n
)
2
≤
N
T
e
(
n
)
∞
≤
C
δ
N
n
.This shows that the error grows at most linearly.The above relation shows that the recurrence is at worst linearly unstable; that is, the error grows at most linearly. The existence is proved in the works of Rota and Strang [59] where the concept of joint spectral radius is introduced. The determination of C is not possible but practice reveals that it is small. In conclusion it is possible to employ even a very high degree polynomial preconditioner with a stable computation.Lemma 19.
Givena
n
′, b
n
′, c
n
′ from (44) or (61), respectively, for the Chebyshev and the Jacobi preconditioner, then the roots z
1 and z
2 of the characteristic polynomial of homogeneous recurrence (69), that is,(75)
z
2
-
a
n
′
λ
+
b
n
′
z
-
c
n
′
,satisfy z
1
<
1 and z
2
<
1 for all n
>
0 and 0
<
λ
≤
1.Proof.
Consider first the coefficients for the Jacobi polynomials defined in (61). If the roots are complex, then they must be conjugate; thus z
1
=
z and z
2
=
z
¯, because the coefficients of the polynomial are real. In that case, the constant term of the polynomial is equal to the square of the modulus of the roots, z
z
¯
=
|
z
|
2
=
-
c
n
′; thus it is easy to see that |
z
|
<
1 for all n
>
0. Suppose now that the two roots z
1 and z
2 are real; multiplying the characteristic polynomial by n
+
2
2
(
2
n
+
1
), yields, after some manipulation,(76)
z
1
=
A
-
B
λ
+
B
2
λ
2
-
2
A
B
λ
+
C
D
,
z
2
=
A
-
B
λ
-
B
2
λ
2
-
2
A
B
λ
+
C
Dfor(77)
A
=
2
n
3
+
6
n
2
+
6
n
+
2
,
B
=
4
n
3
+
12
n
2
+
11
n
+
3
,
C
=
3
n
2
+
6
n
+
2
2
,
D
=
2
n
3
+
9
n
2
+
12
n
+
4
,with A, B, C, and D strictly positive for all n
≥
0. The discriminant Δ
(
λ
) of the equation is Δ
(
λ
)
=
B
2
λ
2
-
2
A
B
λ
+
C and represents a convex parabola because B
2
>
0. Its minimum is obtained for λ
=
A
/
B
∈
(
0,1
), which gives Δ
(
A
/
B
)
<
0 and so complex roots, but this case was already considered. Hence we can set Δ
min
⁡
=
0. The maximum of Δ
(
λ
) is achieved at one of the extrema of the interval of definition of λ. A quick calculation shows that Δ
(
λ
) is maximum for λ
=
0, yielding a value of Δ
max
⁡
=
C. Using Δ
max
⁡ and Δ
min
⁡ it is possible to bound the roots z
1 and z
2:(78)
z
1
<
A
+
Δ
max
⁡
D
z
1
=
2
n
3
+
6
n
2
+
6
n
+
2
+
3
n
2
+
6
n
+
2
2
n
3
+
9
n
2
+
12
n
+
4
=
1
,
z
1
≥
A
-
B
+
Δ
min
⁡
D
=
-
2
n
3
+
6
n
2
+
5
n
+
1
2
n
3
+
9
n
2
+
12
n
+
4
>
-
1
,
z
2
<
A
-
Δ
min
⁡
D
=
2
n
3
+
6
n
2
+
6
n
+
2
2
n
3
+
9
n
2
+
12
n
+
4
<
1
,
z
2
≥
A
-
B
-
Δ
max
⁡
D
=
-
2
n
3
+
9
n
2
+
11
n
+
3
2
n
3
+
9
n
2
+
12
n
+
4
>
-
1
.The previous inequalities prove the lemma for the Jacobi preconditioner. Now consider the case of the coefficients of the Chebyshev polynomials defined in (44). Recall the expression for c
=
(
ϵ
-
1
)
/
(
ϵ
+
1
), and notice that, for ϵ
∈
(
0,1
), c is bounded in -
1
<
c
<
0, so that(79)
ω
n
∶
=
c
n
+
c
-
nis positive for even n and negative for odd n; moreover, ω
n
=
c
n
+
c
-
n
≥
2 is monotone increasing for n
=
1,2
,
3
,
…. In the case of complex roots, it was already shown that z
z
¯
=
|
z
|
2
=
-
c
n
′, with(80)
0
≤
-
c
n
′
=
ω
n
-
1
ω
n
+
1
<
1
,
ϵ
∈
0,1
.In the rest of the proof it is useful to consider also the ratio -
1
<
ω
n
/
ω
n
+
1
<
0, for c
∈
(
-
1,0
), that corresponds to ϵ
∈
(
0,1
). The coefficients of (44), observing that ϵ
=
(
ω
1
+
2
)
/
(
ω
1
-
2
), are simplified in(81)
a
n
′
=
4
ω
n
1
-
ϵ
ω
n
+
1
=
-
ω
n
ω
n
+
1
ω
1
-
2
,
b
n
′
=
-
2
1
+
ϵ
ω
n
1
-
ϵ
ω
n
+
1
=
ω
n
ω
1
ω
n
+
1
,
c
n
′
=
-
ω
n
-
1
ω
n
+
1
.Polynomial (75) is rewritten as(82)
z
2
+
ω
n
ω
n
+
1
λ
ω
1
-
2
-
ω
1
z
+
ω
n
-
1
ω
n
+
1and its roots are (using ω
1
ω
n
=
ω
n
-
1
+
ω
n
+
1)(83)
z
1,2
λ
=
ω
1
ω
n
2
ω
n
+
1
-
ω
n
ω
n
+
1
ω
1
-
2
2
λ
±
Δ
λ
,
Δ
λ
=
λ
2
-
2
ω
1
ω
1
-
2
λ
+
ω
n
-
1
-
ω
n
+
1
2
ω
n
2
ω
1
-
2
2
.Looking at the discriminant Δ
(
λ
), the minimum of the associated convex parabola is for λ
=
ω
1
/
ω
1
-
2
∈
(
1
/
2,1
). The corresponding value is(84)
Δ
ω
1
ω
1
-
2
=
-
4
ω
n
-
1
ω
n
+
1
ω
n
2
ω
1
-
2
2
<
0
.The value at the right extremum is also negative:(85)
Δ
1
=
4
ω
n
2
-
ω
n
-
1
ω
n
+
1
ω
n
2
ω
1
-
2
2
<
0
,and in fact(86)
ω
n
2
-
ω
n
-
1
ω
n
+
1
=
2
-
c
-
2
-
c
2
=
-
c
-
1
2
c
+
1
2
c
2
<
0
.For λ
=
0, with some manipulations, the roots of (82) are(87)
z
1
0
=
1
,
z
2
0
=
ω
n
-
1
ω
n
+
1
<
1
.Moreover, Δ
(
ϵ
)
=
Δ
(
1
); thus there exists 0
<
λ
⋆
<
ϵ such that Δ
(
λ
⋆
)
=
0. Thus for λ
∈
[
λ
⋆
,
1
] the roots are complex conjugate and with modulus less than 1. For λ
∈
[
0
,
λ
⋆
], where the roots z
1,2
(
λ
) are real, z
1
(
λ
) and its derivative satisfy(88)
z
1
λ
=
ω
1
ω
n
2
ω
n
+
1
+
ω
n
ω
n
+
1
ω
1
-
2
2
Δ
λ
-
λ
,
z
1
′
λ
=
ω
n
ω
n
+
1
ω
1
-
2
2
Δ
′
λ
Δ
λ
-
1
,
Δ
′
λ
=
2
λ
-
2
ω
1
ω
1
-
2
,and thus for λ
=
0 we have z
1
′
(
0
)
<
0. Hence in a neighbourhood of λ
=
0, -
1
<
z
1,2
(
λ
)
<
1, and for λ
∈
(
0
,
λ
⋆
] there are no roots equal to 1 or 0; thus the roots of (82) are bounded in (
0,1
). In fact, by contradiction, let z
=
1 be a root; then by (82)(89)
1
+
ω
n
ω
n
+
1
λ
ω
1
-
2
-
ω
1
+
ω
n
-
1
ω
n
+
1
=
0
,
⟹
λ
=
0
.Moreover z
=
0 is never a root of (82).
Thus the roots are bounded in the interval(
-
1,1
) for n
≥
0 and λ
∈
(
0,1
].


## 7. Numerical Tests

In this section a group of tests is proposed for the solution of a complex linear system of the form (1); that is, (
B
+
i
C
)
(
y
+
i
z
)
=
c
+
i
d, where B and C are semi-SPD with B
+
C SPD. The solvers used are COCG (Algorithm 2), COCR (Algorithm 3), and the standard Matlab’s QMR. The first two are tested together with different preconditioners such as ILU/ILU0 and with the present polynomial preconditioners. As a comparison QMR is also added, with and without preconditioning. When feasible, the direct LU factorization is used; then two variants are proposed, an incomplete LU and an ILU with threshold. More in detail, Table 1 presents the preconditioners:(i)
LU, the complete LU decomposition, for the matrixA
=
B
+
i
C,
(ii)
ILU0, the incomplete ILU(0), for the matrixA
=
B
+
i
C, the threshold used is 1
0
-
5,
(iii)
ILU, the incomplete ILU(0), for the matrixA
=
B
+
i
C,and the use of the previous preconditioners with COCG, COCR, and QMR:(i)
COCG-ILU0 and COCG-ILU, COCG solver preconditioned with ILU0 and ILU,
(ii)
COCR-ILU0 and COCR-ILU, COCR solver preconditioned with ILU0 and ILU,
(iii)
QMR, standard QMR implementation of Matlab, with no preconditioning,
(iv)
QMR-ILU0 and QMR-ILU, standard QMR implementation of Matlab, with ILU0 and ILU preconditioning.Table 2 presents the results of COCG and COCR iterative solver with the preconditioner MHSS approximated with(i)
COCG-MHSS-LL
T, COCG-MHSS-IC0, and COCG-MHSS-IC the complete and incomplete Cholesky decompositions for the preconditioner P defined in (16), threshold used for IC being 1
0
-
5;
(ii)
COCR-MHSS-LL
T, COCR-MHSS-IC0, and COCR-MHSS-IC the complete and incomplete Cholesky decompositions for the preconditioner P defined in (16), threshold used for IC being 1
0
-
5.Table 3 presents the results of COCG and COCR iterative solver with the preconditioner of Section 5:(i)
COCG-MHSS-J, COCG-MHSS-C the Jacobi and Chebyshev polynomial approximation of MHSS used with COCG iterative method;
(ii)
COCR-MHSS-J, COCR-MHSS-C the Jacobi and Chebyshev polynomial approximation of MHSS used with COCR iterative method;The Incomplete LU decomposition was computed with the standard Matlab commandilu with parameters ’type’ = ’crout’, ’droptol=1e-5’; the ILU0 was computedilu with parameters ’type’ = ’nofill’; for the Cholesky decomposition was usedichol with parameters ’type’ = ’ict’, ’droptol=1e-5’; for the IC0 the parameter used was ’type’=’nofill’. The degrees used for the polynomial preconditioner are 10, 50, 100, 500, and 1000. Due to the lack of complex symmetric matrices with SPD real and imaginary part, it was decided to combine two real SPD matrices of not too far dimension (eventually padding with zeros to match the size of the biggest one). The real SPD matrices used are summarized in Table 4 and can be found on the NIST “Matrix Market” Sparse Matrix Collection [60] or on University of Florida Sparse Matrix Collection [61]. As usual “NNZ” means number of nonzero elements and it is understood that the matrices are square; hence only the number of rows is reported.Table 1
Numerical results for COCG, COCR, and QMR preconditioned on the basis of the complex matrixA with LU decomposition, with incomplete LU (ILU) and with incomplete LU with threshold 10
-
5 (ILU0). The reported numbers represent the iterations of the corresponding solver with the specified preconditioner. The dash indicates that it was not feasible to compute a particular test. The column “preconditioner” shows the time required to assemble the preconditioner. The stopping tolerance was 10
-
8. The time elapsed in the computation is expressed in seconds.
Test case
Preconditioner
COCG
COCR
QMR
LU
ILU0
ILU
ILU0
ILU
ILU0
ILU
No Preco
ILU0
ILU
T5k
Iter.
86
9
86
9
85
8
Time
0.59
0.027
1.28
0.21
1.34
0.22
1.34
—
0.65
1.73
T16k
Iter.
65
6
65
6
399
64
5
Time
13.9
0.08
315
0.44
317
0.44
317
1.75
1.35
337
T80k
Iter.
204
205
200
Time
21.8
0.26
49.6
—
73
—
74
—
—
266
T150k
Iter.
661
24
634
24
642
23
Time
13.0
0.058
152
12.7
156
13.1
156
—
40.4
174
T500k
Iter.
Time
—
0.79
—
—
—
—
—
—
—
—
S13k
Iter.
1253
1258
2152
1307
Time
—
0.032
6.15
—
36.5
—
36.7
—
29.4
333
S32k
Iter.
Time
—
0.082
28
—
—
—
—
—
—
—
S500k
Iter.
Time
—
—
—
—
—
—
—
—
—
—Table 2
Numerical results for COCG and COCR preconditioned on the basis of the matrix (B + C) with complete Cholesky decomposition (L
L
T), with incomplete Cholesky (IC), and with incomplete Cholesky with threshold 10
-
5 (IC0). The reported numbers represent the iterations of the corresponding solver with the specified preconditioner. The dash indicates that it was not feasible to compute a particular test. The stopping tolerance was 10
-
8. The column “preconditioner” shows the time required to assemble the preconditioner. The time elapsed in the computation is expressed in seconds.
Test case
Preconditioner
COCG-MHSS
COCR-MHSS
L
L
T
IC0
IC
L
L
T
IC0
IC
L
L
T
IC0
IC
T5k
Iter.
25
148
25
25
148
25
Time
0.21
0.01
0.25
0.51
0.3
0.48
0.52
0.32
0.42
T16k
Iter.
24
78
24
24
75
24
Time
25
0.031
22
34.7
0.48
23.4
36
0.64
23.9
T80k
Iter.
22
22
Time
5.3
—
—
11
—
—
11.2
—
—
T150k
Iter.
638
32
634
32
Time
—
0.03
3.0
—
11
5.6
—
11.2
5.8
T500k
Iter.
Time
—
—
—
—
—
—
—
—
—
S13k
Iter.
5
5
Time
0.66
—
—
0.88
—
—
0.85
—
—
S32k
Iter.
Time
—
—
—
—
—
—
—
—
—
S500k
Iter.
Time
—
—
—
—
—
—
—
—
—Table 3
Numerical results for COCG and COCR, the reported numbers represents the iterations of the corresponding solver with the specified preconditioner. The dash indicates that it was not feasible to compute a particular test, while the letters “—” mean “not converged” after 10 000 iterations. The value ofδ used in the Chebyshev preconditioner was 0.2 while the stopping tolerance was 10
-
8. The time elapsed in the computation is expressed in seconds.
Test case
Degree of MHSS-J
Degree of MHSS-C
10
50
100
500
1000
10
50
100
500
1000
COCG
T5k
Iter.
127
32
26
25
25
111
32
32
31
31
Time
0.9
1.0
1.6
7.4
15.0
0.8
1.0
1.9
9.2
18.4
T16k
Iter.
31
24
24
23
23
31
29
28
29
30
Time
0.6
1.9
3.6
16.6
32.9
0.6
2.3
4.2
21.1
43.2
T80k
Iter.
4765
657
298
55
31
3904
596
251
52
33
Time
275.6
162.8
145.4
130.4
146.4
225.8
148.0
121.6
123.8
156.6
T150k
Iter.
323
81
42
19
19
277
65
36
23
23
Time
14.2
15.1
15.5
33.9
67.4
12.9
12.4
13.3
41.1
81.9
T500k
Iter.
4928
829
479
91
47
4805
929
461
87
46
Time
1665
1215
1374
1378
1393
1763
1333
1343
1246
1380
S13k
Iter.
4681
1123
591
80
38
4016
1199
613
93
57
Time
29.6
29.1
30.1
19.8
23.3
25.6
33.2
32.0
24.1
28.6
S32k
Iter.
4665
1102
686
106
49
4012
1156
612
102
64
Time
77.1
91.3
96.6
64.8
56.8
65.1
77.3
80.6
66.9
81.3
S500k
Iter.
4328
2165
394
107
4179
2359
390
185
Time
—
5366
5232
4757
2646
—
5250
5715
4629
4637
COCR
T5k
Iter.
124
32
26
25
25
109
32
32
32
32
Time
0.9
1.0
1.6
7.5
14.8
0.8
1.0
1.9
9.5
18.1
T16k
Iter.
31
24
23
23
23
31
29
28
29
30
Time
0.6
1.9
3.6
16.6
32.9
0.6
2.2
4.2
21.0
43.2
T80k
Iter.
3591
533
226
52
32
2947
413
223
49
34
Time
209.5
132.1
109.3
123.4
151.1
172.0
102.6
107.7
116.1
160.2
T150k
Iter.
303
79
42
19
19
267
62
36
23
23
Time
14.2
15.1
15.5
34.0
67.4
13.0
12.4
13.3
41.1
81.9
T500k
Iter.
4066
828
407
80
41
4067
807
397
88
46
Time
1410
1204
1182
1208
1218
1471
1172
1174
1330
1357
S13k
Iter.
3683
896
481
64
38
3153
742
452
93
53
Time
24.7
23.8
24.2
16.1
20.5
20.8
20.3
20.6
24.2
27.8
S32k
Iter.
3904
918
487
76
43
3352
756
408
90
61
Time
64.1
62.6
65.27
48.1
58.0
54.6
52.8
56.9
57.2
74.6
S500k
Iter.
8504
2025
1017
191
75
7203
1634
856
216
102
Time
2567
2525
2473
2386
1909
2393
2058
2083
2572
2546Table 4
Name
Number of rows
NNZ
s1rmq4m1
5 489
262 411
s3rmt3m3
5 357
207 123
Pres_Poisson
14 822
715 804
Dubcova1
16 129
253 009
nasasrb
54 870
2 677 324
apache1
80 800
542 184
denormal
89 400
726 674
G2_circuit
150 102
726 674
pwtk
217 918
11 524 432
parabolic_fem
525 825
3 674 625These matrices are used paired where the first matrix of the pair corresponds to the real part and the second to the imaginary part. If the dimensions disagree, the smallest matrix is padded with zero rows and columns up to the size of the biggest one. The pairing of the matrices with the name of the corresponding test is resumed in Table5.Table 5
Test name
Matrix pairing
Number of rows
NNZ
T5k
(s1rmq4m1, s3rmt3m3)
5 489
399 907
T16k
(Pres_Poisson, Dubcova1)
16 129
925 819
T80k
(apache1, nasasrb)
80 800
3 072 500
T150k
(G2_circuit, denormal)
150 102
1 616 970
T500k
(pwtk, parabolic_fem)
525 825
14 810 591The right hand side used for all tests, unless explicitly written, is assumed to be(
1
+
i
)
1, that is, 1
+
i for all components. A test from a real application is found in [2], from which the complex symmetric SPD matrices are provided with a specific right hand side. Size and name of the matrices are resumed in Table 6.Table 6
Test name
Matrix name
Number of rows
NNZ
S13k
eddy_6k_gauged
13 067
295 571
S32k
eddy_16k_gauged
31 853
720 689
S500k
eddy_265k_gauged
538 709
11 690 125List of the tests is as follows: for the first group (T tests) the right hand side was(
1
+
i
)
1; for the second group (S tests) the right hand side was the one prescribed in the paper [2].


### 7.1. Analysis of the Experiments

From the comparison of Tables1, 2, and 3 it is clear that there is no preconditioner that is always better than the others. The optimal preconditioner is in general problem dependent; therefore, for general preconditioners, a good quality measure is the robustness over different problems. Table 3 shows that the strategy presented in Section 2 is effective. In fact, when the ILU factorization is available and the iterations converge, incomplete factorization preconditioner is faster than polynomial preconditioner, but the proposed polynomial preconditioner is an effective alternative when ILU is not available or when it is not sufficient as preconditioner. This is true also for the relatively small matrices with 5000 rows, where the number of iterations of the ILU-based methods is comparable with the polynomial preconditioners of low degree, for example, 10. Rising the degree of the polynomial corresponds to lowering the number of iterations needed by COCG/COCR. It is also apparent that it is not possible to go below a certain number of iterations even with a very high degree polynomial; this is evident, for example, in test T16k of Table 3 with both COCG and COCR and with both preconditioners. This is explained from the fact that the condition number of the preconditioned system is independent of the system size. Another behaviour that is common to most of the tests is the better performance of the COCR over the COCG: this can be appreciated looking at Figure 2 and Figure 3. They show the history of the residual for each iteration of both methods with the Jacobi and the Chebyshev preconditioners. In both cases the decrease of the residual is always slower in the COCG than in the COCR; also the number of iterations is higher.Figure 2
The history of the residual for the test T80k with COCR and COCG preconditioned with Jacobi and Chebyshev polynomial.Figure 3
The history of the residual for the test S500k with COCR and COCG preconditioned with Jacobi and Chebyshev polynomial.


## 7.1. Analysis of the Experiments

From the comparison of Tables1, 2, and 3 it is clear that there is no preconditioner that is always better than the others. The optimal preconditioner is in general problem dependent; therefore, for general preconditioners, a good quality measure is the robustness over different problems. Table 3 shows that the strategy presented in Section 2 is effective. In fact, when the ILU factorization is available and the iterations converge, incomplete factorization preconditioner is faster than polynomial preconditioner, but the proposed polynomial preconditioner is an effective alternative when ILU is not available or when it is not sufficient as preconditioner. This is true also for the relatively small matrices with 5000 rows, where the number of iterations of the ILU-based methods is comparable with the polynomial preconditioners of low degree, for example, 10. Rising the degree of the polynomial corresponds to lowering the number of iterations needed by COCG/COCR. It is also apparent that it is not possible to go below a certain number of iterations even with a very high degree polynomial; this is evident, for example, in test T16k of Table 3 with both COCG and COCR and with both preconditioners. This is explained from the fact that the condition number of the preconditioned system is independent of the system size. Another behaviour that is common to most of the tests is the better performance of the COCR over the COCG: this can be appreciated looking at Figure 2 and Figure 3. They show the history of the residual for each iteration of both methods with the Jacobi and the Chebyshev preconditioners. In both cases the decrease of the residual is always slower in the COCG than in the COCR; also the number of iterations is higher.Figure 2
The history of the residual for the test T80k with COCR and COCG preconditioned with Jacobi and Chebyshev polynomial.Figure 3
The history of the residual for the test S500k with COCR and COCG preconditioned with Jacobi and Chebyshev polynomial.


## 8. Conclusions

A polynomial preconditioner was presented for the solution of the linear systemA
x
=
b, for A complex symmetric, that is, such that A
=
B
+
i
C, where B, C are real symmetric positive semidefinite matrices (semi-SPD) and B
+
C is symmetric positive definite (SPD). Typical problems of this form come from the field of electrodynamics, where the involved matrices are complex but not Hermitian and standard methods cannot be used directly. This algorithm is suitable for large matrices, where Cholesky decomposition, or its inexact form, are too costly or infeasible. It works as a polynomial approximation of a single step of the MHSS method, but it is successfully applied as preconditioner of conjugate gradient-like methods; in particular it is showed how to use it together with COCG or COCR. Following the trend of the last years, but being aware of the criticism that arose in the 80s, the proposed new preconditioner is computed as a recurrence of orthogonal polynomials and is proved to be stable. This allows employing polynomials of very high degree and numerical tests confirm the expected theoretical good performances.



---

*Source: 548609-2015-03-19.xml*