# Deflated BiCG with an Application to Model Reduction

Authors: Jing Meng; Pei-Yong Zhu; Hou-Biao Li
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Most calculations in model reduction involve the solutions of a sequence of dual linear systems with multiple right-hand sides. To solve such systems efficiently, a new deflated BiCG method is explored in this paper. The proposed algorithm uses harmonic Ritz vectors to approximate left and right invariant subspaces inexpensively via small descenting direction vectors found by subsequent runs of deflated BiCG and then derives the deflated subspaces for the next pair of dual linear systems. This process leads to faster convergence for the next pair of systems. Numerical examples illustrate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Large scale simulations play an important role in the study of a great variety of complex physical phenomena, leading often to overwhelming demands on computational resources [1–5]. Hence, the common approach is to produce a surrogate model of much smaller dimension which provides a high-fidelity approximation of the original model. For such problems, interpolatory model reduction method combines flexibility and scalability and has proven effectiveness. It transfers function interpolations in the frequency domain to meet various desirable approximation goals. During this process, it requires the solutions of dual linear systems with multiple right-hand sides (RHSs):(1)
A
x
i
=
b
i
,
A
T
x
~
i
=
b
~
i
,
i
=
1,2
,
…
,
s
,where A
∈
R
n
×
n is a sparse, nonsymmetric matrix and RHSs are not available simultaneously.For the solutions of primary linear systemsA
x
i
=
b
i, i
=
1,2
,
…
,
s, deflated Krylov methods have been appearing. This is due to the fact that they take advantage of the fact that several systems share the same matrix. In addition, the convergence of Krylov subspace solvers for a linear system, to a great extent, depends on the spectrum of the matrix. If one could project the eigenvectors corresponding to the smallest eigenvalues out from the initial residual and then solve the deflated system it will converge much faster. The process is referred to as deflation. Variants of deflated Krylov solvers for the primary linear systems have been fully studied in the literature [2, 3, 6–14]. However, deflation for dual linear systems has not yet been fully investigated [15, 16].In this paper we extend this idea to the BiCG algorithm for the dual case. The goal of this paper is to develop a new deflated BiCG method for solving dual linear systems with multiple RHSs. Likewise, for BiCG, it can be shown that if the primal Krylov subspace is deflated with right eigenvectors, the corresponding left eigenvectors are removed from the dual residual. Therefore, while solving a pair of systems, we select approximate left and right invariant subspaces ofA and then use those to accelerate the convergence of the next pair of systems. The proposed algorithm uses harmonic Ritz vectors to approximate left and right invariant subspaces inexpensively via small descenting direction vectors found by subsequent runs of deflated BiCG and then derives the deflated subspaces for the next pair of dual linear systems. Furthermore, we describe a cheap way to build the deflated subspaces.In the next section, we describe the outline of the deflated Lanczos algorithm used in the derivation of the deflated BiCG method. In Section3, we derive a deflated BiCG method using previously computed deflated subspace matrices. How to compute and update such deflated space matrices is given in Section 4. In Section 5, we investigate the deflated BiCG algorithm for solving the dual systems with multiple RHSs. The effectiveness of the proposed method is also demonstrated in Section 6. Finally, some conclusions are summarized in Section 7.Throughout this paper,A
T is referred to as the transpose conjugate operation of matrix A, (
·
,
·
) denotes the inner product, 0 is defined as the zero matrix, and ⊥
b is biorthogonality.


## 2. The Deflated Lanczos Algorithm

