A New Block Krylov Subspace Method for Computing HighAccuracy Solutions
Hiroto TADANO1) , Tetsuya SAKURAI1) , and Yoshinobu KURAMASHI2)
1) Department of Computer Science, University of Tsukuba, 1-1-1 Tennodai, Tsukuba,
Ibaraki 305-8573, Japan
2) Graduate School of Pure and Applied Sciences, University of Tsukuba, 1-1-1 Tennodai,
Tsukuba, Ibaraki 305-8571, Japan
E-mail: {tadano, sakurai}@cs.tsukuba.ac.jp,
kuramasi@het.ph.tsukuba.ac.jp
Abstract
In this paper, the inﬂuence of errors which arise in matrix multiplications on theaccuracy of approximate solutions generated by the Block BiCGSTAB method is analyzed. In order to generate high accuracy solutions, a new Block Krylov subspacemethod is also proposed. Some numerical experiments illustrate that high accuracy
solutions can be obtained by using the proposed method compared with the Block
BiCGSTAB method.
Keywords: Block Krylov subspace methods, Block BiCGSTAB, linear systems withmultiple right hand sides, high accuracy solutions

1

Introduction

Linear systems with multiple right hand sides
AX = B,(1)

where A ∈ Cn×n , X, B ∈ Cn×L , appear in many scientiﬁc applications such as lattice quantum chromodynamics (lattice QCD) calculation of physical quantities [1], an eigensolverusing contour integration [6]. To solve these linear systems, some Block Krylov subspacemethods (e.g. Block BiCG [5], Block BiCGSTAB [3], Block QMR [2]) have been proposed.
Block Krylov subspace methods can compute approximate solutions of linear systems
with multiple right hand sides efﬁciently compared with Krylov subspace methods for single
right hand side [2]. However, the difference between the residual generated by the recursionof the Block BiCGSTAB method and the true residual may arise. In this paper, this difference which arises in the Block BiCGSTAB method is analyzed. Then, a new Block Krylov
subspace method for computing high accuracy approximate solutions of (1) is also proposed.
This paper is organized as follows. In Section 2, a matrix-valued polynomial and an
operation are deﬁned. The Block BiCGSTAB method is brieﬂy described in Section 3. In
Section 4, the inﬂuence of errors which arise in matrix multiplications on the accuracy ofapproximate solutions of the Block BiCGSTAB method. In Section 5, a new Block Krylov
subspace method is proposed for computing high accuracy approximate solutions. Thenthe true residual of the proposed method is also evaluated. In Section 6, the true residual
norm of the Block BiCGSTAB method and the proposed method is veriﬁed by numericalexperiments. The paper is concluded in Section 7.
1

2

Matrix-valued polynomial

Let Mk (z) be a matrix-valued polynomial of degree k deﬁned by
Mk (z) ≡k∑z j Mj ,j=0

where Mj ∈ CL×L and z ∈ C. The operation ◦ is used in this paper for the multiplication
Mk (A) ◦ V ≡k∑

Aj V M j ,j=0

where V ∈ Cn×L . This operation satisﬁes the following properties [3].
Proposition 2.1. Let M(z) and N (z) be matrix-valued polynomials of degree k and let Vand ξ be an n × L matrix and an L × L matrix, respectively. Then, the following properties
are satisﬁed.
(1) (M(A) ◦ V )ξ = (Mξ)(A) ◦ V ,(2) (M + N )(A) ◦ V = M(A) ◦ V + N (A) ◦ V .

3

The Block BiCGSTAB method

The (k+1)-th residual Rk+1 ∈ Cn×L of the Block BiCGSTAB method is deﬁned by
Rk+1 = B − AXk+1 ≡ (Qk+1 Rk+1 )(A) ◦ R0 ,(2)

