A LOOK-BACK-TYPE RESTARTFOR THE RESTARTED KRYLOV SUBSPACE METHODSTO SOLVE NON-HERMITIAN LINEAR SYSTEMSAKIRA IMAKURA ∗ , TOMOHIRO SOGABE † , AND SHAO-LIANG ZHANG ‡Abstract. In this paper, we investigate the restarted Krylov subspace methods, as typiﬁed by
the GMRES(m) method and the FOM(m) method, for solving non-Hermitian linear systems. We
have recently focused on the restart of the GMRES(m) method and proposed the extension of the
GMRES(m) method based on the error equations. The main purpose of this paper is to apply theextension to other restarted Krylov subspace methods, and propose a speciﬁc restart technique forthe restarted Krylov subspace method. The speciﬁc restart technique is named as a Look-Backtype restart, and is based on an implicit residual polynomial reconstruction via the initial guess. The
comparison analysis based on the residual polynomials and some numerical experiments indicate that
the Look-Back-type restart achieves more eﬃcient convergence results than the traditional restarted
Krylov subspace methods.
Key words. non-Hermitian linear systems, restarted Krylov subspace methods, residual polynomials, a Look-Back-type restart
AMS subject classifications. 65F10

1. Introduction. In this paper, we consider solving large and sparse linear systems of the form:(1.1)

Ax = b,A ∈ Cn×n ,x, b ∈ Cn ,where the coeﬃcient matrix A is assumed to be non-Hermitian and nonsingular. Theselinear systems often arise from the discretization of partial diﬀerential equations inthe ﬁelds of the computational science and engineering.
Direct methods such as the LU decomposition are generally selected in terms ofaccuracy and stability for the small and dense coeﬃcient matrix A. On the other hand,for the case where the coeﬃcient matrix A is large and sparse, iterative methods arewidely used because of their computational costs and storage requirements. Recently,
the Krylov subspace methods are recognized as standard algorithms for large, sparseand non-Hermitian linear systems (1.1); for details see [16, 21].
In this paper, we investigate the restarted Krylov subspace methods, as typiﬁed by
the GMRES(m) method [17] and the FOM(m) method [15], for solving non-Hermitian
linear systems (1.1). The Arnoldi-based Krylov subspace methods have some diﬃculties in terms of the computational costs and storage requirements due to the long-termrecurrence. The (two-sided) Lanczos-based Krylov subspace methods also have somediﬃculties in terms of the pseudo convergence due to accumulations of the round-oﬀ
errors. In order to remedy these diﬃculties, the restart is often applied to the Krylov
subspace methods.
(1)
Let m be the restart frequency, and x0 be the initial guess of the 1st restart cycle.
Then the ℓth restart cycle of the restarted Krylov subspace methods are operated as∗ Graduate School of Systems and Information Engineering, University of Tsukuba, 1-1-1, Tennodai Tsukuba-city, Ibaraki 305-8573, Japan (imakura@cs.tsukuba.ac.jp).
† Graduate School of Information Science and Technology,
Aichi Prefectural University, 1522-3 Ibaragabasama, Kumabari, Nagakute-cho, Aichi-gun, Aichi, 480-1198, Japan(sogabe@ist.aichi-pu.ac.jp).
‡ Graduate School of Engineering, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603,
Japan (zhang@na.cse.nagoya-u.ac.jp).

1

2

A. IMAKURA, T. SOGABE AND S.-L. ZHANG

follows:
1. Solve (approximately) Ax = b by m iterations of some Krylov subspace(ℓ)
method with the initial guess x0 , and get the approximate solution x(ℓ) .
(ℓ+1)
2. Update the initial guess of the (ℓ + 1)th restart cycle: x0:= x(ℓ) , and go
to the (ℓ + 1)th restart cycle.
The restart remedies the diﬃculties of the Krylov subspace methods due tothe long-term recurrence and/or accumulations of the round-oﬀ errors; however, therestart generally slows their convergence. Therefore, in order to improve Step 1, several improvement techniques have been proposed, such as the adaptive preconditioning techniques based on the deﬂation [1, 6, 8], the techniques based on the augmented
Krylov subspace [2, 4, 11–13] and the techniques based on adaptively determining therestart frequency m [3, 14, 19, 22].
On the other hand, we focused on the update of the initial guess in each restart
cycle (Step 2), then we recently proposed the extension of the GMRES(m) method [9].
It is also shown that this is a natural extension in terms of the error equations and theiterative reﬁnement scheme, and it has high potential for more eﬃcient convergencethan the traditional GMRES(m) method [9]. In this regard, however, the convergenceof the extension of the GMRES(m) method did not analyzed enough, and speciﬁc
algorithms for eﬃcient convergence still have not proposed.
The main purpose of this paper is to apply the extension to other restarted Krylov
subspace methods, and propose a speciﬁc restart technique for the restarted Krylov
subspace method. The speciﬁc restart technique is named as a Look-Back-type restart,
and is based on an implicit residual polynomial reconstruction via the initial guess.
The performance of the Look-Back-type restart is evaluated by a comparison analysis
based on the residual polynomials and some numerical experiments.
This paper is organized as follows. In Section 2, we brieﬂy describe a generalform of the restarted Krylov subspace methods. In Section 3, from the analysis
based on the residual polynomials, we propose a Look-Back-type restart based on an
implicit residual polynomial reconstruction. The reconstructed residual polynomials
are compared with the traditional restarted Krylov subspace methods in Section 4.
In Section 5, the performance of the Look-Back-type restart is evaluated by somenumerical experiments. Our conclusions are summarized in Section 6.
Throughout this paper, let V , W be the subspaces, then V + W denotes the sumof the subspaces of V , W , i.e., V + W := {v + w|v ∈ V , w ∈ W }.
2. Restarted Krylov subspace methods. We have the following theorem forthe exact solution of the linear systems (1.1); see e.g., [16].
Theorem 2.1. Let x0 be an initial guess and r0 := b − Ax0 be the corresponding
initial residual, respectively. We also let d be the grade of r0 with respect to A, i.e.,d := d(A, r0 ) := min{k|Pk (A)r0 = 0, Pk (λ) ∈ Pk , Pk (0) = 1},
where Pk is the set of k-degree polynomials. Then the exact solution x∗ := A−1 bsatisfies the following relation:x∗ ∈ x0 + Kd (A, r0 ),Kd (A, r0 ) := span{r0 , Ar0 , . . . , Ad−1 r0 }.

For a general property of the Krylov subspaces, we can also derive the following
lemma.

LOOK-BACK-TYPE RESTART FOR KRYLOV SUBSPACE METHODS

3

Lemma 2.2. Let A ∈ Cn×n , u, w ∈ Cn . Then, for any v ∈ Ks+1 (A, u),
w + Kt (A, v − Aw) ⊂ Ks+t (A, u)(2.1)
is satisfied if and only if

w ∈ Ks (A, u).

(2.2)

We assume, however, that s + t < d, where d := d(A, u) is the grade of u with respect
to A.
Proof. First, we prove that (2.1) is a necessary condition for (2.2). If w ∈Ks (A, u), then v − Aw ∈ Ks+1 (A, u) for any v ∈ Ks+1 (A, u). This leads to
Kt (A, v − Aw) ⊂ Ks+t (A, u).
Next, we prove that (2.1) is a suﬃcient condition for (2.2). From (2.1), we have
(2.3)

w ∈ Ks+t (A, u),(2.4)

Kt (A, v − Aw) ⊂ Ks+t (A, u).

Then, we introduce a condition for the vector w from (2.4), then consider the intersection with (2.3).
Since a general property of the Krylov subspaces Kk (A, x), i.e.,
Kk (A, x + y) ⊂ Kk (A, x) + Kk (A, y),we have(2.5)

Kt (A, Aw) ⊂ Kt (A, v) + Kt (A, v − Aw).