In this section, we describe a deflated Lanczos algorithm that builds two sequences{
v
j
}
j
=
1,2
,
…,  {
v
~
j
}
j
=
1,2
,
… of vectors such that (2)
v
j
+
1
⊥
b
s
p
a
n
W
~
,
v
~
1
,
v
~
2
,
…
,
v
~
j
,
v
~
j
+
1
⊥
b
s
p
a
n
W
,
v
1
,
v
2
,
…
,
v
j
,and (
v
~
j
+
1
,
v
j
+
1
)
=
1, where W
=
[
w
1
,
w
2
,
…
,
w
k
] and W
~
=
[
w
~
1
,
w
~
2
,
…
,
w
~
k
] are two sets of linearly independent vectors.We assume that the matrixW
~
T
A
W is nonsingular. Let (3)
A
1
=
I
-
A
W
W
~
T
A
W
-
1
W
~
T
A
,
(4)
A
~
1
=
I
-
A
T
W
~
W
T
A
T
W
~
-
1
W
T
A
T
.Since the matrices A
1 and A
~
1 are each other’s conjugate transpose, we can apply the standard Lanczos procedure to the above auxiliary matrices. Let unit vectors v
1, v
~
1 be biorthogonal to W
~, W, respectively; then we compute the Lanczos vectors which satisfy (5)
A
1
V
j
=
V
j
T
j
+
t
j
+
1
,
j
v
j
+
1
e
j
T
,
(6)
A
~
1
V
~
j
=
V
~
j
T
~
j
+
t
~
j
+
1
,
j
v
~
j
+
1
e
j
T
,where V
j
=
[
v
1
,
v
2
,
…
,
v
j
], V
~
j
=
[
v
~
1
,
v
~
2
,
…
,
v
~
j
], T
j and T
~
j are j
×
j tridiagonal matrices, t
j
+
1
,
j, t
~
j
+
1
,
j are the last element of the last row of T
_
j, and T
~
_
j
∈
R
(
i
+
1
)
×
i, respectively.The Lanczos procedure guarantees that the vectorsv
j
+
1
⊥
b
V
~
j and v
~
j
+
1
⊥
b
V
j. Since W
~
T
A
1
=
0 and W
T
A
~
1
=
0, we have v
j
+
1
T
W
~
=
0 and v
~
j
+
1
T
W
=
0 for j
=
1,2
,
…, via (5) and (6). Hence the sequences {
v
j
}
j
=
1,2
,
… and {
v
~
j
}
j
=
1,2
,
… satisfy the properties of (2).To make these ideas more concrete, we give a deflated Lanczos algorithm by substituting the right-hand sides of (3) and (4).


## 3. The Deflated BiCG Algorithm