where R0 = B − AX0 is an initial residual. The matrix-valued polynomial Rk+1 (z) ofdegree (k+1) which appears in (2) can be computed by the following recursions
R0 (z) = P0 (z) = IL ,Rk+1 (z) = Rk (z) − zPk (z)αk ,Pk+1 (z) = Rk+1 (z) + Pk (z)βk ,where Pk+1 (z) is an auxiliary matrix-valued polynomial of degree (k+1), IL is an L × L
identity matrix, αk and βk are L × L complex matrices. The polynomial Qk+1 (z) of degree(k+1) can be computed by
Q0 (z) = 1,
Qk+1 (z) = (1 − ζk z)Qk (z),where ζk ∈ C. The residual Rk+1 can be computed by the following recursions,
Rk+1 = Tk − ζk ATk ,Pk+1 = Rk+1 + (Pk − ζk APk )βk ,Tk = Rk − APk αk ,

2

(3)
(4)

X0 ∈ Cn×L is an initial guess,
Compute R0 = B − AX0 ,Set P0 = R0 ,Choose R̃0 ∈ Cn×L ,For k = 0, 1, . . . , until ∥Rk ∥F ≤ ε∥B∥F do:
Vk = APk ,Solve (R̃0H Vk )αk = R̃0H Rk for αk ,Tk = Rk − Vk αk ,Zk = ATk ,[][]ζk = Tr ZkH Tk /Tr ZkH Zk ,Xk+1 = Xk + Pk αk + ζk Tk ,Rk+1 = Tk − ζk Zk ,Solve (R̃0H Vk )βk = −R̃0H Zk for βk ,Pk+1 = Rk+1 + (Pk − ζk Vk )βk ,End
Figure 1: Algorithm of the Block BiCGSTAB method.
where matrices Pk+1 and Tk are deﬁned by Pk+1 ≡ (Qk+1 Pk+1 )(A)◦R0 and Tk ≡ (Qk Rk+1 )(A)◦
R0 , respectively. From Eqs. (2), (3), and (4), recursion for the approximate solution Xk+1
can be obtained by
Xk+1 = Xk + Pk αk + ζk Tk .
(5)
The L × L matrices αk and βk are determined so that bi-orthogonal conditions:R̃0H Aj (Rk (A) ◦ R0 ) = OL , j = 0, 1, . . . , k−1,
R̃0H Aj+1 (Pk (A) ◦ R0 ) = OL , j = 0, 1, . . . , k−1,

(6)
(7)

are satisﬁed. Here, R̃0 is an n × L nonzero matrix, OL is an L × L zero matrix, and ∥ · ∥Fdenotes the Frobenius norm of a matrix. The scalar parameter ζk is determined so that
∥Rk+1 ∥F is minimized. Fig. 1 shows the algorithm of the Block BiCGSTAB method. Here,Tr[ · ] denotes the trace of a matrix, and ε > 0 is a sufﬁciently small value for the stopping
criterion.

4

Evaluation of the true residual of the Block BiCGSTAB
methodIn this section, it is assumed that computation errors arise in the multiplications with α0 , α1 , . . . , αkwhich appear in the Block BiCGSTAB method. The inﬂuence of these errors on the true
residual of the Block BiCGSTAB method is considered. A matrix enclosed by a symbol〈 · 〉 denotes the perturbed matrix. Throughout this section, it is assumed that no calculation
errors arise except for multiplications with α0 , α1 , . . . , αk .

3The perturbed matrices 〈Pj αj 〉 and 〈(APj )αj 〉 are required for the computation of Xj+1
and Rj+1 , respectively. These matrices can be written as follows:〈Pj αj 〉 = Pj αj + Fj ,〈(APj )αj 〉 = APj αj + Gj ,

(8)
(9)

where Fj and Gj denote error matrices.
From Eqs. (5) and (8), Xk+1 is written as
Xk+1 = Xk + 〈Pk αk 〉 + ζk Tk= Xk + Pk αk + ζk Tk + Fkkk∑∑
= X0 +(Pj αj + ζj Tj ) +
Fj .
j=0

(10)

j=0

By using the Eq. (9), the residual Rk+1 generated by the recursion (3) is also written as
Rk+1 = Rk − 〈(APk )αk 〉 − ζk ATk= Rk − APk αk − ζk ATk − Gkkk∑∑
= R0 −(APj αj + ζj ATj ) −
Gj .
j=0

(11)

j=0

