Adaptively Preconditioned GMRES Algorithms.
The restarted GMRES algorithm proposed by Saad and Schultz [SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856--869] is one of the most popular iterative methods for the solution of large linear systems of equations Ax=b with a nonsymmetric and sparse matrix. This algorithm is particularly attractive when a good preconditioner is available. The present paper describes two new methods for determining preconditioners from spectral information gathered by the Arnoldi process during iterations by the restarted GMRES algorithm. These methods seek to determine an invariant subspace of the matrix A associated with eigenvalues close to the origin and to  move these eigenvalues so that a higher rate of convergence of the iterative methods is achieved.
Introduction
. Many problems in Applied Mathematics and Engineering give
rise to very large linear systems of equations
with a sparse nonsymmetric nonsingular matrix A. It is often desirable, and sometimes
necessary, to solve these systems by an iterative method. Let x 0 be an initial
approximate solution of (1.1), and let r be the associated residual vector.
Introduce the Krylov subspaces
associated with the matrix A and vector r 0 . Many popular iterative methods determine
the mth iterate, xm , so that xm \Gamma x We refer to such methods
as Krylov subspace iterative methods; see, e.g., Freund et al. [12] for a recent review.
Let the iterate xm be generated by a Krylov subspace iterative method. Then
the residual error r associated with xm satisfies
where the residual polynomial pm is determined by the iterative method, and satisfies
denote the Euclidean norm on R n , as well as the associated
induced matrix norm on R n\Thetan . The restarted Generalized Minimal Residual
algorithm by Saad and Schultz [22], described in Section 3, is one of
the most popular Krylov subspace iterative methods for the solution of linear systems
with a nonsymmetric matrix. The residual polynomial determined by this algorithm
satisfies
Department of Mathematics and Computer Science, Kent State University, Kent, OH 44242.
E-mail: jbaglama@mcs.kent.edu. Research supported in part by NSF grant F377 DMR-8920147
ALCOM.
y Department of Mathematical Sciences, Stevens Institute of Technology, Hoboken, NJ 07030.
E-mail: na.calvetti@na-net.ornl.gov. Research supported in part by NSF grant DMS-9404692.
z Computer Science Department, Stanford University, Stanford, CA 94305.
x Department of Mathematics and Computer Science, Kent State University, Kent, OH 44242.
E-mail: reichel@mcs.kent.edu. Research supported in part by NSF grant DMS-9404706.
J. Baglama et al.
m denotes the set of all polynomials p of degree at most m such that
The analysis and implementation of the restarted GMRES(m) algorithm, and
modifications thereof, continue to receive considerable attention; see, e.g., [4, 5, 7,
These algorithms are particularly attractive when a suitable
preconditioner n\Thetan for the matrix A is available; see, e.g., [2, 15, 21] for
recent discussions on preconditioners. A matrix M \Gamma1 is a good preconditioner if the
application of an iterative method of interest to the preconditioned linear system of
equations
gives a higher rate of convergence of the computed iterates than application of the
iterative method to the original linear system (1.1). Moreover, we would like the
preconditioner M \Gamma1 have the property that for any w 2 R n , the vector M \Gamma1 w can
be rapidly evaluated. The matrix M \Gamma1 in (1.5) is sometimes referred to as a left
preconditioner.
The present paper describes two new adaptive methods for determining preconditioners
during the iterations with the restarted GMRES(m) algorithm. The standard
implementation of the restarted GMRES(m) algorithm [22] is based on the Arnoldi
process [1], described in Section 2, and this allows spectral information of A to be
gathered during the iterations. We use this information to determine an approximation
of an invariant subspace of A associated with eigenvalues close to the origin.
Our preconditioner essentially removes the influence of these eigenvalues on the rate
of convergence. We will focus on the effect of the preconditioner on the spectrum
of A, however, it is known that the rate of convergence of the iterates computed by
the GMRES(m) algorithm also is determined by pseudospectra of A; see Nachtigal
et al. [19]. For ease of presentation, we ignore the effect of the preconditioner on
the pseudospectra in the present paper. Our preconditioners are particularly effective
when there is a cluster of a few eigenvalues of A that have a large influence on the rate
of convergence. A few illustrations can be found in Section 4. The determination as
well as the application of our preconditioners does not require the evaluation of any
matrix-vector products with the matrix A in addition to those needed for the Arnoldi
process and for the evaluation of certain residual errors. The implementations use
the recurrence formulas of the Implicitly Restarted Arnoldi (IRA) method described
by Sorensen [23] and more recently by Lehoucq [17]. Our preconditioners can be
combined with other preconditioners, and are also applicable when no other known
efficient preconditioner is available.
A different method to adaptively determine a preconditioner during iterations
by the restarted GMRES(m) algorithm has recently been described by Erhel et al.
[11]. By utilizing the recurrence formulas of the IRA method, our preconditioning
scheme allows more flexibility in the choice of preconditioner and requires less computer
memory than the method described in [11]. Another adaptive preconditioning
method has been presented by Kharchenko and Yeremin [16]. Their method differs
from our schemes in how approximate invariant subspaces are determined. Morgan
[18] also uses approximate invariant subspaces to improve the rate of convergence
of the restarted GMRES(m) algorithm; instead of constructing a preconditioner, he
appends an approximate invariant subspace to the Krylov subspaces generated by
the Arnoldi process. We feel that our new algorithms are attractive because of their
simplicity, and because the IRA method, on which our algorithms are based, typically
determines adequate approximate invariant subspaces fairly rapidly.
Adaptive preconditioners 3
Let the matrix A have the spectral factorization
Then (1.3) and (1.4) yield the bound
where -(A) denotes the spectrum of A. Note that the bound (1.7) would decrease
if we were able to replace -(A) by a subset. Our preconditioners have roughly this
effect. For definiteness, assume that the eigenvalues of A have been ordered according
to
and let A be scaled so that
A good approximation of such a scaling of A can be determined during the iterations.
This is discussed below.
The Arnoldi process determines a decomposition of the form
where Vm 2 R n\Thetam
m\Thetam is an upper
Hessenberg matrix. We refer to (1.10) as an Arnoldi decomposition of A. Throughout
this paper e j denotes the jth axis vector of appropriate dimension, and I j denotes
the identity matrix of order j. When Vm e the columns of Vm span the
Krylov subspace Km (A; r 0 ) defined by (1.2). For future reference, we define
Let the matrix V k 2 R n\Thetak consist of the first k columns of Vm , and let the
columns of the matrix Wn\Gammak span R n nspanfV k g, where spanfV k g denotes the span
of the columns of V k . Assume that W T
I n\Gammak . Thus, the columns of the
matrix [V k Wn\Gammak ] form an orthogonal basis of R n . Introduce the matrix
We will use the inverse of matrices of the form (1.12) with k - n as left precondi-
tioners. The form of the inverse is given below.
Proposition 1.1. Let Q 2 R n\Thetan be an orthogonal matrix, and partition it
according to where the submatrix V consists of the k first columns of Q,
and the submatrix W consists of the remaining columns. Assume that
nonsingular. Then the matrix
is nonsingular, and its inverse is given by
4 J. Baglama et al.
Proof. The matrix (1.13) can be written as
\Theta
I n\Gammak
and therefore
\Theta
I n\Gammak
This shows (1.14).
When the columns of the matrix V in Proposition 1.1 span an invariant subspace
of A, the eigenvalues of the matrix M \Gamma1 A can be expressed in terms of the eigenvalues
of A.
Corollary 1.2. Let the matrices V , W and H be as in Proposition 1.1, and
assume, moreover, that the columns of the matrix V span an invariant subspace of A
associated with the eigenvalues
where the eigenvalue 1 has multiplicity at least k.
Proof. The matrix A is similar to
~
A
\Theta

A 12
A 22
and -( ~
A 22 -ng. Formula (1.16) and the representation (1.15) yield
\Theta
I
A 12
A 22
Thus, the spectrum of M \Gamma1 A consists of -( ~
A 22 ) and the eigenvalue 1. The multiplicity
of the latter eigenvalue is at least k.
A result analogous to Corollary 1.2 for a right preconditioner is shown by Erhel et
al. [11]. We remark that application of preconditioners of the form (1.14) is simplified
by the fact that
Thus, the matrix W does not have to be computed.
The following example compares bounds for the rate of convergence of iterates
determined by the GMRES(m) algorithm when applied to the original linear system
(1.1) and to the preconditioned linear system (1.5) with the preconditioner (1.14),
where we assume that the conditions of Corollary 1.2 hold.
Example 1.1. Assume that A has a spectral factorization of the form (1.6) with
all eigenvalues real and positive, and let the eigenvalues be ordered according to (1.8).
Then (1.7) yields that
lim sup
min
Adaptive preconditioners 5
where Tm (z) is the Chebyshev polynomial of the first kind of degree m for the interval
and the equality (1.19) follows
from well-known properties of the Chebyshev polynomials; see, e.g., [13, Section
10.1.5].
be the preconditioner (1.14), and assume that the conditions of Corollary
1.2 hold. This preconditioner eliminates the influence of the k smallest eigen-values
of A on the rate of convergence of the GMRES(m) algorithm. Specifically,
the GMRES(m) algorithm applied to the preconditioned linear system (1.5) yields a
sequence of residual vectors that can be bounded by
lim sup
where, as usual, r . The bound (1.20) can be shown by first noting that
lim sup
and then applying the bound (1.18) to the right-hand side of (1.21). 2
In actual computations, we determine a preconditioner from a Krylov subspace
spanfV k g, which is close to an invariant subspace. The computations of Example
1.1 suggest that the GMRES(m) algorithm will require fewer iterations to determine
an accurate approximate solution of (1.1) when applied to the preconditioned linear
system (1.5) with such a preconditioner than when applied to the original unpreconditioned
system (1.1). This is verified by numerical experiments, some of which are
presented in Section 4.
2. Construction of the preconditioners. In this section we describe how to
determine an approximate invariant subspace of A associated with the eigenvalues
by using the recursion formulas of the IRA method of Sorensen [23].
We first present the Arnoldi process [1].
Algorithm 2.1. Arnoldi process
Input: k, m, upper Hessenberg matrix H
Output: upper Hessenberg matrix
n\Thetam ,
do
do h 'j := v T
endfor j;We may assume that all vectors f generated by Algorithm 2.1 are
nonvanishing, because if f then the columns of the matrix V j generated span an
invariant subspace of A, and V j can be used to construct a preconditioner as described
in Example 1.1. When on input to Algorithm 2.1, only the initial vector f 0 has
6 J. Baglama et al.
to be provided. We note that if f m 6= 0, then we can define the matrix
with orthogonal columns. In the sequel, we will use the matrix
This is an (m matrix of Hessenberg-type, whose leading m \Theta m principal
submatrix is Hm , and whose (m 1)st row is
.
Numerical difficulties can arise in the computation of the vectors f j in the al-
gorithm. The computations are done by the modified Gram-Schmidt process, with
one reorthogonalization. Neglecting to enforce orthogonality of each vector f j against
the vectors give rise to spurious eigenvalues of the matrix Hm , i.e.,
eigenvalues which cannot be considered approximations of eigenvalues of A.
Given the Arnoldi decomposition (1.10) with initial vector the recursion
formulas of the IRA method can be used to compute the vector
for any monic polynomial
of degree evaluating any new matrix-vector products with the matrix
A. The coefficient jm\Gammak is a scaling factor chosen so that kv (m\Gammak)
1. We will
discuss the selection of the zeros z j below.
The recursion formulas of the IRA method are truncated versions of the recursion
formulas for the QR algorithm with explicit shifts, with the zeros z j chosen as
shifts; see, e.g., [13, Chapter 7] for a description of the QR algorithm. We therefore
sometimes refer to the zeros z j as shifts. Thus, let the decomposition (1.10) be
given, and determine the QR factorization
I m and R (1) is upper triangular. Putting
is also a Hessenberg matrix.
Multiplication of equation (2:4:2) by e 1 yields
where ae (1)
Equation (2.5) displays the relationship
between the initial Arnoldi vector v 1 and the vector v (1)
1 . After applying shifts
we obtain the decomposition
(2.
Adaptive preconditioners 7
where
Here denotes the orthogonal matrix associated
with the shift z j . Introduce the partitioning
Hm\Gammak
and equate the first k columns on the right-hand side and left-hand side of (2.6). This
gives
where
and V (m\Gammak)
m\Gammak ]. It follows from
and (2.8), that (V (m\Gammak)
Thus, (2.7) is an Arnoldi decomposition of A.
By construction, the vector v (m\Gammak)
can be written as (2.3).
While our description of the IRA method is based on recursion formulas for the
QR algorithm with explicit shifts, our implementation is based on the QR algorithm
with implicit shifts for reason of numerical stability; see [13, Chapter 7] for a description
of this QR algorithm. The use of implicit shifts allows the application of complex
conjugate shifts without using complex arithmetic.
We first apply the Arnoldi process to compute the Arnoldi decomposition (1.10),
and then use the recursion formulas of the IRA method to determine the Arnoldi
decomposition (2.7). The purpose of these computations is to determine an accurate
approximation of an invariant subspace of A associated with the eigenvalues
We would like to choose the zeros z z m\Gammak of /m\Gammak , so that
the first column v (m\Gammak)
k defined by (2.3) is in, or close to, an invariant
subspace of A associated with the eigenvalues
Let f' (m)
denote the eigenvalues of the upper Hessenberg matrix Hm in
(1.10), and order them so that
m AVm is an orthogonal projection of A, we consider the ' (m)
j to be
approximations of eigenvalues of A. In order to force the vector v (m\Gammak)
1 into an
invariant subspace of A associated with the k eigenvalues of A of smallest magnitude,
we choose the zeros
i.e., the z j are chosen to be available approximations of the eigenvalues of A
of largest magnitude. This selection of zeros is discussed by Sorensen [23], Calvetti
8 J. Baglama et al.
et al. [6] and Lehoucq [17], and these zeros are there referred to as "exact shifts".
Numerical experience indicates that the ordering of the shifts according to (2.10) is
adequate in the sense that the computed matrices H (m\Gammak)
very close to the k eigenvalues of Hm of smallest magnitude.
Let f' (k)
j=1 be eigenvalue-eigenvector pairs of H (m\Gammak)
k , and introduce the
vectors
are approximate eigenvalue-eigenvector
pairs of A with residual errors
We accept spanfV k g as an approximate invariant subspace of A if
kf (m\Gammak)
where ffl subspace ? 0 is a parameter. The purpose of the matrix H (m\Gammak)
k in (2.11) is to
make the bound invariant under scaling of A.
If the inequalities (2.11) are not satisfied, then we apply Algorithm 2.1 with the
Arnoldi decomposition (2.7) as input in order to determine a new Arnoldi decomposition
(1.10) with an m \Theta m upper Hessenberg matrix Hm . We then again apply the
recursion formulas of the IRA method with the zeros chosen to be the
of Hm of largest magnitude. This gives an Arnoldi decomposition of the form
(2.7), and we check whether the inequalities (2.11) are satisfied. The computations
are repeated in this fashion until (2.11) holds. We obtain in this way an Arnoldi
decomposition of the form (2.7) with matrices
k and H
k , such
that, in general, spanfV k g is an accurate approximation of an invariant subspace associated
with the k eigenvalues of smallest magnitude of A, and -(H k ) is an accurate
approximation of the set f- j g k
. The accuracy of the approximations depends on
the parameter ffl subspace in (2.11), the distribution of the eigenvalues of A, and the
departure from normality of A. The matrices V k and H k so obtained are used to
define our first preconditioner
where we have used (1.17).
We describe in Section 3 how to combine the IRA process with the restarted
GMRES algorithm and Richardson iteration, so that we can improve an available
approximate solution of (1.1) while determining the preconditioner M \Gamma1
1 .
Having computed the preconditioner M
1 , we apply the method outlined
above to the preconditioned system (1.5) in order to determine an approximation of an
invariant subspace associated with the eigenvalues of smallest magnitude of the matrix
simultaneously improve an available approximate solution of (1.1). This
yields a new preconditioner M \Gamma1
2 for the system M \Gamma1
equivalently,
a new preconditioner M
1 for the system (1.1). The computations are
continued in this manner until we have determined a preconditioner of the form
for some specified integer ff 0 - 1. The form (2.13) of the preconditioner makes it
natural to scale A so that (1.9) holds. An approximation of such a scaling is achieved
by scaling the linear system (1.1) by the factor 1=j' (m)
m is an eigenvalue
Adaptive preconditioners 9
of largest magnitude of one of the matrices Hm generated by Algorithm 2.1 during
our computation of the preconditioner M \Gamma1
1 . We remark that for certain matrices A
other techniques for achieving such a scaling may be available. For instance, one may
be able to use Gershgorin disks or the inequality
matrix norm induced by a vector norm; see [24, Chapter 6] for details on the latter
topics.
3. The iterative methods. This section describes our two algorithms for adaptive
preconditioning in detail. One of them, Algorithm 3.5, combines the IRA process
with Richardson iteration and the GMRES algorithm. The other scheme, Algorithm
3.6, does not apply Richardson iteration. We first recall the restarted GMRES(m)
algorithm by Saad and Schultz [22] for the solution of linear systems of equations
(1.1).
Algorithm 3.1. Restarted GMRES(m) algorithm
Input: m, initial approximate solution x
Output: approximate solution xm , associated residual vector r m .
while krm k=kr solution do
Compute by Algorithm 2.1 with input
. Then the matrices Vm+1 and -
Hm , defined by (2.1) and (2.2),
respectively, are also available.
Compute solution ym
Hm yk.
endwhile;We now describe how to improve an available approximate solution by Richardson
iteration while applying the recursion formulas of the IRA method to an Arnoldi
decomposition. These iterations can be carried out without evaluating matrix-vector
products with the matrix A. Let x 0 be an available approximate solution of (1.1).
Richardson iteration can be written as
where the are relaxation parameters. We would like the parameters ffi j to be
such that the approximate solutions x j converge rapidly to the solution of (1.1) as j
increases. For future reference, we note that the residual vectors (3.2) can be written
as
Y
Theorem 3.2. Let x 0 be an approximate solution of (1.1), and let r
Consider the Arnoldi decomposition
m with the initial vector
k. Apply the recursion formulas of the IRA method
with zeros z m. Then the residual vectors (3.2) associated
with the iterates (3.1) computed by Richardson iteration with relaxation parameters
J. Baglama et al.
are given by
Y
m. Here Q (') denotes
the orthogonal matrix and R (') the upper triangular matrix associated with the zero
z ' in the IRA recursion formulas. Moreover, v (m)
Proof. We first show (3.4) for 1. Substitution of v
The representation (3.3) now shows that
We turn to the case when 2. From (3.3) and (3.5), we obtain
Replace by VmQ (1) in equations (2.4.1)-(2.4.4), and multiply the equation
(2.4.2) so obtained by e 1 . This shows, analogously to (2.5), that
Substitution of (3.7) into (3.6) shows (3.4) for 2. Continuing in this manner yields
m.
The case has to be treated separately. We have the Arnoldi decomposition
and similarly as in [3], we obtain v (m)
1 . Choosing
ae (m)
completes the proof.
Prior to the development of the GMRES(m) algorithm, Saad [20] introduced the
Full Orthogonalization algorithm. This is a Galerkin method for the solution
of (1.1). Let x 0 be an approximate solution of (1.1) and let r 0 be the associated
residual vector. Consider the Arnoldi decomposition (1.10), and let v be the
same as in Theorem 3.2. The FO(m) algorithm determines an improved approximate
solution xm of (1.1) by solving the linear system
and then letting
The following result shows how this approximate solution can be determined by
Richardson iteration.
Theorem 3.3. Let the vectors x 0 and r 0 , and the Arnoldi decomposition (1.10),
be the same as in Theorem 3.2. Assume that the Arnoldi decomposition exists with
Adaptive preconditioners 11
kfm k 6= 0 and that the matrix Hm in the Arnoldi decomposition is nonsingular. Let
in (3.4), and let the relaxation parameters for Richardson iteration be reciprocal
values of the eigenvalues of Hm . Then, in exact arithmetic, the approximate solution
xm determined by Richardson iteration (3.1)-(3.2) equals the approximate solution
computed by the FO(m) algorithm.
Proof. Substitute use the fact that the linear system
(3.8) can be written as Hm
Introduce, for polynomials f and g, the bilinear form
By construction,
where g j is a polynomial of degree j. The
In particular, equations (3.10) and (3.11) yield
which shows that pm
is the residual polynomial of degree m
for the FO(m) algorithm, and therefore satisfies pm Combining formulas
(1.10) and (3.11) yields the identity
which shows that the eigenvalues f' (m)
of Hm are the zeros of gm . In particular,
all ' (m)
therefore pm can be written as pm It follows that
Y
A comparison of (3.12) with (1.3) and (3.3) shows that m steps of Richardson iteration
with relaxation parameters
j and an application of the FO(m) algorithm
correspond to the same residual polynomial, and therefore are equivalent.
The implementation of our iterative method is based on the following observation.
Corollary 3.4. Let x j \Gamma1 be an approximate solution of (1.1), and let r
be the associated residual vector. Let
be an Arnoldi decomposition, with initial vector
j=1 the eigenvalues of H ' , let x be the approximate solution
J. Baglama et al.
obtained by one step of Richardson iteration with relaxation parameter
q , for
let an application of the recursion formulas of the IRA method
to (3.13) with shift ' (')
q yield the Arnoldi decomposition
AV (1)
is the triangular matrix in a QR factorization of H
I ' . Moreover, -(H (1)
.
Proof. The Corollary follows from Theorem 3.2 and the fact that when we use
an exact shift, the eigenvalues of the reduced matrix H (1)
are the eigenvalues of the
original matrix H ' , except for the shift. The latter result is shown by Sorensen [23,
Lemma 3.10].
The corollary above shows that we can apply shifts, one at a time, and
determine the required residual vectors from the first column of the matrices V ' in
the available Arnoldi decompositions. An analogous result can be established for
complex conjugate shifts. In the latter case, the recursion formulas for the IRA
method are implemented by using the QR algorithm with implicit double shifts. This
obviates the need to use complex arithmetic. A double step of Richardson iteration,
with complex conjugate relaxation parameters, also can be carried out without using
complex arithmetic. For notational simplicity, the algorithm below for our iterative
method does not use double shifts and double steps, however, our implementation of
the algorithm used for the computed examples of Section 4 does.
Algorithm 3.5. Adaptively preconditioned GMRES(m) algorithm
with Richardson iteration
Input: tolerance for computed approximate solution ffl solution , tolerance for approximate
invariant subspace ffl subspace , dimension m of largest Krylov subspace
determined, dimension k of approximate invariant subspace to be computed, maximum
number ff 0 of preconditioners M \Gamma1
j to be computed, maximum number fi 0
of Arnoldi decompositions of order m to be determined for each preconditioner
.
Output: computed approximate solution x associated residual vector r j , preconditioner
for do
Compute
m by Algorithm
2.1 with initial vector
for do
Compute eigenvalues f' (m)
of matrix Hm in Arnoldi decomposition
and order them according to (2.9).
scale matrix and right-hand side of linear system by
factor 1=j' (m)
j. Then equation (1.9) is approximately satisfied.
do
Apply shift ' (m)
m+1\Gamma' to Arnoldi decomposition and compute
residual vector M \Gamma1 r j as described by Corollary 3.4. This gives
m\Gamma' .
Adaptive preconditioners 13
endfor ';
if bound (2.11) is satisfied then goto
Use the Arnoldi decomposition M \Gamma1 AV (m\Gammak)
k as input to Algorithm 2.1 and apply the Arnoldi process to
compute the Arnoldi decomposition M
.
endfor fi;
1: Improve approximate solution by GMRES(k) and update preconditioner:
The
k , as well as the matrices
and (2.2), are available.
Compute solution y k 2 R k of min
2: M \Gamma1
3: r j+k :=
solution then done;
endfor ff;
while kr j k=kr solution do
Apply
M by Algorithm 2.1 with initial vector
the matrices Vm+1 and -
Hm defined by (2.1) and (2.2), respectively.
Compute solution ym
4: r j+m :=
endwhile;In Algorithm 3.5, we only have to compute matrix-vector products with the matrix
A when applying the Arnoldi process, and when evaluating residual vectors r ' in the
lines labeled "3:" and "4:".
We now examine the storage requirement of Algorithm 3.5 and count the number
of n-vectors that have to be stored. Storage necessary to represent the matrix A is
ignored, since it is independent of the iterative method used. Each preconditioner
requires the storage of an n \Theta k matrix V k , and we limit the number of these
preconditioners to ff 0 . Thus, the preconditioner M \Gamma1 defined by (2.13) requires the
storage of at most ff 0 k n-vectors. In particular, the matrix M \Gamma1 is not actually
formed. The line marked "2:" in Algorithm 3.5 is to be interpreted symbolically to
mean that the storage for the matrix M \Gamma1 and the formula for evaluating matrix-vector
products with M \Gamma1 are updated. The GMRES(m) algorithm in the while-loop
of Algorithm 3.5 requires additional storage for the vectors x j and r j , and for the
matrix Vm+1 2 R n\Theta(m+1) . This is equivalent to the storage of m+ 3 n-vectors. The
vector Algorithm 3.5 is up to a scaling factor stored in the first column of
the matrix Vm+1 . The last column of Vm+1 contains the vector fm up to a scaling
factor. The right-hand side vector b also has to be stored. Therefore, the total storage
requirement of Algorithm 3.5 is at most ff
Algorithm 3.6 below is obtained by replacing Richardson iteration in Algorithm
3.5 by the GMRES algorithm. This replacement makes the the residual error decrease
more smoothly as the iterations proceed. However, the iterates and preconditioners
generated by Algorithms 3.5 and 3.6 are not the same, and we have found that the
14 J. Baglama et al.
former algorithm not seldom gives faster convergence. This is illustrated in Section
4. We therefore feel that both algorithms are of interest. The storage requirement
of Algorithm 3.6 is essentially the same as of Algorithm 3.5. For notational simplic-
ity, Algorithm 3.6 does not use double shifts, however, our implementation of the
algorithm used for the computed examples of Section 4 does.
Algorithm 3.6. Adaptively preconditioned GMRES(m) algorithm
Input: tolerance for computed approximate solution ffl solution , tolerance for approximate
invariant subspace ffl subspace , dimension m of largest Krylov subspace
determined, dimension k of approximate invariant subspace to be computed, maximum
number ff 0 of preconditioners M \Gamma1
j to be computed, maximum number fi 0
of Arnoldi decompositions of order m to be determined for each preconditioner
.
Output: computed approximate solution x associated residual vector r j , preconditioner
for do
Compute
m by Algorithm
2.1 with initial vector
Apply GMRES(m): determine the matrices Vm+1 and -
Hm defined by (2.1)
and (2.2), respectively.
Compute solution ym
for do
Compute eigenvalues f' (m)
of matrix Hm in Arnoldi decomposition
and order them according to (2.9).
scale matrix and right-hand side of linear system by
factor 1=j' (m)
j. Then equation (1.9) is approximately satisfied.
do
Apply shift ' (m)
m+1\Gamma' to Arnoldi decomposition by using the IRA
formulas (2.4)-(2.8). This gives Arnoldi decomposition
m\Gamma' .
endfor ';
if bound (2.11) is satisfied then goto
Use the Arnoldi decomposition M \Gamma1 AV (m\Gammak)
k as input to Algorithm 2.1 and apply the Arnoldi process
to compute the Arnoldi decomposition M
.
Apply GMRES(m): determine the matrices Vm+1 and -
Hm defined by
and (2.2), respectively.
Compute solution ym
endfor fi;
1: Improve approximate solution by GMRES(k) and update preconditioner:
The
k , as well as the ma-
Adaptive preconditioners 15
trices V k+1 and -
and (2.2), are available.
Compute solution y k 2 R k of min
2: M \Gamma1
3: r j+k :=
solution then done;
endfor ff;
while kr j k=kr solution do
Apply
M by Algorithm 2.1 with initial vector
the matrices Vm+1 and -
Hm defined by (2.1) and (2.2), respectively.
Compute solution ym
4: r j+m :=
endwhile;The comments regarding the lines with labels "2:", "3:" and "4:" for Algorithm
3.5 also apply to Algorithm 3.6.
4. Numerical experiments. All the numerical experiments presented in this
section were carried out on an HP 9000/735 computer using MATLAB. In all examples
we chose the initial approximate solution x b. The vector
b had randomly generated uniformly distributed entries in the open interval (0; 1).
The purpose of the experiments was to compare Algorithms 3.5 and 3.6 to a restarted
where the parameter m 0 is chosen so that the latter algorithm
is allowed at least as much computer storage as the former two algorithms. We also
compare Algorithms 3.5 and 3.6 to the GMRES algorithm without restarts, and refer
to the latter scheme as "Full GMRES". We terminated the iterations with these
iterative methods as soon as a residual vector r j was determined, such that
with ffl solution . For Algorithms 3.5 and 3.6, we chose the input parameter
values ffl subspace 20. The storage
requirement for both Algorithms 3.5 and 3.6 with this choice of parameters is at most
54 n-vectors. We compare these schemes with the restarted GMRES(60) algorithm,
which requires the storage of 62 n-vectors for V 61 and xm ; see Algorithm 3.1. This
storage count assumes that the residual vector r m in Algorithm 3.1 up to a scaling
factor is stored in the first column of the matrix V 61 .
Example 4.1. Let the matrix A 2 R 200\Theta200 be partitioned according to
A T
1;2 A 2;2
where A 1;1 2 R 30\Theta30 is a circulant matrix with first row [\Gamma3=2; 2]. The entries
of the diagonal matrix A 2;2 2 R 170\Theta170 are uniformly distributed random numbers in
the interval (1; 10). The matrix A 1;2 is a zero matrix of appropriate order. Thus, the
J. Baglama et al.
matrix A has 30 eigenvalues on a circle with center \Gamma3=2 and radius 2. The remaining
eigenvalues are uniformly distributed in the open interval (1; 10). Figure 4.1 shows
denotes the last preconditioner of the
computed by Algorithm 3.5 with shifts (2.10). The eigenvalues are shown
for the unscaled matrix A, and the eigenvalues for M \Gamma1 A are also for the unscaled
matrix A and the associated preconditioner. The unscaled preconditioner maps the
eigenvalues of A of smallest magnitude to approximately sign(Re(- n ))j- n j. This is
illustrated by Figure 4.1. Figure 4.2 shows that the iterates converge rapidly when the
preconditioner has removed many of the eigenvalues on the circle fz 2g.
We remark that the plot of -(M determined by Algorithm 3.6
looks roughly the same as the plot of the eigenvalues of the preconditioner shown in

Figure

4.1.
The graph for Algorithm 3.5 in Figure 4.2 (continuous curve) was generated by
evaluating kr j k for every value of j for which the residual vector r j is defined, i.e.,
after every step of Richardson iteration in, and after every minimization of the residual
error by the GMRES algorithm. The graph for Algorithm 3.6 in Figure 4.2 (dashed
curve) was generated by evaluating kr j k after every minimization of the residual error
by the GMRES algorithm. The number of matrix-vector products with the matrix A
reported in Table 4.1, however, is only the number actually required by Algorithms
3.5 and 3.6. The piecewise linear graph for GMRES(60) in Figure 4.2 is obtained
by linear interpolation between the nodes
The nodes are marked with circles. The column "size of Krylov subspace" in Table
4.1 displays the parameter m used for Algorithms 3.1, 3.5 and 3.6. The column "#
preconditioners" shows the number of preconditioners M \Gamma1
used before a sufficiently
accurate solution was found. This number is bounded by ff 0 . The column "# vectors
in each preconditioner" is the parameter k in Algorithms 3.5 and 3.6. The column
labeled "total # vectors used" counts the number of n-vectors in storage.
The graph in Figure 4.2 (dash-dotted curve) for "Full GMRES" is obtained by
applying GMRES(m) to the solution of (1.1) for increasing values of m in order to
improve the initial approximate solution x 0 until an approximate solution xm with
a sufficiently small residual error kr m k has been determined. Figure 4.2 shows the
10-logarithm of the relative residual error kr k k=kr 0 k for all
Example 4.2. Consider the 200 \Theta 200 block bidiagonal matrix
\Gammay
\Gammay
\Gammay
Its eigenvalues are given by - 2j
\Gamma1.

Figures

4.3 and 4.4 are analogous to Figures 4.1 and 4.2, respectively, and Table
4.2 is analogous to Table 4.1. The distribution of eigenvalues of M \Gamma1 A in Figure
4.3 indicates that the tolerance ffl subspace used in the computations is too
large to determine an accurate approximate invariant subspace of A. Nevertheless,
Adaptive preconditioners 17
the eigenvalues of A closest to the origin were removed, and Algorithms 3.5 and 3.6
yield faster convergence than the restarted GMRES(60) algorithm; see Figure 4.4. 2
Example 4.3. Let is the Pores3 matrix of the Harwell-Boeing
matrix collection. The matrix A 0 is nonsymmetric, of order
3474 non-zero entries. The purpose of the shift was to obtain a matrix with
some positive eigenvalues. Figures 4.5 and 4.6 are analogous to Figures 4.1 and 4.2,
respectively, and Table 4.3 is analogous to Table 4.1. We can see that some eigenvalues
of the matrix A are very close to the origin and others are of large magnitude. Figure
4.5 illustrates how the the preconditioner moves eigenvalues of A away from the origin
to approximately sign(Re(- n ))j- n j, which is negative. Figure 4.6 shows the rate of
convergence. 2
Example 4.4. Let A be a diagonal matrix of order 200 with diagonal entries
a

Figures

4.7 and 4.9 are analogous to Figures 4.1 and 4.2, respectively, and Table 4.4
is analogous to Table 4.1. Figure 4.8 illustrates that the preconditioner moved all the
smallest eigenvalues of A, except for one, away from the origin. Figure 4.9 shows the
rate of convergence. 2
Example 4.5. In all examples above, we chose the shifts according to (2.10),
i.e., we determined approximations of subspaces associated with a few eigenvalues
of smallest magnitude. The present example illustrates that the Algorithms 3.5 and
3.6 easily can be modified to determine approximations of other invariant subspaces.
Specifically, we used Algorithm 3.6 to solve the same linear system of equations as in
Example 4.1, and chose as shifts the eigenvalues with largest real part of the
matrices Hm generated during the iterations. Thus, we sought to determine invariant
subspaces associated with a few of the eigenvalues with smallest real part. Figure
4.10 is analogous to Figure 4.1 and shows -(A) (dots) and -(M \Gamma1
eigenvalues of A on the circle were removed, and the number of matrix-vector products
required before the stopping criterion was satisfied was 311, which is less than the
numbers of matrix-vector products reported in Table 4.1. 2
5. Conclusion. This paper describes new preconditioning methods that are well
suited for use with the restarted GMRES(m) algorithm. Numerous computed examples
indicate that iterates generated by our methods can converge significantly faster
than iterates determined by a restarted GMRES algorithm that requires more computer
storage. Algorithms 3.5 and 3.6 describe versions of our preconditioning method
in which the eigenvalues - j of smallest magnitude of the matrix A are mapped to approximately
that it is easy to modify our
preconditioners so that other eigenvalues are mapped.

Acknowledgements

. Work on this paper was carried out while the last three authors
visited the Computer Science Department and IPS at the ETH. They would like to
thank Walter Gander and Martin Gutknecht for making these visits possible. We
would like to thank Marcus Grote for discussions and code for extracting matrices
from the Harwell-Boeing matrix collection, and Richard Lehoucq for providing us
with reference [16].
J. Baglama et al.



--R

The principle of minimized iterations in the solution of the matrix eigenvalue problem
Cambridge University Press
Iterative methods for computing a few eigenvalues of a large symmetric matrix
A Newton basis GMRES implementation
A parallel implementation of the GMRES al- gorithm
An implicitly restarted Lanczos method for large symmetric eigenvalue problems
Deflated and augmented Krylov subspace techniques
A parallel implementation of the restarted GMRES iterative algorithm for nonsymmetric systems of linear equations
Numerical stability of GMRES
A parallel GMRES version for general sparse matrices
Restarted GMRES preconditioned by deflation
Iterative solution of linear systems
Matrix Computations
GMRES/CR and Arnoldi/Lanczos as matrix approximation problems
Parallel preconditioning with sparse approximate inverses
Eigenvalue translation based preconditioners for the GMRES(k) method
Analysis and implementation of an implicitly restarted Arnoldi iteration
A restarted GMRES method augmented with eigenvectors
A hybrid GMRES algorithm for nonsymmetric linear systems
Krylov subspace methods for solving large unsymmetric systems of linear equations
Preconditioned Krylov subspace methods for CFD applications
GMRES: a generalized minimum residual algorithm for solving non-symmetric linear systems
Implicit application of polynomial filters in a k-step Arnoldi method
Introduction to Numerical Analysis
The superlinear convergence behaviour of GMRES

--TR

--CTR
D. Loghin , D. Ruiz , A. Touhami, Adaptive preconditioners for nonlinear systems of equations, Journal of Computational and Applied Mathematics, v.189 n.1, p.362-374, 1 May 2006
Ronald B. Morgan, Restarted block-GMRES with deflation of eigenvalues, Applied Numerical Mathematics, v.54 n.2, p.222-236, July 2005
Paul J. Harris , Ke Chen, On efficient preconditioners for iterative solution of a Galerkin boundary element equation for the three-dimensional exterior Helmholtz problem, Journal of Computational and Applied Mathematics, v.156 n.2, p.303-318, 15 July
arnoldi process;iterative method;preconditioner;nonsymmetric linear system