In this section, we derive the deflated BiCG method from Algorithm1 in exactly the same way as BiCG was derived from the Lanczos biorthogonalization procedure [17]. For convenience, we drop the superscript i in (1) and refer to(7)
A
x
=
b
,
A
T
x
~
=
b
~as the primary and the dual system, respectively.Algorithm 1:Deflated Lanczos algorithm.
DefineW
=
[
w
1
,
w
2
,
…
,
w
k
] and W
~
=
[
w
~
1
,
w
~
2
,
…
,
w
~
k
].
Choose initial unit vectorsv
1 and v
~
1 such that v
1
T
W
~
=
0 and v
~
1
T
W
=
0. Set
v
0
=
0
,
v
~
0
=
0.
for   l
=
1,2
,
…   do
SolveW
~
T
A
W
μ
j
=
W
~
T
A
v
j; W
T
A
T
W
~
μ
~
j
=
W
T
A
T
v
~
j for μ
j
,
μ
~
j
z
j
=
A
v
j
-
A
W
μ
j;    z
~
j
=
A
T
v
~
j
-
A
T
W
~
μ
~
j
ρ
j
=
(
v
~
j
,
z
j
);    ρ
~
j
=
(
v
j
,
z
~
j
)
v
^
j
+
1
=
z
j
-
σ
j
v
j
-
1
-
ρ
j
v
j;    v
~
^
j
+
1
=
z
~
j
-
σ
~
j
v
~
j
-
1
-
ρ
~
j
v
~
j
σ
j
+
1
=
(
v
~
^
j
+
1
,
v
^
j
+
1
)
1
/
2;    σ
~
j
+
1
=
v
~
^
j
+
1
,
v
^
j
+
1
σ
j
+
1
v
j
+
1
=
v
^
j
+
1
σ
j
+
1;   v
~
j
+
1
=
v
~
^
j
+
1
σ
~
j
+
1
end forAssume initial guessesx
0 and x
~
0 are given such that r
0
=
b
-
A
x
0
⊥
b
W
~ and r
~
0
=
b
~
-
A
T
x
~
0
⊥
b
W. Let v
1
=
r
0
/
σ
1, v
~
1
=
r
~
0
/
σ
~
1. Define(8)
K
k
,
j
A
,
W
,
r
0
≡
s
p
a
n
W
,
V
j
,
K
k
,
j
A
~
,
W
~
,
r
~
0
≡
s
p
a
n
W
~
,
V
~
j
.The jth solutions that are updated in deflated BiCG algorithm become (9)
x
j
=
x
0
+
W
ζ
j
+
V
j
η
j
∈
K
k
,
j
A
,
W
,
r
0
,
(10)
x
~
j
=
x
~
0
+
W
~
ζ
j
+
V
~
j
η
j
∈
K
k
,
j
A
~
,
W
~
,
r
~
0
.With deflation, the biorthogonality conditions V
j
⊥
b
W
~, V
~
j
⊥
b
W, and V
~
j
⊥
b
V
j define the Petrov Galerkin conditions (11)
r
j
=
r
0
-
A
W
ζ
j
-
A
V
j
η
j
⊥
b
W
~
,
V
~
j
,
(12)
r
~
j
=
r
~
0
-
A
T
W
~
ζ
j
-
A
T
V
~
j
η
j
⊥
b
W
,
V
j
.Lemma 1.
Ifx
j, x
~
j, r
j, and r
~
j satisfy (9), (10), (11), and (12), respectively, then r
j
=
c
j
v
j
+
1 and r
~
j
=
c
~
j
v
~
j
+
1 with c
j, c
~
j
∈
R. Thus the residual r
j is orthogonal to K
k
,
j
(
A
~
,
W
~
,
r
~
0
)
≡
s
p
a
n
{
W
~
,
r
~
0
,
r
~
1
,
…
,
r
~
j
-
1
} and r
~
j is orthogonal to K
k
,
j
(
A
,
W
,
r
0
)
≡
s
p
a
n
{
W
,
r
0
,
r
1
,
…
,
r
j
}.Proof.
From (3) and (5), we have (13)
A
V
j
=
A
W
Υ
j
+
V
j
T
j
+
t
j
+
1
,
j
v
j
+
1
e
j
T
,where Υ
j
∶
=
W
~
T
A
W
-
1
W
~
T
A
V
j. Using (9), we get (14)
r
j
=
r
0
-
A
W
ζ
j
-
A
V
j
η
j
,
=
r
0
-
A
W
ζ
j
-
A
W
Υ
j
+
V
j
T
j
+
t
j
+
1
,
j
v
j
+
1
e
j
T
η
j
.
From the orthogonality conditionsr
j
⊥
b
[
W
~
,
V
~
j
] leads to the following equations: (15)
ζ
j
=
-
Υ
j
η
j
,
η
j
=
σ
1
T
j
-
1
e
1
.Substituting (15) into (14), we get r
j
=
c
j
v
j
+
1 for some scalar c
j. Similarly, it is able to prove that r
~
j
=
c
~
j
v
~
j
+
1 and r
~
j
⊥
b
[
W
,
V
j
].Now we defineT
j
=
L
j
D
j
U
j, where D
j, L
j, and U
j are diagonal, lower triangular, and upper triangular, respectively. Since a pivotless LDU decomposition of the tridiagonal matrix T
j will lead to a breakdown in BiCG, it can be avoided by performing the LDU decomposition with 2 × 2 block diagonal elements [18]. However, such breakdown rarely happens in practice; hence we will not discuss it for the sake of brevity. Define(16)
P
j
≡
p
0
,
p
1
,
…
,
p
j
-
1
=
-
W
Υ
j
+
V
j
U
j
-
1
Λ
j
,
P
~
j
≡
p
~
0
,
p
~
1
,
…
,
p
~
j
-
1
=
-
W
~
Υ
~
j
+
V
~
j
L
j
-
T
Λ
~
j
,where (17)
Υ
~
j
=
W
T
A
T
W
~
-
1
W
T
A
T
V
~
j
,
Λ
j
=
d
i
a
g
c
0
,
c
1
,
…
,
c
j
-
1
,
Λ
~
j
=
d
i
a
g
c
~
0
,
c
~
1
,
…
,
c
~
j
-
1
.Using the above information, two-term recurrences for dual linear systems can be easily derived from the deflated Lanczos procedure.Theorem 22.
The solutionsx
j, x
~
j, the residuals r
j, r
~
j, and the descent directions p
j, p
~
j for the dual linear systems satisfy the following recurrence relations:(18)
x
j
=
x
j
-
1
+
α
j
-
1
p
j
-
1
;
x
~
j
=
x
~
j
-
1
+
α
j
-
1
p
~
j
-
1
,
r
j
=
r
j
-
1
+
α
j
-
1
p
j
-
1
;
r
~
j
=
r
~
j
-
1
+
α
j
-
1
A
T
p
~
j
-
1
,
p
j
=
β
j
-
1
p
j
-
1
+
r
j
-
W
μ
j
;
p
~
j
=
β
j
-
1
p
~
j
-
1
+
r
~
j
-
W
~
μ
~
j
.Proof.
It is analogous to the proof of Theorem2.1 [7].Theorem 33.
The descent directionp
~
j is A-orthogonal to p
i
i
=
1,2
,
…
j. In addition, it is also A-orthogonal to w
i
i
=
1,2
,
…
j.Proof.
It is equivalent to prove thatP
~
j
T
A
P
j is a diagonal matrix and W
~
j
T
A
P
j
=
0. Consider the following:(19)
P
~
j
T
A
P
j
=
Λ
~
j
T
L
j
-
1
-
W
~
Υ
~
j
+
V
~
j
T
A
-
W
Υ
j
+
V
j
U
j
-
1
Λ
j
=
Λ
~
j
T
L
j
-
1
V
~
j
T
A
V
j
-
V
~
j
T
A
W
Υ
j
U
j
-
1
Λ
j
=
Λ
~
j
T
L
j
-
1
V
~
j
T
I
-
A
W
W
~
T
A
W
-
1
W
~
T
A
V
j
U
j
-
1
Λ
j
=
Λ
~
j
T
L
j
-
1
V
~
j
T
A
1
V
j
U
j
-
1
Λ
j
=
Λ
~
j
T
D
J
Λ
jis diagonal and (20)
W
~
j
T
A
P
j
=
-
W
~
T
A
W
W
~
T
A
W
-
1
W
~
T
A
V
j
+
W
~
T
A
V
j
-
W
~
T
A
W
W
~
T
A
W
-
1
W
~
T
A
V
j
U
j
-
1
Λ
j
=
0
.The proof of W
j
T
A
T
P
~
j
=
0 goes along the same lines as above. We next find expressions for α
j
′
s, α
~
j
′
s, β
j
′
s, and β
~
j
′
s with iteration vectors at hand.Theorem 4.
The coefficients in the deflated BiCG method satisfy the following relations:(21)
α
j
-
1
=
r
~
j
-
1
,
r
j
-
1
p
~
j
-
1
,
A
p
j
-
1
;
μ
j
=
W
T
A
W
-
1
W
T
A
r
j
;
β
j
-
1
=
r
~
j
,
r
j
r
~
j
-
1
,
r
j
-
1
;
μ
~
j
=
W
T
A
T
W
~
-
1
W
T
A
T
r
~
j
,where α
~
j
-
1
=
α
-
j
-
1 and β
~
j
-
1
=
β
-
j
-
1.Proof.
The proof follows from [7].Algorithm2 provides an outline of the deflated BiCG method.Algorithm 2:Deflated BiCG algorithm.
(1) GivenW and W
~. If W and W
~ are not available; then initialize W and W
~ to empty matrices.
(2) Choosex
-
1, x
~
-
1 and compute x
0, x
~
0, r
0, r
~
0 using (22).
(3) if(
r
0
,
r
~
0
)
=
0 then x
~
-
1 to a random vector.
(4) SolveW
~
T
A
W
μ
0
=
W
~
T
A
r
0; W
T
A
T
W
~
μ
~
0
=
W
T
A
T
r
~
0
(5) Setp
0
=
r
0
-
W
μ
0; p
~
0
=
r
~
0
-
W
~
μ
0.
(6)for   l
=
1,2
,
…   do
(7)α
j
-
1
=
r
~
j
-
1
,
r
j
-
1
p
~
j
-
1
,
A
p
j
-
1
(8)x
j
=
x
j
-
1
+
α
j
-
1
p
j
-
1;    x
~
j
=
x
~
j
-
1
+
α
-
j
-
1
p
~
j
-
1
(9)r
j
=
r
j
-
1
+
α
j
-
1
p
j
-
1;    r
~
j
=
r
~
j
-
1
+
α
-
j
-
1
A
T
p
~
j
-
1
(10)if   |
|
r
j
|
|
≤ tol and |
|
r
~
j
|
|
≤ tol  then
(11)   break;
(12)end if
(13)β
j
-
1
=
r
~
j
,
r
j
r
~
j
-
1
,
r
j
-
1
(14)  SolveW
~
T
A
W
μ
j
=
W
~
T
A
r
j; W
T
A
T
W
~
μ
~
j
=
W
T
A
T
r
~
j
(15)p
j
=
β
j
-
1
p
j
-
1
+
r
j
-
W
μ
j;    p
~
j
=
β
-
j
-
1
p
~
j
-
1
+
r
~
j
-
W
~
μ
~
j
(16)end forIn order to guaranteer
0
⊥
b
W
~ and r
~
0
⊥
b
W, we take the following forms. Take x
-
1 and x
-
1 as the initial guesses and r
-
1
=
b
-
A
x
-
1, r
~
-
1
=
b
-
A
T
x
~
-
1 as the corresponding initial residuals. Let (22)
x
0
=
W
W
~
T
A
W
-
1
W
~
T
r
-
1
;
r
0
=
r
-
1
-
A
W
W
~
T
A
W
-
1
W
~
T
r
-
1
,
x
~
0
=
W
~
W
T
A
T
W
~
-
1
W
T
r
~
-
1
;
r
~
0
=
r
~
-
1
-
A
T
W
~
W
T
A
T
W
~
-
1
W
T
r
~
-
1
.