By using Eqs. (10) and (11), the true residual B − AXk+1 of the Block BiCGSTAB
method is given by
B − AXk+1 = R0 −k∑(APj αj + ζj ATj ) −j=0= Rk+1 +

k
∑

AFj

j=0

k
∑

Ej ,

(12)

j=0

where the matrix Ej is deﬁned by Ej ≡ Gj − AFj . From Eqs. (8) and (9), the matrix Ej canbe represented by
Ej = 〈(APj )αj 〉 − A〈Pj αj 〉.
The error matrices E0 , E1 , . . . , Ek appear in (12) when the computation errors occur in themultiplications with α0 , α1 , . . . , αk . The Eq. (12)∑implies that the true residual B − AXk+1of the Block BiCGSTAB method approaches to kj=0 Ej when the residual norm ∥Rk+1 ∥F
is sufﬁciently small.

5

Derivation of a new Block Krylov subspace method

The matrices Fj and Gj give adverse effects to the true residual of Block BiCGSTAB. To
negate the inﬂuence of these matrices, a condition Gj = AFj should be satisﬁed. In this section, a new Block Krylov subspace method is proposed to compute high accuracy solutions.

45.1 Construction of an algorithm
There are two ways of constructing the recursion for the residual Rk+1 = (Qk+1 Rk+1 )(A) ◦R0 . In the Block BiCGSTAB method, ﬁrstly, the polynomial Qk+1 (z) is expanded. In thiscase, as shown in the Eq. (12), the true residual B − AXk+1 is not equal to the residual
Rk+1 generated by the recursion. In the proposed method, ﬁrstly, the polynomial Rk+1 is
expanded for computing Qk+1 Rk+1 . The recursion of this polynomial is given by
Qk+1 Rk+1 = Qk Rk − ζk zQk Rk − zQk+1 Pk αk .
The polynomials Qk+1 Rk and Qk+1 Pk+1 are computed by the following recursions:
Qk+1 Pk = Qk Pk − ζk zQk Pk ,Qk+1 Pk+1 = Qk+1 Rk+1 + Qk+1 Pk βk .
From the above recursions, the residual Rk+1 and auxiliary matrices can be computed by
Rk+1 = Rk − ζk ARk − AUk ,Pk+1 = Rk+1 + Uk αk−1 βk ,Sk = Pk − ζk APk ,

(13)
(14)

where Sk ≡ (Qk+1 Pk )(A) ◦ R0 and Uk ≡ Sk αk . By Eqs. (2) and (13), the recursion for
Xk+1 can be obtained by
Xk+1 = Xk + ζk Rk + Uk .
(15)
Matrices αk and βk are determined so that bi-orthogonality conditions (6) and (7) aresatisﬁed. From the Eq. (6), the matrix αk can be computed by
αk = (R̃0H APk )−1 R̃0H Rk .

(16)

By the bi-orthogonality condition (7) and the relation
R̃0H Rk+1 = −ζk R̃0H ATk ,the matrix βk can be obtained by
βk = (R̃0H APk )−1 R̃0H Rk+1 /ζk .

(17)

The matrix γk ≡ αk−1 βk is appeared in the Eq. (14). By using Eqs. (16) and (17), γk can beobtained by
γk = (R̃0H Rk )−1 R̃0H Rk+1 /ζk .
If the parameter ζk is determined so that ∥Rk+1 ∥F is minimized, then extra multiplications
with A are required in the proposed method. To avoid the multiplications with A, the parameter ζk is computed by[][]ζk = Tr (ARk )H Rk /Tr (ARk )H ARk .
In the proposed method, the three multiplications with A are required in each iteration. To
reduce the number of multiplications with A, the matrix APk+1 is computed by the following
recursion
APk+1 = ARk+1 + AUk γk .
The proposed method can be summarized in Fig. 2.
5

