# MINRES Seed Projection Methods for Solving Symmetric Linear Systems with Multiple Right-Hand Sides

Authors: Xin Li; Hao Liu; Jingfu Zhu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We consider the MINRES seed projection method for solving multiple right-hand side linear systemsA
X
=
B, where A
∈
R
n
×
n is a nonsingular symmetric matrix, B
∈
R
n
×
p. In general, GMRES seed projection method is one of the effective methods for solving multiple right-hand side linear systems. However, when the coefficient matrix is symmetric, the efficiency of this method would be weak. MINRES seed projection method for solving symmetric systems with multiple right-hand sides is proposed in this paper, and the residual estimation is analyzed. The numerical examples show the efficiency of this method.


---

## Body

## 1. Introduction

Consider the multiple right-hand side linear systems(1)
A
X
=
B
,
where A
∈
R
n
×
n is a nonsingular symmetric matrix, and B
∈
R
n
×
p.Equation (1) plays an important role in chemistry, electronics, structures, control, and other problems; see [1, 2] for detail.In the last few years, the block methods have been developed to solve (1), such as the block conjugate gradient algorithm (BCG) [3], the block generalized minimal residual method [4–6], the block BiCGSTAB method [7], the block QMR method [8], the block least squares method [9], the block Lanczos method [10], and the block IDR(s) method [11] which have been proposed recently. In general, the block methods are faster than solving each one separately.The meshless methods are extensively used for solving (1); these meshless methods show to be efficient and accurate in terms of their numerical results; see [12–18]. The global methods [19, 20] are also a class of important methods. Following the work [20], many other global methods have been developed, including the global BiCG and global BiCGSTAB methods [21], the global Hessenberg and global CMRH methods [22], and the polynomial preconditioned global CMRH method [23, 24]. Generally, the global methods are more appropriate for large and sparse systems.In many practical applications, the right-hand sides are not arbitrary and are very close; then the seed projection methods are often used to solve (1); see [1, 2, 25, 26]. The main idea of this method is selecting one system to be the seed systems firstly, then solving the seed systems by some Krylov subspace method, and creating a Krylov subspace K, then projecting the residual of the other systems, called nonseed systems, onto this Krylov subspace K to get the approximate solutions. The process is repeated with other seed systems until all the systems are solved. The seed projection methods were proposed by Smith et al. [1] for the CG method firstly. When A is unsymmetric, Simoncini and Gallopulos [2] proposed the GMRES seed projection method. Later, a seed method which uses Morgan's Krylov subspace augmented with eigenvectors was presented in [25]. Moreover, the seed method can be used to solve unsymmetric shifted systems with multiple right-hand sides [26]. However, if A is a symmetric matrix, the efficiency of these methods would be weak.In this paper, we propose the MINRES seed projection method for solving symmetric systems with multiple right-hand sides, and the residual estimation is analyzed.The paper is organized as follows. In Section2, we give a quick overview of the GMRES seed projection method. In Section 3, we present the MINRES seed projection method and the residual estimation. In the last section, we show the efficiency of our method by numerical experiments.


## 2. GMRES Seed Projection Method