## 4. Computing Approximate Invariant Subspaces

The matricesW and W
~ are defined as the primary and dual system deflated spaces, respectively. The deflated spaces are fixed to solve the current linear systems; however, the bases of the deflated spaces are updated periodically for the next pair of linear systems. In the literature, approximate invariant subspace of A is taken as deflated space. Since we focus on the solutions of the dual linear systems, it is necessary to build the left invariant subspace for the dual and the right invariant for the primary. We use harmonic Ritz vectors [19] to approximate left and right invariant subspaces cheaply. Hence it is needed to solve a small generalized eigenvalue problem whose solution gives the desired approximate invariant subspaces. Although it is cheap to solve the generalized eigenvalue problem, it would be expensive to set up the corresponding matrices in a straightforward manner. For the less matrix vector products, we use the descent directions p
j, p
~
j to update the deflated subspaces following [7].LetW
(
s
)
=
[
w
1
(
s
)
,
w
2
(
s
)
,
…
,
w
k
(
s
)
], W
~
(
s
)
=
[
w
~
1
(
s
)
,
w
~
2
(
s
)
,
…
,
w
~
k
(
s
)
]
. For each new system, we define the matrices Z,  Z
~ as follows: (23)
Z
s
=
W
s
,
P
l
s
,
Z
~
s
=
W
~
s
,
P
~
l
s
,where (24)
P
l
s
=
p
1
s
,
p
2
s
,
…
,
p
l
s
,
P
~
l
s
=
p
~
1
s
,
p
~
2
s
,
…
,
p
~
l
s
.The harmonic projection method computes the k left and right approximate eigenvectors corresponding to the k smallest eigenvalues by solving the following generalized eigenproblems: (25)
G
s
y
i
-
θ
i
F
s
y
i
=
0
,
i
=
1
,
…
,
k
,
G
(
s
)
T
y
~
i
-
θ
~
i
F
(
s
)
T
y
~
i
=
0
,
i
=
1
,
…
,
k
,where G
(
s
)
=
A
~
T
Z
~
s
T
A
Z
(
s
) and F
(
s
)
=
A
T
Z
~
s
T
Z
(
s
). Let Y
(
s
)
=
[
y
1
(
s
)
,
y
2
(
s
)
,
…
,
y
k
(
s
)
], Y
~
(
s
)
=
[
y
~
1
(
s
)
,
y
~
2
(
s
)
,
…
,
y
~
k
(
s
)
]; then the new deflated subspaces for the next system are defined as (26)
W
s
+
1
=
Z
s
Y
s
,
W
~
s
+
1
=
Z
~
s
Y
~
s
.For avoiding the matrix multiplication and reducing the number of the matrix vector products in G
(
s
) and F
(
s
), we take the same way as [7].Lemma 55.
LetR
l
(
s
)
=
[
r
1
(
s
)
,
r
2
(
s
)
,
…
,
r
l
(
s
)
], R
~
l
(
s
)
=
[
r
~
1
(
s
)
,
r
~
2
(
s
)
,
…
,
r
~
l
(
s
)
]. One has Δ
l
+
1
(
s
)
=
μ
0
(
s
)
,
μ
1
(
s
)
,
…
,
μ
l
(
s
), Δ
~
l
+
1
(
s
)
=
μ
~
0
(
s
)
,
μ
~
1
(
s
)
,
…
,
μ
~
l
(
s
), and(27)
L
l
s
=
1
α
0
s
-
1
α
0
s
1
α
1
s
-
1
α
1
s
⋱
⋱
1
α
l
-
1
s
-
1
α
l
-
1
s
,
U
l
s
=
1
-
β
0
s
1
-
β
1
s
⋱
⋱
1
-
β
l
-
1
s
1
,where L
l
(
s
), U
l
(
s
) take α
-
j
(
s
) and -
β
-
,
(
s
) respectively. Then the matrices A
P
l
(
s
), A
T
P
~
l
(
s
) can be computed by (28)
A
P
l
s
=
W
s
Δ
l
+
1
+
P
l
+
1
s
U
l
s
L
l
,
A
T
P
~
l
s
=
W
~
s
Δ
~
l
+
1
+
P
~
l
+
1
s
U
~
l
s
L
~
l
.Proof.
The proof follows from [7].Lemma 66.
Define(29)
D
^
l
s
=
d
i
a
g
d
0
s
,
d
1
s
,
…
,
d
l
-
1
s
,
w
h
e
r
e
d
j
s
=
p
~
j
s
T
A
p
j
s
,
G
^
s
i
,
j
=
d
j
-
1
s
α
j
-
1
s
1
+
β
j
-
1
s
,
i
=
j
,
G
^
s
i
,
j
=
-
d
i
s
α
i
-
1
s
,
i
=
j
+
1
,
G
^
s
i
,
j
=
G
^
s
j
,
i
,
i
=
j
-
1
,
G
^
s
i
,
j
=
0
,
e
l
s
e
.
Then the matricesG
(
s
), F
(
s
), A
W
(
s
), and A
T
W
~
(
s
+
1
) are given by (30)
G
s
=
A
T
W
~
s
T
A
W
s
W
~
s
T
A
W
s
Δ
l
+
1
L
l
Δ
~
l
+
1
L
~
l
T
A
W
s
G
^
s
,
(31)
F
s
=
W
~
s
T
A
W
s
0
0
D
^
s
,
(32)
A
W
s
+
1
=
A
W
s
A
P
l
s
Y
s
,
(33)
A
T
W
~
s
+
1
=
A
T
W
~
s
A
T
P
~
l
s
Y
~
s
.Proof.
The proof is analogous to the proofs as given in [7].