X0 ∈ Cn×L is an initial guess,
Compute R0 = B − AX0 ,Set P0 = R0 and V0 = W0 = AR0 ,Choose R̃0 ∈ Cn×L ,For k = 0, 1, . . . , until ∥Rk ∥F ≤ ε∥B∥F do:
Solve (R̃0H Vk )αk = R̃0H Rk for αk ,[][]
ζk = Tr WkH Rk /Tr WkH Wk ,Sk = Pk − ζk Vk ,Uk = Sk αk ,Yk = AUk ,Xk+1 = Xk + ζk Rk + Uk ,Rk+1 = Rk − ζk Wk − Yk ,Wk+1 = ARk+1 ,Solve (R̃0H Rk )γk = R̃0H Rk+1 /ζk for γk ,Pk+1 = Rk+1 + Uk γk ,Vk+1 = Wk+1 + Yk γk ,EndFigure 2: Algorithm of the proposed method.

5.2 Evaluation of the true residual
Similar to the previous section, assume that no calculation errors arise except for the multiplications with α0 , α1 , . . . , αk . The multiplication with αj is appeared in the computation of
Uj = Sj αj . By using the symbol 〈 〉, the perturbed matrix 〈Sj αj 〉 is represented by〈Sj αj 〉 = Sj αj + Hj ,(18)

where Hj is an error matrix. From Eqs. (15) and (18), the approximate solution Xk+1 iswritten as
Xk+1 = Xk + ζk Rk + 〈Sk αk 〉= Xk + ζk Rk + Sk αk + Hkkk∑∑
= X0 +(ζj Rj + Sj αj ) +
Hj .
j=0

j=0

By using Eqs. (13) and (19), Rk+1 is represented by
Rk+1 = Rk − ζk ARk − A〈Sk αk 〉= Rk − ζk ARk − A[Sk αk + Hk ]kk∑∑
= R0 −(ζj ARj + ASj αj ) −
AHj
j=0

j=0

6

(19)

[
= B − A X0 +k∑

(ζj Rj + Sj αj ) +

j=0

k
∑

]
Hj

j=0

= B − AXk+1 .
By regarding the matrices Hj and AHj as Fj and Gj , it is conﬁrmed that the proposed
method satisﬁes the condition Ej = Gj − AFj = O.

6

Numerical experiments

Test matrices used in numerical experiments were PDE900, JPWH991, and CONF5.4-00L8X81000 [4]. The size and the number of nonzero elements of these matrices are shown in Table
1. The coefﬁcient matrix of CONF5.4-00L8X8-1000 is constructed by In − κD, where D isan n × n non-Hermitian matrix and κ is a real valued parameter. This parameter was set to0.1782.
The initial solution X0 was set to the zero matrix. The shadow residual R̃0 was given by
a random number generator. The right hand side B of (1) was given by B = [e1 , e2 , . . . , eL ],where ej is a j-th unit vector. The value ε for the stopping criterion was set to 1.0 × 10−14 .
All experiments were carried out in double precision arithmetic on CPU: Intel Core
2 Duo 2.4GHz, Memory: 4GBytes, Compiler: Intel Fortran ver. 10.1, Compile option:-O3 -xT -openmp. The multiplication with the coefﬁcient matrix was parallelized by
OpenMP.
The results of the Block BiCGSTAB method are shown in Table 2. In this Table, #Iter.,Res, and True Res. denote the number of iterations, the relative residual norm ∥Rk ∥F /∥B∥F ,and the true relative residual norm ∥B − AXk ∥F /∥B∥F , respectively.
As shown in Table 2, the relative residual norms of the Block BiCGSTAB method weresatisﬁed the convergence criterion. However, the true residual norms did not reach 10−14
when L = 2, 4.
∑
The relation between the true relative residual norm and ∥ kj=0 Ej ∥F /∥B∥F for JPWH991
with L ∑= 4 is shown in Fig. 3. The true relative residual norm became almost equal to thevalue ∥ kj=0 Ej ∥F /∥B∥F . The Eq. (12) was veriﬁed through this numerical example.
The results of the proposed method are shown in Table 3. The true relative residual
norms computed by the proposed method reached 10−14 except for JPWH991 with L = 1.
By using the proposed method, high accuracy solutions can be obtained compared with theBlock BiCGSTAB method.
Table 1: The size and the number of nonzero elements of test matrices. NNZ denotes thenumber of nonzero elements.
Matrix name
Size
NNZ
PDE900900
4,380
JPWH991
991
6,027
CONF5.4-00L8X8-1000 49,152 1,916,928