On the other hand, for any v ∈ Ks+1 (A, u), we also have
Kt (A, v) ⊂ Ks+t (A, u).

(2.6)

Then, substituting (2.4) and (2.6) in (2.5), the relation
Kt (A, Aw) ⊂ Ks+t (A, u)
is derived. This leads to the following condition for the vector w:
At w ∈ Ks+t (A, u).

(2.7)

We then consider the intersection of (2.3) and (2.7). From (2.3), the vector w
can be written by the following polynomial form:
( s−1
)
s+t−1
s+t−1
∑
∑
∑
i
i
i
w=
ρi A u =ρi A u +ρi A u .
i=0

i=0

i=s

Substituting the polynomial form in (2.7), we getρs = ρs+1 = · · · = ρs+t−1 = 0.
This means the vector w satisﬁes
w=

s−1
∑

ρi Ai u ∈ Ks (A, u).

i=0Therefore, it is true that (2.1) is a suﬃcient condition for (2.2). Thus Lemma 2.2 isproved.
In what follows, based on Theorem 2.1 and Lemma 2.2, we consider a relationshipamong the sequence of the iterative solutions constructed from the Krylov subspaces
and describe a general form of the restarted Krylov subspace methods.

4

A. IMAKURA, T. SOGABE AND S.-L. ZHANG

2.1. A relationship among the sequence of the iterative solutions con(ℓ)
structed from the Krylov subspaces. For ℓ = 1, 2, . . . , let x0 be a sequence of(ℓ)(ℓ)
initial guesses, and r0 := b − Ax0 be the sequence of the corresponding residual
vectors, respectively. Applying the sequence of m(ℓ) dimensional Krylov subspaces(ℓ)Km(ℓ) (A, r0 ) for iterative solutions of the linear systems (1.1), iterative solutions(ℓ)
x can be expressed as
(2.8)

(ℓ)

x(ℓ) = x0 + z (ℓ) ,(ℓ)

z (ℓ) ∈ Km(ℓ) (A, r0 ),ℓ = 1, 2, . . . ,where the vectors z (ℓ) are designed by some condition, e.g., the minimal residual
condition, the Ritz-Galerkin condition or the Petrov-Galerkin condition.
Note that in the case of
(ℓ+1)

(2.9)

x0

= x(ℓ) ,ℓ = 1, 2, . . . ,if one designs the vector z (ℓ) using the minimal residual condition and set m(ℓ) = m,then the sequence of x(ℓ) is just the sequence of the approximate solutions of the
GMRES(m) method. In this case, the vectors x(ℓ) have
(1)

(1)

(2)

(ℓ)

x(ℓ) ∈ x0 + Km (A, r0 ) + Km (A, r0 ) + · · · + Km (A, r0 )
(1)

(1)

= x0 + Km×ℓ (A, r0 ).
This is a fortunate property for solving the linear systems because of Theorem 2.1.
Here, let us consider a more general case of
(ℓ+1)

x0

(2.10)= x(ℓ) + y (ℓ+1) ,ℓ = 1, 2, . . . .

Then in this section, we discuss a condition for the vector y (ℓ+1) in which all thevectors x(ℓ) on the aﬃne space:
(1)

(1)

x(ℓ) ∈ x0 + KM (ℓ) (A, r0 ),M (ℓ) =ℓ∑

m(i) ,ℓ = 1, 2, . . . ,i=1