## 5. Systems with Multiple Right-Hand Sides

In this section, we describe the deflated BiCG algorithm for solving the dual linear systems with multiple RHSs. The algorithm uses descent direction vectorsP
l, P
~
l to improve approximate eigenvectors found by subsequent runs of deflated BiCG and uses deflation to accelerate convergence. The new vectors P
l, P
~
l are appended to the current deflated subspaces, W
(
s
) and W
~
(
s
). These incrementally built spaces are then used to generate approximate subspaces for deflating the next systems. The resulting algorithm is given in Algorithm 3.Algorithm 3:Deflated BiCG algorithm for multiple right-hand sides.
(1) Choosel and k with l
≥
k.
(2) Solve the first dual linear systems of (1) with the standard BiCG.
(3) ComputeG
(
1
) and F
(
1
) using (30) and (31). Set W
(
1
)
=
[
⌀
] and W
~
(
1
)
=
[
⌀
].
(4)for   s
=
2,3
,
…,   do
(5)  Solve fork left and right eigenvectors. Set W
(
s
+
1
)
=
Z
(
s
)
Y
(
s
), W
~
(
s
+
1
)
=
Z
~
(
s
)
Y
~
(
s
).
(6)  Choosel
≥
0.
(7)  Solve thesth system of (1) by deflated BiCG with W
=
W
(
s
), W
~
=
W
~
(
s
+
1
).
(8)  ComputeG
(
s
) and F
(
s
). Note that G
~
(
s
)
=
G
(
s
)
T and F
~
(
s
)
=
F
(
s
)
T
(9)end forThe algorithm requires the storage ofW, A
W, W
~, A
T
W
~, P, and P
~ for approximating eigenvectors, amounting to 4
k
+
2
l vectors of length n. In addition, we need to save Δ
l
+
1, Δ
~
l
+
1, α
i
(
s
), β
i
(
s
), and d
i
(
s
) of the first l steps of the algorithm.