In this section, we recall the GMRES seed projection method for solving (1). Details of the algorithm can be found in [2, 25]. We summarize it in the following algorithm.Algorithm 1 (X
=
Gseed
(
A
,
B
,
X
(
0
)
,
ε
,
m
)).
We have the following:(1)
X
=
X
(
0
);
(2)
R
=
B
-
A
X;
(3)
forl
=
1,2
,
…
,
p until all the systems are solved
(4)
[
σ
,
r
σ
]
=
S
E
E
D
(
R
),  β
=
∥
r
σ
∥
2;
(5)
fork
=
1,2
,
…, until convergence;
(6)
[
V
m
+
1
,
H
]
=
A
r
n
o
l
d
i
(
A
,
r
σ
);
(7)
b
σ
^
=
β
e
1, where e
1
=
(
1,0
,
…
,
0
)
T
∈
R
m
+
1;
(8)
b
j
^
=
V
m
+
1
Τ
(
b
j
-
A
x
j
),  j
=
1
,
…
,
p,  j
≠
σ;
(9)
computey
j by minimizing ∥
H
y
-
b
j
^
∥
2,  y
∈
R
m,  j
=
1
,
…
,
p;
(10)
X
-
=
X
+
V
m
Y, where Y
=
[
y
1
,
…
,
y
p
]
Τ;
(11)
R
-
=
B
-
A
X
-,  r
-
(
j
)
=
∥
R
-
(
:
,
j
)
∥
2;
(12)
if∥
r
-
(
j
)
∥
2
<
ε, then delete the jth systems, let w be the number of the j conforming to this condition, and then p
:
=
p
-
w;
(13)
if∥
r
-
σ
∥
2
<
ε, then delete these seed systems, set p
:
=
p
-
1, and go to 3;
(14)
end(k);
(15)
end(l).We now make a few descriptions of Algorithm1. Firstly, we give an initial approximation to the solutions X
(
0
), compute the initial residuals, and select seed systems by a function S
E
E
D. In this algorithm, S
E
E
D applied to the R returns σ and r
σ, where σ is the index of column of R having the maximum norm. Secondly, we apply restarted GMRES method for solving the seed systems and function A
r
n
o
l
d
i applies the Arnoldi procedure to generate an orthogonal basis V
m
+
1
=
[
v
1
,
…
,
v
m
] for the Krylov subspace K
m
+
1
(
A
,
r
σ
). Meanwhile, the nonseed solutions are approximated by projecting the residual r
j
=
b
j
-
A
x
j on K
m
+
1
(
A
,
r
σ
) and solving the least square problem min
⁡
y
j
∈
R
m
⁡
∥
H
y
j
-
V
m
+
1
Τ
r
j
∥
2. Thirdly, after the seed systems are solved to desired accuracy, new seed systems are selected from the unsolved systems and then the whole procedure is repeated until all the systems are solved. Some theoretical analysis about the above algorithm can be found in [2, 25].


## 3. MINRES Seed Projection Method