7

Table 2: Results of the Block BiCGSTAB method.
L #Iter.
1
53
2
46
4
41
L #Iter.
1
56
2
49
4
43
L #Iter.
1
555
2
452
406
4

PDE900Time/L [s]Res.
True Res.
−15
0.0096 4.8 × 10
4.8 × 10−15−150.0067 1.1 × 10
2.0 × 10−130.0031 4.8 × 10−15 1.8 × 10−12
JPWH991
Time/L [s]Res.
True Res.
0.0159 5.7 × 10−15 1.2 × 10−140.0083 8.3 × 10−15 4.1 × 10−130.0034 6.3 × 10−15 5.9 × 10−12
CONF5.4-00L8X8-1000
Time/L [s]Res.
True Res.
−15
13.9408 8.9 × 10
9.5 × 10−15−15
7.5609 7.3 × 10
2.5 × 10−13
6.1544 8.7 × 10−15 2.8 × 10−13

Table 3: Results of the proposed method.

L #Iter.
1
53
2
46
4
45
L #Iter.
1
52
2
51
4
44
L #Iter.
1
555
2
456
4
386

PDE900Time/L [s]Res.
True Res.
0.0107 3.2 × 10−15 3.3 × 10−15
0.0051 1.1 × 10−15 1.4 × 10−15
0.0031 5.5 × 10−15 5.6 × 10−15
JPWH991
Time/L [s]Res.
True Res.
−15
0.0134 8.4 × 10
1.3 × 10−14−150.0082 3.7 × 10
6.1 × 10−15
0.0035 1.5 × 10−15 2.3 × 10−15
CONF5.4-00L8X8-1000
Time/L [s]Res.
True Res.
14.2714 7.4 × 10−15 8.5 × 10−15
8.1093 5.6 × 10−15 6.7 × 10−15
6.0348 7.4 × 10−15 8.6 × 10−158Figure∑k 3: Relation between the true relative residual norm of Block BiCGSTAB and∥ j=0 Ej ∥F /∥B∥F (JPWH991, L = 4). — : true relative residual norm ∥B −∑AXk ∥F /∥B∥F , — : relative residual norm ∥Rk ∥F /∥B∥F , — : ∥ kj=0 Ej ∥F /∥B∥F .

7

Conclusions

In this paper, we have evaluated the true residual of the Block BiCGSTAB method when
the computation errors arise in the multiplications with α0 , α1 , . . . , αk . We have shown thatthe true residual of the Block BiCGSTAB method approaches to the sum of error matrices
when the residual norm is sufﬁciently small. Then, we have proposed the new method forgenerating high accuracy approximate solutions. Through some numerical experiments, wehave veriﬁed that the proposed method generates the high accuracy solutions compared with
the Block BiCGSTAB method.

Acknowledgments
This work was supported by Grant-in-Aid for Young Scientists (Start-up) (No. 20800009).

References[1] PACS-CS Collaboration, S. Aoki et al., 2 + 1 Flavor Lattice QCD toward the PhysicalPoint, arXiv:0807.1661v1 [hep-lat].
[2] R. W. Freund, A Block QMR Algorithm for Non-Hermitian Linear Systems with Multiple Right-Hand Sides, Lin. Alg. Appl., Vol. 254, pp. 119–157, 1997.
[3] A. El Guennouni, K. Jbilou, and H. Sadok, A Block Version of BiCGSTAB for Linear
Systems with Multiple Right-Hand Sides, Elec. Trans. Numer. Anal., Vol. 16, pp. 129–
142, 2003.
[4] Matrix Market, http://math.nist.gov/MatrixMarket/[5] D. P. O’Leary, The Block Conjugate Gradient Algorithm and Related Methods, Lin.
Alg. Appl., Vol. 29, pp. 293–322, 1980.

9

[6] T. Sakurai, H. Tadano, T. Ikegami, and U. Nagashima, A Parallel Eigensolver Using
Contour Integration for Generalized Eigenvalue Problems in Molecular Simulation,
Tech. Rep. CS–TR–08–14, University of Tsukuba, 2008.

10