## 6. Numerical Examples

In this section we give numerical results which indicate the potential effectiveness of our approach. In all of our runs we use a zero initial guess for the first dual system and take (22) as initial guesses for the remaining systems. Also, we keep the data in first l
=
20 steps and k
=
20 approximate eigenvectors associated with smallest eigenvalues. The stopping criterion is taken as b
-
A
x
k
2
/
b
2
≤
1
0
-
6. All the numerical experiments were performed in MATLAB 2011b. The machine we have used is a PC Pentium(R)4, CPU 2.50 GHz, 2.00 GB of RAM.


### 6.1. Convection Diffusion Problem

In this subsection, we take two different approaches BiCG and deflated BiCG (DBiCG) to solve the same problem for highlighting the efficiency of the deflation. The example is the convection diffusion problem with Dirichlet boundary conditions in the unit square0,1
2
=
Ω
∩
∂
Ω, where Ω
=
0,1
2, (34)
-
ε
Δ
u
+
c
u
x
+
d
u
y
=
g
i
n
Ω
,
u
=
0
o
n
∂
Ω
.This problem is discretized with a second-order finite difference scheme for a vertex centered location of unknowns. We adopt h
/
ε
m
a
x
(
c
,
d
)
=
2 to satisfy the Péclet condition, where h
=
1
/
N
-
1 is the mesh size and N is the number of points per direction. Here, h
=
1
/
32, c
=
d
=
128, and ε
=
1 are taken; then the problem size is 1024 and has 4992 nonzero entries, which is real, sparse, and nonsymmetric. The primary system right-hand side comes from PDE. For convenience, we take the same vector as the dual system right-hand side. Their convergence behaviors are plotted in Figure 1.Figure 1
(a) Primary system and (b) dual system.
(a)
(b)Figure1 shows the convergence histories of two methods. Using the initial guesses (22), DBiCG (without any deflated subspaces) converges rather slowly in the first run. The numerical behavior of DBiCG will be the same as the ones by applying BiCG to solve the first dual linear systems. The only difference is that the deflated spaces for the next run are generated during the first process of DBiCG. For the second run, we solve the same systems by DBiCG-sh with deflated spaces. Compared with the first run, the reduction in iteration is around 25
% in the second run.