Based on the MINRES method for solving the symmetric linear systems, in this section, we combine the GMRES seed projection method and the MINRES method and propose the MINRES seed projection method for solving (1). The Arnoldi procedure in GMRES seed projection method is exchanged by Lanczos procedure, and applying the seed projection idea to MINRES method, the MINRES seed projection method is proposed as the following algorithm.Algorithm 2 (X
=
Mseed
(
A
,
B
,
X
(
0
)
,
ε
,
m
)).
We have the following:(1)
X
=
X
(
0
);
(2)
R
=
B
-
A
X;
(3)
forl
=
1,2
,
…
,
p until all the systems are solved;
(4)
[
σ
,
r
σ
]
=
S
E
E
D
(
R
),  β
=
∥
r
σ
∥
2;
(5)
fork
=
1,2
,
…, until convergence;
(6)
[
Q
m
+
1
,
T
~
m
]
=
L
a
n
c
z
o
s
(
A
,
r
σ
);
(7)
b
σ
^
=
β
e
1, where e
1
=
(
1,0
,
…
,
0
)
T
∈
R
m
+
1;
(8)
b
j
^
=
V
m
+
1
Τ
(
b
(
j
)
-
A
x
(
j
)
),  j
=
1
,
…
,
p,  j
≠
σ;
(9)
computed
(
j
) by minimizing ∥
T
~
m
d
-
b
j
^
∥
2,  d
∈
R
m,  j
=
1
,
…
,
p;
(10)
X
~
=
X
+
Q
m
D, where D
=
[
d
(
1
)
,
…
,
d
(
p
)
]
Τ;
(11)
R
~
=
B
-
A
X
~,  r
~
(
j
)
=
∥
R
~
(
:
,
j
)
∥
2;
(12)
if∥
r
~
(
j
)
∥
2
<
ε, then delete jth system, let w be the number of the j conforming to this condition, and then p
:
=
p
-
w;
(13)
if∥
r
~
σ
∥
2
<
ε, then delete this seed systems, set p
:
=
p
-
1, go to 3;
(14)
end(k);
(15)
end(l).We now make a few descriptions about Algorithm2. Firstly, we give an initial approximation to the solutions X
(
0
), compute the initial residuals, and select seed systems by a function S
E
E
D. In this algorithm, S
E
E
D applied to the R returns σ and r
σ, where σ is the index of column of R having the maximum norm. Secondly, we apply restarted Lanczos method for solving the seed systems and function L
a
n
c
z
o
s applies the Lanczos procedure to generate an orthogonal basis Q
m
+
1
=
[
q
1
,
…
,
q
m
] for the Krylov subspace K
m
+
1
(
A
,
r
σ
). Meanwhile, the nonseed solutions are approximated by projecting the residual r
(
j
)
=
b
(
j
)
-
A
x
(
j
) on K
m
+
1
(
A
,
r
σ
) and solving the least square problem min
⁡
d
(
j
)
∈
R
m
⁡
∥
T
~
m
d
(
j
)
-
Q
m
+
1
Τ
r
(
j
)
∥
2. Thirdly, after the seed systems are solved to desired accuracy, new seed systems are selected from the unsolved systems and then the whole procedure is repeated until all the systems are solved.According to Algorithm2, we can get
(2)
r
~
(
j
)
=
b
(
j
)
-
A
x
~
(
j
)
=
b
(
j
)
-
A
(
x
(
j
)
+
Q
m
d
(
j
)
)
=
r
(
j
)
-
A
Q
m
d
(
j
)
=
r
(
j
)
-
Q
m
+
1
T
~
m
d
(
j
)
and Q
m
+
1
T
r
~
(
j
)
=
Q
m
+
1
T
r
(
j
)
-
T
~
m
d
(
j
). In addition, let P
=
Q
m
+
1
Q
m
+
1
T be the orthogonal operator on K
m
+
1
(
A
,
r
(
s
)
). It follows that
(3)
r
~
(
j
)
=
r
(
j
)
-
A
Q
m
d
(
j
)
=
(
I
-
P
)
r
(
j
)
+
(
P
r
(
j
)
-
A
Q
m
d
(
j
)
)
=
(
I
-
P
)
r
(
j
)
+
Q
m
+
1
(
Q
m
+
1
T
r
(
j
)
-
T
~
m
d
(
j
)
)
.
The following property can be attained.Property 1.
The residualr
~
(
j
) in nonseed systems in Algorithm 2 satisfy
(4)
∥
Q
m
+
1
T
r
~
(
j
)
∥
2
=
min
⁡
x
∈
x
(
j
)
+
K
m
∥
Q
m
+
1
T
(
b
(
j
)
-
A
x
)
∥
2
,
∥
r
~
(
j
)
∥
2
2
=
∥
(
I
-
P
)
r
(
j
)
∥
2
2
+
min
⁡
d
∈
R
m
∥
Q
m
+
1
T
r
(
j
)
-
T
~
m
d
∥
2
2
.Theorem 3.
SetK
m
≡
K
m
(
A
,
r
(
s
)
) as m-dimension Krylov subspace; one has(1)
for the seed systems:x
~
(
s
)
∈
{
x
(
s
)
}
+
K
m, r
~
(
s
)
∈
K
m
+
1, and r
~
(
s
)
⊥
A
K
m
⊆
K
m
+
1;
(2)
for the nonseed systems:x
~
(
j
)
∈
{
x
(
j
)
}
+
K
m, r
~
(
j
)
∈
{
r
(
j
)
}
+
K
m
+
1, and r
~
(
j
)
⊥
A
K
m.Proof.
(
1
) According to the Lanczos method, the conclusion of (
1
) is right obviously;
(
2
) since x
~
(
j
)
=
{
x
(
j
)
}
+
Q
m
d
(
j
)
∈
{
x
(
j
)
}
+
K
m, it follows that
(5)
r
~
(
j
)
=
b
(
j
)
-
A
x
~
(
j
)
=
r
(
j
)
-
A
Q
m
d
(
j
)
=
r
(
j
)
-
Q
m
+
1
T
~
m
d
(
j
)
.
Using A
Q
m as inner product, we can get
(6)
(
A
Q
m
)
T
r
~
(
j
)
=
T
~
m
T
Q
m
+
1
T
(
r
(
j
)
-
Q
m
+
1
T
~
m
d
(
j
)
)
=
T
~
m
T
(
Q
m
+
1
T
r
(
j
)
-
T
~
m
d
(
j
)
)
.
Since d
(
j
)
=
arg
min
⁡
d
∈
R
m
⁡
∥
Q
m
+
1
T
r
(
j
)
-
T
~
m
d
∥
2 is the least squares solution and satisfies the normal equation T
~
m
T
T
~
m
d
(
j
)
=
T
~
m
T
Q
m
+
1
T
r
(
j
), then (
A
Q
m
)
T
r
~
(
j
)
=
0.