as well as the case of (2.9).
(ℓ+1)
(ℓ)
),
Lemma 2.3. For any z (ℓ) ∈ Km(ℓ) (A, r0 ) and z (ℓ+1) ∈ Km(ℓ+1) (A, r0
(ℓ)

(ℓ)

x(ℓ+1) ∈ x0 + Km(ℓ) +m(ℓ+1) (A, r0 )is satisfied if and only if(ℓ)

y (ℓ+1) ∈ Km(ℓ) (A, r0 ).
(ℓ)

We assume, however, that m(ℓ) + m(ℓ+1) < d, where d := d(A, r0 ) is the grade of(ℓ)
r0 with respect to A.
(ℓ+1)
Proof. From the deﬁnition of the vector x(ℓ+1) and r0= r (ℓ) − Ay (ℓ+1) ,
(ℓ)

(ℓ)

x(ℓ+1) ∈ x0 + Km(ℓ) (A, r0 ) + y (ℓ+1) + Km(ℓ+1) (A, r (ℓ) − Ay (ℓ+1) ),(ℓ)

for any z (ℓ) ∈ Km(ℓ) (A, r0 ). Here, from Lemma 2.2,
(ℓ)

(ℓ)

y (ℓ+1) +Km(ℓ+1) (A, r (ℓ) −Ay (ℓ+1) ) ⊂ Km(ℓ) +m(ℓ+1) (A, r0 ) ⇔ y (ℓ+1) ∈ Km(ℓ) (A, r0 ),LOOK-BACK-TYPE RESTART FOR KRYLOV SUBSPACE METHODS

5

Algorithm 1 A general form of the restarted Krylov subspace method(1)

1: Choose the initial guess x0
2: For ℓ = 1, 2, . . . , until convergence Do:

Set the restart frequency m(ℓ) and the KS(ℓ) method
Solve (approximately) Ax = b by m(ℓ) iterations of the KS(ℓ) method with the(ℓ)
initial guess x0 , and get the approximate solution x(ℓ)
∑ℓ
(1)
5:
Set the vector y (ℓ+1) ∈ KM (ℓ) (A, r0 ), where M (ℓ) = i=1 m(i)(ℓ+1)
6:Update the initial guess x0:= x(ℓ) + y (ℓ+1)
7: End For
3:
4:

(ℓ)

(ℓ)

because of r (ℓ) ∈ Km(ℓ) +1 (A, r0 ) for any z (ℓ) ∈ Km(ℓ) (A, r0 ). Therefore, also using
the following obvious relation:
(ℓ)

(ℓ)

Km(ℓ) (A, r0 ) ⊂ Km(ℓ) +m(ℓ+1) (A, r0 ),Lemma 2.3 is proved.
(ℓ)
Proposition 2.4. For any z (ℓ) ∈ Km(ℓ) (A, r0 ), ℓ = 1, 2, . . .,
(1)

(1)

x(ℓ) ∈ x0 + KM (ℓ) (A, r0 ),ℓ = 1, 2, . . . ,is satisfied if and only if(1)y (ℓ+1) ∈ KM (ℓ) (A, r0 ),ℓ = 1, 2, . . . ,∑ℓ(1)where M (ℓ) = i=1 m(i) . We assume, however, that M (ℓ) < d, where d := d(A, r0 )(1)
is the grade of r0 with respect to A.
Proof. Proposition 2.4 is proved by using Lemma 2.3 recursively, starting fromℓ = 1.
2.2. A general form of the restarted Krylov subspace method. We consider the sequential process which reﬁnes the approximate solution of the linear system(1.1) by Eqs. (2.8) and (2.10). This sequential process can be regarded as a generalform of the restarted Krylov subspace method; see Algorithm 1, which can dynamically set the restart frequency m(ℓ) and the using Krylov subspace method KS(ℓ) in
each restart cycle, and which also update the initial guess by Eq. (2.10). Note thatthe algorithms of the extension of the GMRES(m) method proposed in [9] can be
introduced from the Algorithm 1 [9, Algorithm 3.2].
In this regard, however it is expected that the approximate solution x(ℓ) obtainedby Algorithm 1 satisﬁes
(1)

(1)

x(ℓ) ∈ x0 + KM (ℓ) (A, r0 ),because of Theorem 2.1. Therefore, in this paper, we impose(1)y (ℓ+1) ∈ KM (ℓ) (A, r0 )on the vector y (ℓ+1) in each restart cycle of Algorithm 1; see Proposition 2.4.

6

A. IMAKURA, T. SOGABE AND S.-L. ZHANG

3. A Look-Back-type restart. In this section, from the analysis based on
the residual polynomials, we propose a Look-Back-type restart based on an implicit
residual polynomial reconstruction via the initial guess.
In Section 3.1, we analyze residual polynomials of the general form of the restarted
Krylov subspace method (Algorithm 1) and describe an implicit residual polynomialreconstruction via the initial guess. In Section 3.2, we propose a Look-Back-typerestart based on an implicit residual polynomial reconstruction, and its implementation details are shown in Section 3.3.
3.1. An implicit residual polynomial reconstruction via the initial guess.
Recall that for the general form of the restarted Krylov subspace method (Algorithm 1), the initial guess and the corresponding residual are updated by(ℓ+1)(3.1)x0= x(ℓ) + y (ℓ+1) ,

(ℓ+1)

r0

= r (ℓ) − Ay (ℓ+1) ,(1)y (ℓ+1) ∈ KM (ℓ) (A, r0 ),in each restart cycle. Here, let us focus on the residual vector r (ℓ+1) of the (ℓ + 1)th
restart cycle, and consider the corresponding residual polynomial.
(1)
Because of Proposition 2.4, we have r (ℓ) ∈ KM (ℓ) +1 (A, r0 ). Thus, there exists a(ℓ)
polynomial QM (ℓ) (λ) ∈ PM (ℓ) such that
(ℓ)

(1)

r (ℓ) = QM (ℓ) (A)r0 ,(ℓ+1)

QM (ℓ) (0) = 1.
(ℓ+1)

From this observation, let Pm(ℓ+1) (λ) ∈ Pm(ℓ+1) , Pm(ℓ+1) (0) = 1 be the residual polynomial corresponding to m(ℓ+1) iterations of the KS(ℓ+1) method of the (ℓ + 1)th restart(ℓ)
cycle, then the residual vector r (ℓ+1) can be shown by using QM (ℓ) (λ) like
(ℓ+1)

(ℓ+1)

r (ℓ+1) = Pm(ℓ+1) (A)r0(ℓ+1)= Pm(ℓ+1) (A)(r (ℓ) − Ay (ℓ+1) )
(ℓ+1)

(ℓ)

(1)

= Pm(ℓ+1) (A)(QM (ℓ) (A)r0 − Ay (ℓ+1) ).
(1)

For any vectors y (ℓ+1) ∈ KM (ℓ) (A, r0 ), from Proposition 2.4, there is also a polynoe (ℓ)(ℓ) (λ) ∈ PM (ℓ) , Q
e (ℓ)(ℓ) (0) = 1 such thatmial QM
M(3.2)

e (ℓ)(ℓ) (A)r (1) = Q(ℓ)(ℓ) (A)r (1) − Ay (ℓ+1) .
Q
0
0
M
M

Therefore, the residual vector r (ℓ+1) can be rewritten as follows:e (ℓ) (A)r .
r (ℓ+1) = QM (ℓ+1) (A)r0 = Pm(ℓ+1) (A)Q
0
M
(ℓ+1)

(1)

(ℓ+1)

(ℓ)

(1)

As a consequence, we can observe that the vector y (ℓ+1) of the general form of therestarted Krylov subspace method reconstructs the residual polynomial by Eq. (3.2).
In this paper we call this an implicit residual polynomial reconstruction via the initial
guess.
3.2. Proposal for a Look-Back-type restart. Determining the reconstructed
e (ℓ)(ℓ) (λ) particularly, the corresponding vector y (ℓ+1) can be computed
polynomial Q
M
uniquely as follows:()(ℓ)
e (ℓ)(ℓ) (A) r (1) ,(3.3)
y (ℓ+1) = A−1 QM (ℓ) (A) − Q
0MLOOK-BACK-TYPE RESTART FOR KRYLOV SUBSPACE METHODS

7

from Eq. (3.2).
e (ℓ)(ℓ) (λ) should be determined
In terms of computational cost, the polynomial Q
Msuch that Eq. (3.3) is easy to compute. From this observation, in this paper, wee (ℓ)(ℓ) (λ) as the linear combination
express the reconstructed residual polynomial QM(ℓ)

of the current residual polynomial QM (ℓ) (λ) and at most M (ℓ) -degree polynomial(ℓ)(ℓ)
RM (ℓ) (λ) ∈ PM (ℓ) , RM (ℓ) (0) = 1, i.e.,e (ℓ)(ℓ) (λ) := τ (ℓ) Q(ℓ)(ℓ) (λ) + (1 − τ (ℓ) )R(ℓ)(ℓ) (λ),
Q
M
M
M
(ℓ)

where τ (ℓ) ∈ C. Then, we deﬁne the polynomial RM (ℓ) (λ) by using the dth previous(ℓ)
polynomial before QM (ℓ) (λ) as follows:


(ℓd )

 QM (ℓd ) (λ)(ℓ)(3.4) RM (ℓ) (λ) := Qe (ℓd(ℓ) ) (λ)
Md(d : even)
,

ℓd :=

(d : odd) ℓ − d2

(d : even)(d : odd)ℓ − d+1
2

with the (ﬁxed) positive integer parameter d ∈ N. Here we note that one can naturally
extend it by adaptively determining the parameter d in each restart cycle.
In the case of the deﬁnition (3.4), Eq. (3.3) can be rewritten by
(
)
(ℓ)
(ℓ)
(1)
y (ℓ+1) = (1 − τ (ℓ) )A−1 QM (ℓ) (A) − RM (ℓ) (A) r0 ,(3.5)
and we have
(

)

(ℓ)
(ℓ)
QM (ℓ) (A) − RM (ℓ) (A)Let ∆x(ℓ) and ∆r (ℓ) be{ (ℓ)
x − x(ℓd )∆x(ℓ) :=(ℓ )x(ℓ) − x0 d

{
(1)
r0 =

(d : even),(d : odd)

r (ℓ) − r (ℓd )(ℓ )r (ℓ) − r0 d{

∆r (ℓ) :=(d : even).
(d : odd)

r (ℓ) − r (ℓd )(ℓ )r (ℓ) − r0 d

(d : even),(d : odd)

respectively, then from ∆r (ℓ) = −A∆x(ℓ) , Eq. (3.5) can be computed byy (ℓ+1) = µ(ℓ) ∆x(ℓ) ,(3.6)

where µ(ℓ) := τ (ℓ) − 1 ∈ C. In this paper, we name the restart technique based on
Eqs. (3.1) and (3.6) as a Look-Back-type restart.
Based on the Look-Back-type restart, we can eﬃciently achieve the implicit residual polynomial reconstruction without explicit computations for the polynomials such(ℓ)(ℓ)
as QM (ℓ) (λ) and RM (ℓ) (λ).
3.3. Implementation details. For any µ(ℓ) ∈ C the Look-Back-type restart
does not satisfy the monotonic decrease of the residual 2-norm ∥r (ℓ) ∥2 even if it isapplied to the GMRES(m) method. We now have the following proposition for themonotonic decrease of the residual 2-norm of the Look-Back GMRES(m) method.
Proposition 3.1. The Look-Back GMRES(m) method satisfies the monotonic
decrease of the residual 2-norm:
(3.7)

(ℓ+1)

∥r (ℓ+1) ∥2 ≤ ∥r0∥2 = ∥r (ℓ) − Ay (ℓ+1) ∥2 ≤ ∥r (ℓ) ∥2 ,ℓ = 1, 2, . . . ,8

A. IMAKURA, T. SOGABE AND S.-L. ZHANG

Algorithm 2 A restarted Krylov subspace method with a Look-Back-type restart(1)

1: Choose the parameter d ≥ 2 and the initial guess x0
2: For ℓ = 1, 2, . . . , until convergence Do:

Set the restart frequency m(ℓ) and the KS(ℓ) method
Solve (approximately) Ax = b by m(ℓ) iterations of the KS(ℓ) method with the(ℓ)
initial guess x0 , and get the approximate solution x(ℓ)
5:Compute the vector y (ℓ+1) as follows:
If ℓ = 1 then y (ℓ+1) = 0
If ℓ ≥ 2 thenIf (ℓ = d = 2) or (d : even, ℓ ≤ d2 ) or (d : odd, ℓ ≤ d−12 ) then
(1)
(ℓ)
(ℓ)
∆x := x − x0
Else
{ (ℓ)
{
x − x(ℓd ) (d : even)(d : even)ℓ − d2
(ℓ)
∆x :=
,
ℓ
:=
d
(ℓ
)
d−1
d
(ℓ)
ℓ− 2
(d : odd)x − x0(d : odd)
End Ify (ℓ+1) = µ(ℓ) ∆x(ℓ) , µ(ℓ) = arg minµ∈C ∥r (ℓ) − µA∆x(ℓ) ∥2
End If(ℓ+1)
6:
Update the initial guess x0:= x(ℓ) + y (ℓ+1)
7: End For
3:4:

as well as the GMRES(m) method if and only if µ(ℓ) is set by
(ℓ)

(ℓ)

|µ(ℓ) − µmin | ≤ |µmin |,(ℓ)µmin = arg min ∥r (ℓ) − µA∆x(ℓ) ∥2 ,µ∈C

ℓ = 1, 2, . . . .

Proof. The ﬁrst inequality of (3.7) is satisﬁed for any µ(ℓ) because of the minimal
residual condition of the GMRES(m) method. Then we prove the necessary andsuﬃcient condition of the second inequality of (3.7).
Since the deﬁnition of y (ℓ+1) (3.6), the inequality ∥r (ℓ) − Ay (ℓ+1) ∥2 ≤ ∥r (ℓ) ∥2 canbe rewritten as ∥r (ℓ) ∥2 − ∥r (ℓ) − µ(ℓ) A∆x(ℓ) ∥2 ≥ 0, and this is equivalence toµ(ℓ) (r (ℓ) , A∆x(ℓ) ) + µ(ℓ) (A∆x(ℓ) , r (ℓ) ) − µ(ℓ) µ(ℓ) (A∆x(ℓ) , A∆x(ℓ) ) ≥ 0.
(ℓ)

Then, using µmin = (A∆x(ℓ) , r (ℓ) )/(A∆x(ℓ) , A∆x(ℓ) ), we have
(ℓ)

(ℓ)

|µmin | − |µ(ℓ) − µmin | ≥ 0.
Therefore Proposition 3.1 is proved.
The algorithm of the restarted Krylov subspace method with the Look-Back-type(ℓ)
restart is shown in Algorithm 2, where µ(ℓ) is set at µ(ℓ) = µmin . The incremental
operations per restart cycle and storage requirements for the Look-Back-type restartcompared with the traditional restart (2.9) are also shown in Table 3.1, where MatVec and AXPY denote the incremental number of matrix-vector multiplications andadditions of scaled vectors respectively, and Storage means main part of incremental
storage requirements.
We note that these incremental operations and storage requirements correspond to
our current implementation. One can decrease the incremental Mat-Vec of the LookBack-type restart; however, such implementation may increase the eﬀect of round-oﬀ
error.

LOOK-BACK-TYPE RESTART FOR KRYLOV SUBSPACE METHODS

9

Table 3.1: The incremental operations per restart cycle and storage requirements for
the Look-Back-type restart compared with the traditional restart (2.9).
Mat-Vec

AXPY

Inner-Product

1

3

2

Storaged:even d:odd
d
d+1
2n
2 n
(ℓ)

The parameter d depends on the degree of the polynomial RM (ℓ) (λ). It is expectedthat quite large parameter d leads to decrease the eﬀect of the Look-Back-type restart,(ℓ)
because the degree of the corresponding polynomial RM (ℓ) (λ) is very small, and also
leads to increase the storage requirements; see Table 3.1. Thus the quite large d is
not recommended. In this paper, we set d = 3.
4. Comparison analysis based on the residual polynomials. Let the coefﬁcient matrix A be diagonalizable: A = XΛX −1 , where Λ = diag{λ1 , . . . , λn } is the
diagonal matrix of the eigenvalues of A, and r0 and rk be the initial residual and theresidual vector obtained by k iterations of the Krylov subspace method, respectively.
We also let Pk (λ) ∈ Pk be the residual polynomial of the Krylov subspace method,i.e., rk = Pk (A)r0 , Pk (0) = 1. Then we have∥rk ∥2 ≤max

i=1,2,...,n|Pk (λi )|∥X∥2 ∥X −1 ∥2 ∥r0 ∥2 .

Thus, the convergence rate of the residual 2-norm depends strongly on the magnitudeof the residual polynomial corresponding to each eigenvalue λi , and it is expected that
Pm (λ) converge uniformly with respect to each eigenvalue λi for eﬃcient convergence.
In the following this section, we apply the typical restarted Krylov subspacemethods and their Look-Back-type methods: the GMRES(m) method; the FOM(m)method; the Look-Back GMRES(m) method; and the Look-Back FOM(m) method,
to the 50 dimensional symmetric positive deﬁnite linear system of the form:(4.1)

Ax = b,A = diag{0.022 , 0.042 , . . . , 1.002 },b = [1, 1, . . . , 1]T .

Then we analyze comparison between the magnitude of the residual polynomial |Pm×ℓ (λi )|with respect to each eigenvalue of A and the convergence rate of each method. For
solving the model problem (4.1), we set m = 5 as the restart frequency, d = 3 as(1)
the parameter for the Look-Back-type restart and x0 = [0, 0, . . . , 0]T for the initial
guess.
The log-plot of |Pm×ℓ (λi )| for 10th, 20th, . . . , 50th restart cycle of each method
are shown in Fig. 4.1, and the residual 2-norm histories are also shown in Fig. 4.2.
We ﬁrstly consider the residual polynomials and the convergence rate of the
GMRES(m) method and the FOM(m) method. From the white square □ in Fig. 4.1,FG(λ) have at most 5 (mul(λ), Pm×ℓwe can observe that the residual polynomials Pm×ℓ
G
F
tiple) roots even ℓ ≥ 2. Therefore |Pm×ℓ (λi )|, |Pm×ℓ (λi )| corresponding to few eigenG
values near the roots are quite small even at early iterations; however, |Pm×ℓ(λi )|,
F|Pm×ℓ (λi )| corresponding to other eigenvalues show relatively large values even ℓ = 50.
As the results of this non-uniform convergence manner of the residual polynomialsGFPm×ℓ(λi ), Pm×ℓ(λi ) with respect to each eigenvalue, the poor convergence of the
GMRES(m) method and the FOM(m) method are caused; see Fig. 4.2.

A. IMAKURA, T. SOGABE AND S.-L. ZHANG

0

0

-2

-2

log10 | P(λ) |log10 | P(λ) |

10

-4
-6
-8
-10
0

0.2

0.4

0.6

-6
-8
-10

GMRES(m)Look-Back GMRES(m)

-12

-4

FOM(m)
Look-Back FOM(m)

-12
0.8

1

0

0.2

0.4

λ

0

0

-2

-2

-4
-6
-8
-10

-4
-6
-8

-12
0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

λ

0.6

0.8

1

λ

(c) (Look-Back) GMRES(m), ℓ = 20(d) (Look-Back) FOM(m), ℓ = 20

0

0

-2

-2

log10 | P(λ) |log10 | P(λ) |

1

-10

-12

-4
-6
-8
-10

-4
-6
-8
-10

-12

-12
0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

λ

0.6

0.8

1

λ

(e) (Look-Back) GMRES(m), ℓ = 30(f) (Look-Back) FOM(m), ℓ = 30

0

0

-2

-2

log10 | P(λ) |log10 | P(λ) |0.8(b) (Look-Back) FOM(m), ℓ = 10

log10 | P(λ) |log10 | P(λ) |(a) (Look-Back) GMRES(m), ℓ = 10

-4
-6
-8
-10

-4
-6
-8
-10

-12

-12
0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

λ

0.6

0.8

1

λ

(g) (Look-Back) GMRES(m), ℓ = 40(h) (Look-Back) FOM(m), ℓ = 40

0

0

-2

-2

log10 | P(λ) |log10 | P(λ) |

0.6
λ

-4
-6
-8
-10

-4
-6
-8
-10

-12

-12
0

0.2

0.4

0.6

0.8

λ

(i) (Look-Back) GMRES(m), ℓ = 50

1

0

0.2

0.4

0.6

0.8

1

λ

(j) (Look-Back) FOM(m), ℓ = 50

LB−G
G
F
Fig. 4.1: The graphs of log10 |Pm×ℓ(λi )|, log10 |Pm×ℓ(λi )|, log10 |Pm×ℓ(λi )| and
LB−F
log10 |Pm×ℓ (λi )| for the restart cycle ℓ = 10, 20, . . . , 50.

0

GMRES(m)Look-Back GMRES(m)

-2
-4
-6
-8
-10
-12
0

200

400

600

800 1000 1200 1400

Number of restart cycle(a) (Look-Back) GMRES(m)log10 of relative residual 2-norm

log10 of relative residual 2-norm

LOOK-BACK-TYPE RESTART FOR KRYLOV SUBSPACE METHODS

0

11

FOM(m)
Look-Back FOM(m)

-2
-4
-6
-8
-10
-12
0

200

400

600

800 1000 1200 1400

Number of restart cycle(b) (Look-Back) FOM(m)Fig. 4.2: The relative residual 2-norm history of the GMRES(m) method, theFOM(m) method, Look-Back GMRES(m) method and the Look-Back FOM(m)method.

Next we consider the residual polynomials and the convergence rate of the LookBack GMRES(m) method and the Look-Back FOM(m) method. Unlike the nonuniform convergence manner for the GMRES(m) method and the FOM(m) method,
LB−G
the Look-Back-type restart reconstructs the residual polynomial such that |Pm×ℓ(λi )|
LB−F
and |Pm×ℓ (λi )| show almost the uniform convergence manner in all eigenvalues; see
black square ■ in Fig. 4.1. As the results of the uniform convergence manner ofthe residual polynomials, the Look-Back GMRES(m) method and the Look-Back
FOM(m) method have the better convergence for the model problem (4.1) than the
GMRES(m) method and the FOM(m) method shown in Fig. 4.2.
We note that we also have almost the same results for other symmetric positivedeﬁnite linear systems with the diﬀerent eigenvalue distribution.
5. Numerical experiments and results. In this section, we test the performance of the GMRES(m) method and the Look-Back GMRES(m) method with nopreconditioners in Section 5.1 and with the ILU(0) preconditioner [10] in Section 5.2
respectively. We also examine, in Section 5.3, the relationship between the symmetry
property of the coeﬃcient matrix and the convergence rate of these methods.
All the numerical experiments were carried out in double precision arithmetic
on OS: CentOS 64bit, CPU: Intel Xeon X5550 2.67GHz (1 core), Memory: 48GB,
Compiler: GNU Fortran ver. 4.1.2, Compile option: -O3.
5.1. Numerical experiments I. In this section, we test the performance of the
GMRES(m) method and the Look-Back GMRES(m) method without preconditioners. The performance of these methods is evaluated by the test problems from The
University of Florida Sparse Matrix Collection [7].
The characteristics of the coeﬃcient matrices of the test problems are shown in
Table 5.1. (R) or (C) denotes the matrix type: Real nonsymmetric or Complex nonHermitian respectively. n, N nz and Ave.N nz denote the number of dimension, thenumber of nonzero elements and the average nonzero elements per row or column
respectively.
(1)
We set b = [1, 1, . . . , 1]T as the right-hand side, x0 = [0, 0, . . . , 0]T for the initial
guess. We also set m = 30, 100 as the restart frequency, d = 3 as the parameter for
the Look-Back-type restart and stopping criterion was set as ∥rk ∥2 /∥b∥2 ≤ 10−10 .

12

A. IMAKURA, T. SOGABE AND S.-L. ZHANG

Table 5.1: Characteristics of the coeﬃcient matrices of the test problems for the
GMRES(m) method and the Look-Back GMRES(m) method.
(Type) Matrix name(R) CAVITY10(R) CHIPCOOL0(R) COUPLED(R) EPB1(R) EX28(R) FLOWMETER5(C) KIM1(C) LIGHT IN TISSUE(R) NS3DA(R) RAEFSKY1(R) RAJAT03(R) RDB5000(C) WAVEGUIDE3D(R) XENON1(R) XENON2(C) YOUNG1C

n
2597
20082
11341
14734
2603
9669
38415
29282
20414
3242
7602
5000
21036
48600
157464
841

N nz
76367
281150
98523
95053
77781
67391
933195406084167959929427632653
29600
3034681181120
3866688
4089

Ave.N nz

29.41
14.00
8.69
6.45
29.88
6.97
24.29
13.87
82.28
90.77
4.30
5.92
14.43
24.30
24.56
4.86

Application area
Computational fluid dynamics
Model reduction problem
Circuit simulation
Thermal problem
Computational fluid dynamics
Model reduction problem
2D/3D problem
Electromagnetics problem
Computational fluid dynamics
Computational fluid dynamics
Circuit simulation
Computational fluid dynamics
Electromagnetics problem
Materials problem
Materials problem
Acoustics problem

Numerical results. The numerical results are presented in Tables 5.2–5.3. In
these tables, a † denotes that the methods did not converge within 100000 iterations.
We also present in Fig. 5.1 the relative residual 2-norm histories of both methods with
m = 30 for KIM1, LIGHT IN TISSUE, NS3DA, RAJAT03, RDB5000 and XENON2.
We analyze the results in terms of three aspects: convergence rate; computation time
per restart cycle (m iterations); and total computation time.
We ﬁrstly consider the convergence rate of both methods. In terms of the number of iterations (Iter), the Look-Back GMRES(m) method shows almost the sameor lower Iter than the GMRES(m) method in most cases. Especially, for CAVITY10(m = 30), CHIPCOOL0 (m = 30, 100), COUPLED(m = 30), EX28 (m = 30, 100),
FLOWMETER5 (m = 30, 100) and RAJAT03 (m = 30), the GMRES(m) method did
not converge within 100000 iterations; on the other hand, the Look-Back GMRES(m)method converged to the solution satisfying the required accuracy ∥rk ∥2 /∥b∥2 ≤10−10 ; see TRR of Tables 5.2–5.3. We can see from these results that the LookBack-type restart has high potential to improve signiﬁcantly the convergence rate of
the GMRES(m) method. We can also see from the comparison between the numerical results for m = 30, 100 that the smaller restart frequency m leads to the larger
improvements in terms of Iter.
We can also see from Fig. 5.1 that the Look-Back GMRES(m) method shows themonotonic decrease in the residual 2-norm as well as the GMRES(m) method. Thismeans that Proposition 3.1 is experimentally supported by these results. For KIM1,
NS3DA and RDB5000, both methods show the same level of convergence throughout
the whole iteration; see Fig. 5.1 (a), (c), (e). On the other hand, the Look-Back-typerestart well played for LIGHT IN TISSUE, RAJAT03 and XENON2, and then the
Look-Back GMRES(m) method shows a very well convergence property comparedwith the GMRES(m) method throughout the whole iteration; see Fig. 5.1 (b), (d),(f).
Next, we consider the computation time per restart cycle (m iterations); seeLOOK-BACK-TYPE RESTART FOR KRYLOV SUBSPACE METHODS

13

Table 5.2: Convergence results (Iter : number of iterations, tTotal : total computation
time, tRestart : computation time per restart cycle, TRR : explicitly computed relative
residual 2-norm) of the GMRES(m) method and the Look-Back GMRES(m) method
for m = 30.
MatrixCAVITY10
CHIPCOOL0
COUPLED
EPB1
EX28
FLOWMETER5KIM1
LIGHT IN TISSUE
NS3DA
RAEFSKY1
RAJAT03
RDB5000
WAVEGUIDE3D
XENON1XENON2
YOUNG1C

MethodIter

GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)

†
27961
†
34321
†
78392
2826
1993
†
25321
†
78091
2820
3963
2964
938
2330
2317
5036
2551
†
16621
962
1021
34991
29685
11521
1879
15691
2372
6026
5564

Time[sec.]
ttotaltrestart†
7.79 × 10−307.40 × 10
7.94 × 10−3†
5.00 × 10−2
1
5.80 × 10
5.05 × 10−2†
2.43 × 10−2
1
6.43 × 10
2.46 × 10−202.53 × 10
2.69 × 10−201.85 × 10
2.81 × 10−2†
7.85 × 10−306.93 × 10
8.23 × 10−3†
1.78 × 10−2
1
4.68 × 10
1.80 × 10−21
2.79 × 10
3.02 × 10−1
1
4.02 × 10
3.06 × 10−111.80 × 10
1.82 × 10−105.81 × 10
1.86 × 10−1
11.49 × 10
1.92 × 10−1
1
1.51 × 10
1.96 × 10−103.18 × 10
1.89 × 10−201.62 × 10
1.90 × 10−2†
1.30 × 10−207.41 × 10
1.32 × 10−2−1
2.98 × 10
9.27 × 10−3
−13.08 × 10
9.05 × 10−3
1.85 × 102
1.58 × 10−1
1.57 × 102
1.59 × 10−1
1
6.20 × 10
1.61 × 10−1
11.06 × 10
1.70 × 10−12
2.81 × 10
5.37 × 10−1
1
4.26 × 10
5.41 × 10−1
1
7.16 × 10
4.67 × 10−3
1
6.66 × 10
4.33 × 10−3

TRR
-9.37
-10.00-0.43-10.00-6.36-10.00-10.00-10.00-0.82-10.01-0.38-10.00-10.01-10.01-10.00-10.00-10.00-10.00-10.00-10.04-0.55-10.02-10.00-10.02-10.00-10.00-10.00-10.00-10.00-10.00-10.00-10.00tRestart of Tables 5.2–5.3. In terms of computation time per restart cycle, the LookBack GMRES(m) method requires at most 10% more time than the GMRES(m)method for both restart frequency m, and the diﬀerence becomes smaller with increasing the restart frequency m. This depends on the fact that, in our implementation, the required additional operations for the Look-Back-type restart are just one
matrix-vector multiplication and few AXPY and inner-products per restart cycle; see
Table 3.1.

14

A. IMAKURA, T. SOGABE AND S.-L. ZHANG

Table 5.3: Convergence results (Iter : number of iterations, tTotal : total computation
time, tRestart : computation time per restart cycle, TRR : explicitly computed relative
residual 2-norm) of the GMRES(m) method and the Look-Back GMRES(m) method
for m = 100.
MatrixCAVITY10
CHIPCOOL0
COUPLED
EPB1
EX28
FLOWMETER5KIM1
LIGHT IN TISSUE
NS3DA
RAEFSKY1
RAJAT03
RDB5000
WAVEGUIDE3D
XENON1XENON2
YOUNG1C

MethodIter

GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)GMRES(m)LB-GMRES(m)

56098
5001
†
25301
13333
6115
1781
1479
†
19601
†
62901
519
514
888
828
1983
1959
2775
1703
93476
5304
387
380
25845
23826
3888
1912
4933
2273
1636
1493

Time[sec.]
ttotaltrestart
12.56 × 10
4.55 × 10−202.26 × 10
4.52 × 10−2†
3.56 × 10−1
1
9.10 × 10
3.58 × 10−11
2.36 × 10
1.77 × 10−1
1
1.09 × 10
1.78 × 10−103.91 × 10
2.21 × 10−103.24 × 10
2.21 × 10−1†
4.63 × 10−208.96 × 10
4.57 × 10−2†
1.42 × 10−1
1
8.67 × 10
1.38 × 10−109.38 × 10
1.85 × 100
09.29 × 10
1.83 × 100
11.10 × 10
1.24 × 100
11.02 × 10
1.25 × 100
11.64 × 10
8.26 × 10−1
1
1.58 × 10
8.08 × 10−102.41 × 10
8.72 × 10−201.49 × 10
8.71 × 10−2
1
9.56 × 10
1.02 × 10−105.43 × 10
1.03 × 10−1
−1
2.53 × 10
6.55 × 10−2−1
2.48 × 10
6.54 × 10−2
2.58 × 102
9.99 × 10−1
2.39 × 102
1.00 × 100
13.84 × 10
9.91 × 10−111.90 × 10
1.00 × 100
21.58 × 10
3.20 × 100
17.28 × 10
3.24 × 100
14.59 × 10
2.80 × 10−2
1
4.20 × 10
2.80 × 10−2

TRR
-10.00-10.01-1.62-10.00-10.00-10.00-10.00-10.01-2.19-10.02-6.08-10.00-10.03-10.00-10.02-10.00-10.01-10.00-10.01-10.01-10.00-10.00-10.00-10.03-10.00-10.00-10.00-10.00-10.00-10.00-10.00-10.00In terms of the total computation time (tTotal ), from the results of the smaller
Iter and almost the same tRestart , we can see that the Look-Back GMRES(m) method
can obtain the solution within much smaller computation time than the GMRES(m)method except the case for KIM1 (m = 30), NS3DA (m = 30) and RDB5000 (m =30).
5.2. Numerical experiments II. In this section, we test the performance of
the GMRES(m) method and the Look-Back GMRES(m) method with the ILU(0)150log10 of relative residual 2-norm

log10 of relative residual 2-norm

LOOK-BACK-TYPE RESTART FOR KRYLOV SUBSPACE METHODSGMRES(m)Look-Back GMRES(m)

-2
-4
-6
-8
-10
0

1000

2000

3000

0

GMRES(m)Look-Back GMRES(m)

-2
-4
-6
-8
-10
0

4000

-4
-6
-8
-10
0

log10 of relative residual 2-norm

log10 of relative residual 2-norm

GMRES(m)Look-Back GMRES(m)

-2

0

500

1000

1500

2000

-6
-8

GMRES(m)Look-Back GMRES(m)

-10
0

20000

40000

60000

80000 100000(d) RAJAT03

-8
-10
600

3000

-4

(c) NS3DA

400

2500

-2

Number of Iterations

-6

200

2000

0

2500

-4

0

1500

Number of Iterations

GMRES(m)Look-Back GMRES(m)

-2

1000

(b) LIGHT IN TISSUE

log10 of relative residual 2-norm

log10 of relative residual 2-norm(a) KIM1
0500

Number of Iterations

Number of Iterations

800

1000

1200

0

GMRES(m)Look-Back GMRES(m)

-2
-4
-6
-8
-10
0

5000

10000

15000

Number of Iterations

Number of Iterations(e) RDB5000(f) XENON2

20000

Fig. 5.1: The relative residual 2-norm history of GMRES(m) and Look-Back
GMRES(m) of m = 30 without preconditioners for KIM1, LIGHT IN TISSUE,
NS3DA, RAJAT03, RDB5000 and XENON2.

preconditioner, where the right preconditioning was used. The performance of thesemethods is evaluated by the test problems from the discretization of two types ofPartial Diﬀerential Equations (PDEs).
For the ﬁrst example, the test problems are obtained from the discretization of
the following PDE:−(Aux )x − (Auy )y + α exp(2(x2 + y 2 ))ux = F (x, y),with α = 0.0, 0.5, 1.0, 2.0, over the unit square (x, y) ∈ [0, 1] × [0, 1] with the Dirichletboundary condition. The boundary condition and functions A and F are shown in16

A. IMAKURA, T. SOGABE AND S.-L. ZHANG

Case 1

u= 0

u= 1

u= 1

F=102

A=104

A=10-5

A=102

u= 1

Fig. 5.2: The boundary condition and the functions A and F of van der Vorst’s PDE
for numerical experiments II. (case 1)

Fig. 5.2. These problems were taken from van der Vorst’s paper [20], and some authors
have also used to test the performance of the Krylov subspace methods, e.g., [18]. We
discretized them by ﬁve point central diﬀerences over the square grid with size 1/201
in each directions, and the obtained linear systems have 2002 unknowns and 199598
nonzero elements.
For the second example, the test problems are obtained from the discretizationof the Helmholtz equation [5]uxx + uyy + σ 2 u = 0
with σ = 0.69, 1.39, over the square region (x, y) ∈ [0, π] × [0, π]. The boundary
conditions are deﬁned by()1/2
ux |x=0 = i σ 2 − 41
cos y2 ,( 2 1 )1/2
ux − i σ − 4u|x=π = 0,
uy |y=0 = 0,u|y=π = 0,

(Neumann condition)(Radiation condition)(Neumann condition)(Dirichlet condition)

where i denotes i2 = −1. Note that the exact solution is written by u(x, y) =exp{i(σ 2 − 14 )1/2 x} cos y2 . We discretized them by ﬁve point central diﬀerences overthe square grid with size 1/201 in each directions, and the obtained linear systems
have (200 × 201) unknowns and 200596 nonzero elements.
(1)
For ﬁrst and second examples, we set x0 = [0, 0, . . . , 0]T for the initial guess, m =30 as the restart frequency, d = 3 as the parameter for the Look-Back GMRES(m)method and stopping criterion was set as ∥rk ∥2 /∥b∥2 ≤ 10−12 .
Numerical results. The relative residual 2-norm histories of the ILU(0) preconditioned GMRES(m) method and the ILU(0) preconditioned Look-Back GMRES(m)method for the ﬁrst example are presented in Fig. 5.3 and for the second example areshown in Fig. 5.4.
We can see from Fig. 5.3 for the van der Vorst’s PDE that the residual 2-norm
histories of the ILU(0) preconditioned GMRES(m) method stagnated at ∥rk ∥2 /∥b∥2 ≈10−5 ; on the other hand, the ILU(0) preconditioned Look-Back GMRES(m) method

17

GMRES(m)Look-Back GMRES(m)

-4
-6
-8
-10
-12
0

log10 of relative residual 2-norm

log10 of relative residual 2-norm

0
-2

0

10000

20000

30000

40000

-4
-6
-8
-10
-12
0

10000

20000

30000

Number of Iterations(a) α = 0.0(b) α = 0.5

-4
-6
-8
-10
-12
0

GMRES(m)Look-Back GMRES(m)Number of Iterations

GMRES(m)Look-Back GMRES(m)

-2

0
-2

50000

log10 of relative residual 2-norm

log10 of relative residual 2-norm

LOOK-BACK-TYPE RESTART FOR KRYLOV SUBSPACE METHODS

10000

20000

30000

40000

50000

0

40000

50000

GMRES(m)Look-Back GMRES(m)

-2
-4
-6
-8
-10
-12
0

10000

20000

30000

Number of Iterations

Number of Iterations(c) α = 1.0(d) α = 2.0

40000

50000

Fig. 5.3: The relative residual 2-norm history of the GMRES(m) method and the
Look-Back GMRES(m) method with the ILU(0) preconditioner for van der Vorst’s
PDE (case 1) with α = 0.0, 0.5, 1.0, 2.0.

converged smoothly for all parameters α. These results show that the Look-Backtype restart has high potential to remedy the stagnation of the residual 2-norm of the
GMRES(m) method.
For the Helmholtz equations, the ILU(0) preconditioned Look-Back GMRES(m)method also shows much higher convergence rate than the ILU(0) preconditioned
GMRES(m) method; see Fig. 5.4. Especially for the case of σ = 1.39, the convergence
rate of the ILU(0) preconditioned Look-Back GMRES(m) method is more than twotimes faster as compared with the ILU(0) preconditioned GMRES(m) method.
5.3. Numerical experiments III. In this section, we examine the relationshipbetween the symmetry property of the coeﬃcient matrix and the convergence rate of
the GMRES(m) method and the Look-Back GMRES(m) method without preconditioner.
For the test problems, we recall van der Vorst’s PDE:−(Aux )x − (Auy )y + α exp(2(x2 + y 2 ))ux = F (x, y),with α = 0.0, 0.1, 0.2, . . . , 10.0, over the unit square (x, y) ∈ [0, 1] × [0, 1], where we
consider two type of the boundary conditions and functions A and F shown in Fig. 5.5.
The test problems are obtained from ﬁve point central diﬀerences over the square gridwith size 1/51 in each directions, and obtained linear systems have 502 unknowns andA. IMAKURA, T. SOGABE AND S.-L. ZHANG

0log10 of relative residual 2-norm

log10 of relative residual 2-norm

18

GMRES(m)Look-Back GMRES(m)

-2
-4
-6
-8
-10
-12
0

5000

10000

15000

20000

0

GMRES(m)Look-Back GMRES(m)

-2
-4
-6
-8
-10
-12
0

5000

10000

15000

Number of Iterations

Number of Iterations(a) σ = 0.69(b) σ = 1.39

20000

Fig. 5.4: The relative residual 2-norm history of the GMRES(m) method and the
Look-Back GMRES(m) method with the ILU(0) preconditioner for the Helmholtz
equation with σ = 0.69, 1.39.

Case 2

Case 3

u= 0

u= 1

F=1

u= 1

u= 1

F=5

u= 1

A=1
A=0.2

A=2
A=10

A=0.5

u= 0

A=20

u= 1

u= 1

Fig. 5.5: The boundary condition and the functions A and F of van der Vorst’s PDE
for the numerical experiments III (case 2 and case 3).

12398 nonzero elements. Notice that the obtained linear systems are real symmetricfor α = 0.0, and the larger α leads to the stronger non-symmetry for the obtained
linear systems.
(1)
For this experiment, we set x0 = [0, 0, . . . , 0]T for the initial guess, m = 10, 30
as the restart frequency, d = 3 as the parameter for Look-Back GMRES(m) andstopping criterion was set as ||rk ||2 /||b||2 ≤ 10−12 .
Numerical results. We introduce the index θ(A) for the symmetry property ofthe matrix A ∈ Rn×n deﬁned by
θ(A) :=||(A − AT )/2||2,||A||2where 0 ≤ θ(A) ≤ 1, θ(A) ∈ R. We note that θ(A) = 0 for the symmetric matrices:
A = AT and θ(A) = 1 for the skew-symmetric matrices: A = −AT .
For the numerical results, the graphs of the comparison between the symmetry19

LOOK-BACK-TYPE RESTART FOR KRYLOV SUBSPACE METHODS

Number of iterations

1000

800
600
400
200
full GMRES

0

Number of iterations

0

0.05

0.1

0.15

800
600
400
200

0

0.025

0.05

0.075

Index θ(A) of symmetry property

Index θ(A) of symmetry property(a) case 1: full GMRES method(b) case 2: full GMRES method

20000

20000

GMRES(10)Look-Back GMRES(10)

15000
10000
5000

GMRES(10)Look-Back GMRES(10)

10000
5000
0

0

0.05

0.1

0.15

0.2

0

0.025

0.05

0.075

Index θ(A) of symmetry property

Index θ(A) of symmetry property(c) case 1: GMRES(10) method(d) case 2: GMRES(10) method

7000

7000

Number of iterations

GMRES(30)
Look-Back GMRES(30)

6000

0.1

15000

0

Number of iterations

full GMRES

00.2

Number of iterations

Number of iterations

1000

5000
4000
3000
2000
1000
0

0.1

GMRES(30)
Look-Back GMRES(30)

6000
5000
4000
3000
2000
1000
0

0

0.05

0.1

0.15

0.2

0

0.025

0.05

0.075

Index θ(A) of symmetry property

Index θ(A) of symmetry property(e) case 1: GMRES(30) method(f) case 2: GMRES(30) method

0.1

Fig. 5.6: The graphs of the number of iterations versus the index θ(A) of the symmetry property for the full GMRES method , the GMRES(m) method and the
Look-Back GMRES(m) method for van der Vorst’s PDE (case 2 and case 3) with
α = 0, 0.1, 0.2, . . . , 10.0.

property θ(A) of the coeﬃcient matrix and the number of iterations of both methods
are presented in Fig. 5.6 (c)–(f). The white square □ denotes the results of the
GMRES(m) method and the black square ■ denotes the results of the Look-Back
GMRES(m) method respectively. We also show the numerical results for the full
GMRES method (without restart) in Fig. 5.6 (a), (b), for comparison. The left
column of Fig. 5.6 shows the results for the case 2 of the van der Vorst’s PDE andthe right column shows for the case 3 respectively.
It is generally expected that the smaller θ(A) (near symmetric matrices) leads

20

A. IMAKURA, T. SOGABE AND S.-L. ZHANG

to the better convergence rate. This is experimentally supported by the numericalresults of the full GMRES method shown in Fig. 5.6 (a), (b). The number of iterationsof the full GMRES method increases almost monotonically with increasing θ(A).
However, the convergence rate of the GMRES(m) method becomes worse especially for the case of the smaller θ(A) which leaded to the better convergence for thefull GMRES method; see □ of Fig. 5.6 (c)–(f). This means that the restart slows theconvergence rate especially for near symmetric cases.
On the other hand, we can obtained the results that the Look-Back-type restart
remedies this diﬃculty. The number of iterations of the Look-Back GMRES(m)method decreases in many cases especially for the smaller θ(A) as compared with
the GMRES(m) method, then the Look-Back GMRES(m) method shows the better
convergence with decreasing θ(A) as well as the full GMRES method; see ■ of Fig. 5.6(c)–(f).
6. Conclusion. In this paper, we investigated the restarted Krylov subspacemethods, as typiﬁed by the GMRES(m) method and the FOM(m) method, for solving
non-Hermitian linear systems. We proposed the Look-Back-type restart based on
the implicit residual polynomial reconstruction via the initial guess for the restarted
Krylov subspace method.
From our comparison analysis and numerical experiments, we leaned that the
Look-Back-type restart can reconstruct the residual polynomial uniformly and it modiﬁes the convergence property of the restarted Krylov subspace methods in many
problems. Especially for the near symmetric cases, the Look-Back-type restart significantly improves the convergence rate of the restarted Krylov subspace methods.
Acknowledgments. This work is supported in part by Strategic Programs forInnovative Research (SPIRE) Field 5 “The origin of matter and the universe” andKAKENHI (Grant No. 25870099).
REFERENCES[1] James Baglama, Daniela Calvetti, Gene H. Golub, and Lothar Reichel, Adaptively
preconditioned GMRES algorithms, SIAM Journal on Scientiﬁc Computing, 20 (1999),
pp. 243–269.
[2] James Baglama and Lothar Reichel, Augmented GMRES-type methods, Numerical Linear
Algebra with Applications, 14 (2007), pp. 337–350.
[3] Allison H. Baker, Elizabeth R. Jessup, and Tzanio V. Kolev, A simple strategy for varying the parameter in GMRES(m), Journal of Computational and Applied Mathematics,
230 (2009), pp. 751–761.
[4] Allison H. Baker, Elizabeth R. Jessup, and Thomas Manteuffel, A technique for accelerating the convergence of restarted GMRES, SIAM Journal on Matrix Analysis andApplications, 26 (2005), pp. 962–984.
[5] Alvin Bayliss, Charles I. Goldstein, and Eli Turkel, An iterative method for theHelmholtz equation, Journal of Computational Physics, 49 (1983), pp. 443–457.
[6] Kevin Burrage and Jocelyne Erhel, On the performance of various adaptive preconditioned
GMRES strategies, Numerical Linear Algebra with Applications, 5 (1998), pp. 101–121.
[7] Timothy A. Davis, The University of Florida Sparse Matrix Collection.
[8] Jocelyne Erhel, Kevin Burrage, and Bert Pohl, Restarted GMRES preconditioned bydeﬂation, Journal of Computational and Applied Mathematics, 69 (1996), pp. 303–318.
[9] Akira Imakura, Tomohiro Sogabe, and Shao-Liang Zhang, An eﬃcient variant of the
GMRES(m) method based on the error equations, East Asian Journal on Applied Mathematics, 2 (2012), pp. 1–22.
[10] J. A. Meijerink and Henk A. van der Vorst, An iterative solution method for linear systemsof which the coeﬃcient matrix is a symmetric M -matrix, Mathematics of Computation,
31 (1977), pp. 148–162.

LOOK-BACK-TYPE RESTART FOR KRYLOV SUBSPACE METHODS

21[11] Ronald B. Morgan, A restarted GMRES method augmented with eigenvectors, SIAM Journal
on Matrix Analysis and Applications, 16 (1995), pp. 1154–1171.
[12], Implicitly restarted GMRES and Arnoldi methods for nonsymmetric systems of equations, SIAM Journal on Matrix Analysis and Applications, 21 (2000), pp. 1112–1135.
[13], GMRES with deﬂated restarting, SIAM Journal on Scientiﬁc Computing, 24 (2002),
pp. 20–37.
[14] Kentaro Moriya and Takashi Nodera, The DEFLATED-GMRES(m, k) method withswitching the restart frequency dynamically, Numerical Linear Algebra with Applications,
7 (2000), pp. 569–584.
[15] Yousef Saad, Krylov subspace methods for solving large unsymmetric linear systems, Mathematics of Computation, 37 (1981), pp. 105–126.
[16], Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, PA, second ed., 2003.
[17] Yousef Saad and Martin H. Schultz, GMRES: a generalized minimal residual algorithm
for solving nonsymmetric linear systems, SIAM Journal on Scientiﬁc and Statistical Computing, 7 (1986), pp. 856–869.
[18] Gerard L. G. Sleijpen and Diederik R. Fokkema, BiCGStab(l) for linear equations involving unsymmetric matrices with complex spectrum, Electronic Transactions on Numerical
Analysis, 1 (1993), pp. 11–32.
[19] Maria Sosonkina, Layne T. Watson, Rakesh K. Kapania, and Homer F. Walker, A new
adaptive GMRES algorithm for achieving high accuracy, Numerical Linear Algebra with
Applications, 5 (1998), pp. 275–297.
[20] Henk A. van der Vorst, Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for
the solution of nonsymmetric linear systems, SIAM Journal on Scientiﬁc and StatisticalComputing, 13 (1992), pp. 631–644.
[21], Iterative Krylov Methods for Large Linear Systems, Cambridge University Press, New
York, NY, 2003.
[22] Lingie Zhang and Takashi Nodera, A new adaptive restart for GMRES(m) method, Australian & New Zealand Industrial and Applied Mathematics Journal, 46 (2005), pp. 409–
426.