### 6.2. Lattice QCD Problem

The second numerical experiments stem from a quantum chromodynamics (QCD) problem. QCD is the gauge theory that describes the strong nuclear force between quarks and gluons [20, 21]. Quark propagators are obtained by solving the inhomogeneous lattice Dirac equation A
x
i
=
b
i, where A is a large but sparse complex non-Hermitian matrix representing a periodic nearest neighbour coupling on a four-dimensional Euclidean space time lattice. Following [22], we take a preconditioning process which transforms from the original system to the odd even reduced system. That preconditioning process has one-half of the dimension and a smaller condition number than the original system; for details see [23, 24]. The number of right-hand sides is 4 and the b
i is taken as a vector of all ones and three independent random vectors with entries form N
(
0,1
).Figure2 is helpful for seeing how the deflation works in practise. It is observed that DBiCG leads to less iterations than BiCG.Figure 2
(a) Primary system and (b) dual system.
(a)
(b)


## 6.1. Convection Diffusion Problem

In this subsection, we take two different approaches BiCG and deflated BiCG (DBiCG) to solve the same problem for highlighting the efficiency of the deflation. The example is the convection diffusion problem with Dirichlet boundary conditions in the unit square0,1
2
=
Ω
∩
∂
Ω, where Ω
=
0,1
2, (34)
-
ε
Δ
u
+
c
u
x
+
d
u
y
=
g
i
n
Ω
,
u
=
0
o
n
∂
Ω
.This problem is discretized with a second-order finite difference scheme for a vertex centered location of unknowns. We adopt h
/
ε
m
a
x
(
c
,
d
)
=
2 to satisfy the Péclet condition, where h
=
1
/
N
-
1 is the mesh size and N is the number of points per direction. Here, h
=
1
/
32, c
=
d
=
128, and ε
=
1 are taken; then the problem size is 1024 and has 4992 nonzero entries, which is real, sparse, and nonsymmetric. The primary system right-hand side comes from PDE. For convenience, we take the same vector as the dual system right-hand side. Their convergence behaviors are plotted in Figure 1.Figure 1
(a) Primary system and (b) dual system.
(a)
(b)Figure1 shows the convergence histories of two methods. Using the initial guesses (22), DBiCG (without any deflated subspaces) converges rather slowly in the first run. The numerical behavior of DBiCG will be the same as the ones by applying BiCG to solve the first dual linear systems. The only difference is that the deflated spaces for the next run are generated during the first process of DBiCG. For the second run, we solve the same systems by DBiCG-sh with deflated spaces. Compared with the first run, the reduction in iteration is around 25
% in the second run.