## 4. Numerical Experiments

In order to prove the efficiency of our algorithm, we compare the following methods.(
1
) MINRES1, it uses MINRES method for solving the multiple right-hand side linear systems one by one and uses the convergent solution of (j
-
1)th systems as the iterative initial vector of jth systems when the (j
-
1)th systems converge. (
2
) MINRES2, it uses MINRES method for solving the multiple right-hand side linear systems one by one, and the initial vector is zero vector. (
3
) Mseed is MINRES seed projection method.All numerical experiments are implemented in MATLAB 2009 and run in Intel Pentium Dual T2390 computer. We setm
=
30; all the tests are stopped as soon as ∥
r
(
j
)
∥
2
<
ε
≡
1
0
-
6.Example 1.
A is a 1024
×
1024 symmetric matrix:
(7)
A
=
[
1
0.1
0.5
0.5
⋯
0.5
0.1
2
0.1
0.5
0.1
3
0.1
0.5
0.1
4
⋱
⋮
⋱
⋱
0.1
0.5
0.1
1024
]
.
In numerical experiment, the right-hand sides of systems are constructed by two forms as the following, respectively.Form 1:(8)
B
1
=
[
b
(
1
)
,
…
,
b
(
p
)
]
,
b
(
j
)
≡
b
(
j
)
(
t
i
)
=
-
cos
⁡
(
5
cos
⁡
(
t
i
-
2
(
j
-
1
)
π
128
)
)
,where j
=
1
,
…
,
p;  t
i
=
1
+
0.1
(
i
-
1
),  i
=
1
,
…
,
n.Form 2:(9)
B
2

=

[
b
(
1
)
,
…
,
b
(
p
)
]
,
where
b
(
j
)

=

A
u
(
j
)
,
u
(
1
)

=