## 6.2. Lattice QCD Problem

The second numerical experiments stem from a quantum chromodynamics (QCD) problem. QCD is the gauge theory that describes the strong nuclear force between quarks and gluons [20, 21]. Quark propagators are obtained by solving the inhomogeneous lattice Dirac equation A
x
i
=
b
i, where A is a large but sparse complex non-Hermitian matrix representing a periodic nearest neighbour coupling on a four-dimensional Euclidean space time lattice. Following [22], we take a preconditioning process which transforms from the original system to the odd even reduced system. That preconditioning process has one-half of the dimension and a smaller condition number than the original system; for details see [23, 24]. The number of right-hand sides is 4 and the b
i is taken as a vector of all ones and three independent random vectors with entries form N
(
0,1
).Figure2 is helpful for seeing how the deflation works in practise. It is observed that DBiCG leads to less iterations than BiCG.Figure 2
(a) Primary system and (b) dual system.
(a)
(b)


## 7. Conclusions

In this paper, we have derived a deflated BiCG method for dual linear systems with multiple RHSs. The proposed algorithm uses harmonic Ritz vectors to approximate left and right invariant subspaces inexpensively via small descent direction vectors found by subsequent runs of deflated BiCG and then derives the deflated subspaces for the next pair of dual linear systems. This process leads to faster convergence for the next pair of systems. Numerical examples illustrate the effectiveness of the proposed method, which will encourage us to apply the deflation technique to the BiCGSTAB algorithm in the near future.



---

*Source: 372109-2015-09-16.xml*