(
1
,
…
,
1
)
T
,
u
(
j
)
=
j
·
cos
⁡
(
(
2
j
+
i
)
×
1
0
6
)
·
sin
(
(
3
(
4
-
j
)
+
i
)
×
1
0
(
6
)
)
,
i
=
1
,
…
,
n
,
j
=
2
,
…
,
p
.The numerical results are shown in Tables 1 and 2 and the data in bracket is the sum of iterative steps.Table 1
Iterative steps of convergent solution of each system.
Right-hand
Mseed
MINRES1
MINRES2
B
1
(81) 29, 16, 8, 9, 19
(155) 29, 20, 40, 45, 21
(286) 29, 41, 57, 73, 86
B
2
(57) 25, 15, 15, 2, 0
(125) 25, 27, 25, 21, 27
(181) 25, 36, 34, 46, 40Table 2
Convergent time (second) of each method.
Right-hand
Mseed
MINRES1
MINRES2
B
1
17.1720
19.2810
20.8430
B
2
17.0870
18.5620
19.7832In Table1, the iterative steps of each system are listed and the data in bracket are the sum of iterative steps of all systems. From Table 1, we can know that Mseed can do better than MINRES1 and MINRES2.Figures1, 2, and 3 are the convergent curve of three methods, respectively, when the right side is B
1. And from them, we can see that Mseed can do better than MINRES1 and MINRES2.Figure 1
MINRES seed projection method.Figure 2
MINRES1.Figure 3
MINRES2.Example 2.
Next example comes from Matrix Market, and all the matrices are symmetric as shown in Table3:
(10)
B
1
=
[
b
(
1
)
,
…
,
b
(
p
)
]
,
j
=
1
,
…
,
p
,
where t
i
=
1
+
0.1
(
i
-
1
),  i
=
1
,
…
,
n,  B
2
=
[
b
(
1
)
,
…
,
b
(
p
)
], and  b
(
j
)
=
(
1
,
…
,
1
,
0
,
…
,
0
︸
j
-
1
)
T.Table 3
Matrix name
Order
If diagonal dominance
Conditions number
bcsstm19
817
×
817
Y
2.3
e
+
05
bcsstk27
1224
×
1224
N
7.7
e
+
04
bcsstk22
138
×
138
N
1.7
e
+
05Letp
=
5, and the dimensions of projection space m
=
25 and results of calculation are shown in Table 4.Table 4
Convergent iterative steps of each system and CPU time (second).
B
A
Minres seed projection method
Minres2
Iterative steps
T1 (s)
Iterative steps
T2 (s)
B
1
bcsstm9
(718) 171, 146, 179, 129, 93
3.13
(980) 171, 201, 241, 191, 176
5.77
bcsstk27
(781) 211, 187, 154, 124, 105
11.82
(1033) 211, 206, 203, 203, 210
24.98
bcsstk22
(1833) 984, 574, 142, 119, 14
4.32
(4705) 984, 935, 850, 1014, 922
16.82
B
2
bcsstm9
(475) 175, 124, 62, 68, 46
2.24
(898) 175, 201, 175, 176, 171
5.37
bcsstk27
(748) 215, 156, 143, 140, 94
11.34
(1079) 215, 216, 216, 216, 216
25.85
bcsstk22
(2312) 1019, 590, 393, 162, 148
5.33
(5231) 1019, 986, 1079, 1059, 1088
18.85From Table4, we can see that Mseed is superior to MINRES method.Example 3.
A is a 1000
×
1000 symmetric matrix:
(11)
A
=
[
1
0.1
0.1
2
0.1
0.1
3
0.1
0.1
4
⋱
⋱
⋱
0.1
0.1
1024
]
.
Form 1:(12)
B
1
=
[
b
1
,
b
2
,
…
,
b
p
]
,
(
b
j
)
i
=
-
cos
⁡
(
5
cos
⁡
(
i
-
2
(
j
-
1
)
π
128
)
)
,where j
=
1
,
…
,
p;  t
i
=
1
+
0.1
(
i
-
1
),  i
=
1
,
…
,
n.Form 2:(13)
B
2

=

[
b
(
1
)
,
…
,
b
(
p
)
]
,
where
b
(
j
)

=

A
u
(
j
)
,
u
(
1
)

=

(
1
,
…
,
1
)
T
,
u
(
j
)
=
j
·
cos
⁡
(
(
2
j
+
i
)
×
1
0
3
)
·
sin
(
(
3
(
4
-
j
)
+
i
)
×
10
3
)
,
i
=
1
,
…
,
n
,
j
=
2
,
…
,
p
.The numerical results are shown in Tables 5 and 6 and the data in bracket is the sum of iterative steps.Table 5
Iterative steps of convergent solution of each system.
Right-hand
Mseed
MINRES1
MINRES 2
B
1
(34) 11, 8, 6, 4, 5
(48) 11, 10, 9, 9, 9
(54) 11, 11, 11, 10, 11
B
2
(36) 11, 13, 11, 0, 1
(60) 11, 11, 13, 12, 13
(60) 11, 11, 12, 13, 13Table 6
Convergent time (second) of each method.
Right-hand
Mseed
MINRES1
MINRES2
B
1
4.785072
5.967727
6.053866
B
2
4.934010
5.739067
5.823134


## 5. Conclusion

In this paper, we propose the MINRES seed projection method for solving symmetric systems with multiple right-hand sides, and the residual estimation is analyzed. The numerical examples show that our method is effective.



---

*Source: 357874-2014-06-12.xml*