Using Generalized Cayley Transformations within an Inexact Rational Krylov Sequence Method.
The rational Krylov sequence (RKS) method is a generalization of Arnoldi's method. It constructs an orthogonal reduction of a matrix pencil into an upper Hessenberg pencil. The RKS method is useful when the matrix pencil may be efficiently factored. This paper considers approximately solving the resulting linear systems with iterative methods. We show that a Cayley transformation leads to a more efficient and robust eigensolver than the usual shift-invert transformation when the linear systems are solved inexactly within the RKS method. A relationship with the recently introduced Jacobi--Davidson method is also established.
Introduction
. Suppose that a few eigenvalues near a complex number - and
possibly corresponding eigenvectors of the generalized matrix eigenvalue problem
are needed. Assume that both A and B are large complex matrices of order n:
Also suppose that at least one of A or B is nonsingular so that equation (1.1) has
eigenvalues. Without loss of generality, assume that B is invertible. Following
standard convention, we refer to (A; B) as a matrix pencil. For us, n is considered
large when it is prohibitive to compute all the eigenvalues as a dense algorithm in
would attempt to do.
A standard approach is to perform inverse iteration [17, p.386] with the matrix
-B: The sequence of iterates
is produced. Under some mild assumptions, the sequence converges toward the desired
eigenvector with eigenvalue closest to -, and a Rayleigh quotient calculation gives an
estimate of the eigenvalue. Another approach is to extract the approximate eigenpair
by using the information from the subspace defined by joining together m iterates
of the sequence (1.2). This leads to a straightforward extension [22] of the ideas
introduced by Ericsson and Ruhe [13] for the spectral (shift-invert) transformation
Lanczos method. Starting with the vector v, Arnoldi's method [2] builds, step by
step, an orthogonal basis for the Krylov subspace
Km
One improvement to the inverse iteration scheme given is to possibly vary the
at every step. For example, - j may be set to the Rayleigh quotient
z H Az=z H Bz, where z is an unit vector in the direction of (T SI
elegantly shows how to build an orthogonal basis for the rational Krylov subspace
The work of R. B. Lehoucq was supported by the Mathematical, Information, and Computational
Sciences Division subprogram of the Office of Computational and Technology Research, U.S.
Department of Energy, under Contract W-31-109-Eng-38. The work by Karl Meerbergen was supported
by the project Iterative Methods in Scientific Computing , contract number HCM network
CHRC-CT93-0420, coordinated by CERFACS, Toulouse, France.
y Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL 60439
USA. Current address Sandia National Laboratories, MS 1110, P.O.Box 5800 Albuquerque, NM
87185-1110, rlehoucq@cs.sandia.gov.
z LMS Numerical Technologies, Interleuvenlaan 70, 3001 Heverlee, Belgium, km@lmsnit.be
R. B. LEHOUCQ AND K. MEERBERGEN
The resulting algorithm is called a rational Krylov sequence (RKS) method and is
a generalization of the shift-invert Arnoldi method where the shift is possibly varied
during each step.
All the methods considered require the solution of By for x: This
is typically accomplished by factoring A \Gamma -B: For example, when A \Gamma -B is sparse,
a direct method [7, 8, 9, 10, 12, 11] may be employed. If the shifts - j are not varied,
then use of one of these direct methods in conjunction with ARPACK [21] is a powerful
combination for computing a few solutions of the generalized eigenvalue problem (1.1).
However, for large eigenvalue problems (n ? 10; 000), direct methods using the
RKS method may not provide an efficient solution because of the potentially prohibitive
storage requirements. The motivation for the current study is to investigate
the use of iterative methods for the linear systems of equations arising in the RKS
method. One benefit is that for the many eigenvalue problems arising from a discretization
of partial differential equations, an intelligent preconditioner may often
be constructed. We shall call these methods inexact RKS ones because we no longer
have a rational Krylov space. In particular, we shall demonstrate that a Cayley transformation
performs more robustly than a shift-invert
transformation T SI
using iterative methods for the linear
solves.
Before we continue, some remarks are in order. Although, combining an eigensolver
(using one of the methods discussed previously) with an iterative method for
the linear solves is not a new (or even novel) idea, what is generally not appreciated
is that residuals of the linear systems must be small. To be precise, the matrix vector
product
must be applied so that kBu \Gamma
where v is the approximate solution of the linear system and ffl M is machine precision.
This is a necessary requirement for the correct representation of the underlying Krylov
subspace. If the linear systems are not solved with the above accuracy, there is no
guarantee that a Krylov space for T SI
j has been generated. For example, if Arnoldi's
method is used, there is no reason to expect that the Hessenberg matrix generated
represents the orthogonal projection of T SI
onto the Arnoldi vectors generated. If the
above assumption of accuracy is violated (as is often the case), any results produced
by such an eigensolver should be taken with caution.
Fittingly, the literature on approaches for finding a few solutions to the generalized
eigenvalue problem (1.1), where only approximate solutions to the linear systems
are available, is sparse. Bramble et. al. [4], Knyazev [18], Knyazev et. al. [19], Morgan
[25], and Szyld [42] each consider the situation where the matrix pencil is symmetric
positive definite. (The papers [18, 19, 4] also contain numerous citations to the
Russian literature.) Algorithms based on Jacobi-Davidson method [38] introduced by
Sleijpen and van der Vorst are discussed in [14, 39]. In a recent report, Sorensen [41]
discusses methods based on truncating a QZ iteration. The recent paper by Meerbergen
and Roose [23] provided motivation for the current article. They demonstrate the
superior numerical performance of a Cayley transformation over that of a shift-invert
transformation within an Arnoldi method when using an iterative linear solver.
Our article is organized as follows. We introduce the RKS method in x2. The
inexact RKS method is introduced in x3 along with a connection with inverse iteration
and some examples illustrating our ideas are presented. In x4, we illustrate our
method for a generalized eigenvalue problem. In x5, we show that an appropriate
(approximate) shift-invert transformation could be used. We compare inexact RKS
and Jacobi-Davidson methods in x6. We conclude the paper in x7 with a summary
INEXACT RATIONAL KRYLOV SEQUENCE METHOD 3
ffl Choose a starting vector v1 with
ffl For
1. Select a pole - j , a zero - j
2. (continuation vector).
3. Form w /
4. Orthogonalize w
5.
\Theta
and ~ k
7.
and ~
8. Compute approximate eigenpairs of interest.
9. Check whether the approximate eigenpairs satisfy the convergence criterion.
Fig. 2.1. Computing the Rational Krylov Sequence (RKS) for the matrix pencil (A,B)
of the main ideas and some remaining questions.
In this article, matrices are denoted by upper-case Roman characters. Vectors are
denoted by lower-case Roman characters. The range of the matrix V is denoted by
R(V The Hermitian transpose of the vector x is denoted by x
will be introduced and employed in the next few sections. The norm k \Delta k used is the
Euclidean one.
2. The Rational Krylov Sequence Method. The method is outlined by
the algorithm listed in Figure 2.1. For the practical RKS algorithm given in [32],
Ruhe considers the shift-invert transformation T SI
In exact arithmetic, both transformations lead to the
same rational Krylov space, because
However, in finite-precision arithmetic and/or in conjunction with iterative methods
for linear systems, substantial differences may exist (see [23] for examples). We call
the - j 's the poles, the - j 's the zeros, and the 's the continuation vectors. A
discussion on some possible choices is postponed until x3.2. This section will discuss
some relationships among quantities in Steps 1-7, the form of Gram-Schmidt orthogonalization
we employ, and finally the computation of approximate eigenpairs and
their convergence.
By eliminating w from Steps 2-5, we obtain the relationship
\Theta h 1;j h 2;j
\Theta t T
. Rearranging Equation
(2.2) results in
By putting together the relations for we have that
4 R. B. LEHOUCQ AND K. MEERBERGEN
where ~ h j and ~ t j are associated with the jth columns of ~
Hm and ~
Tm , respectively, and
A final simplification is to rewrite Equation (2.3) as
where
~
Tm and ~
We remark that as long as the sub-diagonal elements (the h j+1;j 's) are nonzero, both
~
Hm and ~
Lm are unreduced upper Hessenberg (rectangular) matrices and hence of full
rank.
2.1. Orthogonalization. The orthogonalization of Step 3 of the algorithm in

Figure

2.1 is performed using an iterative classical Gram-Schmidt algorithm. This is
the same approach used by Sorensen [40] based on the analysis [5] of reorthogonaliza-
tion in the Gram-Schmidt algorithm.
2.2. Computing Eigenvalue Estimates. We now consider the calculation of
approximate eigenpairs for the RKS method and first discuss how to compute Ritz
pairs. The main purpose of this article is to study the use of iterative linear system
solvers in RKS and not the various ways to extract eigenvalues. Therefore, we
use standard Ritz values throughout, though the theory can easily be extended to
harmonic [32, 39] Ritz values.
Consider a matrix C and a subspace R(X), where X 2 C n\Thetak is of full rank. The
pair ('; y jXz) is called a Ritz pair of C with respect to the subspace R(X) if and
only if
This is referred to as a Galerkin projection. Two important properties of a Galerkin
projection are the following. First, if R(X) j C n , the Ritz pairs are exact eigenpairs
of C. Second, if C is normal, the Ritz values lie in the convex hull of the eigenvalues
of C. For example, if C is Hermitian, the Ritz values lie between the smallest and
largest eigenvalue of C.
The following theorem shows how Ritz pairs may be computed from the RKS
method outlined by the algorithm listed in Figure 2.1.
Theorem 2.1. ('; y j Vm+1 ~
Lm z) is a Ritz pair for B \Gamma1 A with respect to
only if
~
Lm z:
Proof. Following the definition (2.6) and Equation (2.4), ('; y) is a Ritz pair when
Lm
Thus, (Vm+1 ~
and the desired equivalence with (2.7)
follows.
We denote by ' (m)
i the i-th Ritz value available after m steps of the RKS algorithm
of

Figure

2.1. Unless otherwise stated, we assume that the Ritz values are in increasing
distance from -m , that is, j' (m)
j: The associated
Ritz vector is denoted by y (m)
The sub- and superscripts are omitted whenever their
meaning from the context is clear.
INEXACT RATIONAL KRYLOV SEQUENCE METHOD 5
2.2.1. Computing Ritz Pairs. The generalized eigenvalue problem (2.7) may
be solved as a standard one. Since ~
Lm is an unreduced upper Hessenberg matrix, ~
Lm
is of full rank, and hence ~
Lm is invertible. Thus, the standard eigenvalue problem
~
is solved giving the 1 -
We remark that ~
m is the
Moore-Penrose generalized inverse of ~
: The explicit formation of the inverse of
~
Lm is not required. Instead, ~
Km may be computed by least squares methods,
for example with the LAPACK [1] software. The Ritz vector is y (m)
Lm z (m)
where
1: Sub- and superscripts are omitted when the context is clear.
2.3. Stopping Criterion. The accuracy of a Ritz pair (';
typically estimated by the residual norm kAy \Gamma By'k: From Equation (2.4), it follows
that
Lm )z (m)
where g (m)
Km z (m)
Lm z (m)
simple check for convergence of a Ritz pair
in the algorithm in Figure 2.1 is when
for a user-defined error tolerance tol.
For any Ritz pair ('; y) it follows that (A+E)y = By', where Hence
small relative to kB \Gamma1 Ak, then ('; y) is an
eigenpair for a nearby problem. If ' is not a poorly conditioned eigenvalue of the
matrix pencil and kB \Gamma1 k is not large, then the size of kgk indicates the accuracy of
the computed Ritz value.
This conclusion motivates us to say that the sequence of Ritz pairs (' (m)
(fixed i) converges toward an eigenpair of Equation (1.1) if and only if kg (m)
to zero as m increases toward n: Although this convergence is not rigorously defined
(we necessarily have kg (n)
does allow us to track the progress of a Ritz pair
after step m of algorithm in Figure 2.1.
3. The Inexact RKS Method. At Steps 3-5 of the RKS algorithm in Figure
2.1 the Cayley transformation
is computed by a two step process. First, the linear system
is solved for w: Next, w is orthogonalized against V j , and the solution V j+1 ~ h j results.
These two steps account for the largest source of errors arising when computing in
floating-point arithmetic. Since our interest is in using a (preconditioned) iterative
method for the solution of Equation (3.1), we neglect the errors in the Gram-Schmidt
orthogonalization phase (but we assume that the columns of V j+1 are orthogonal to
machine precision).
6 R. B. LEHOUCQ AND K. MEERBERGEN
Let us formally analyze the errors arising from the solution of Equation (3.1). Let
denote the computed solution and the
associated residual. Thus,
Here, ks j x H
ks ks is a modest multiple of
machine precision, we say that the direct method computes a backward stable solution.
A robust implementation of a direct method gives a backward stable solution to a
linear system. Note that even if a backward stable solution x j is in hand, it may
share few, if any, digits of accuracy with w: Moreover, achieving such a backward
stable solution with an iterative method may be prohibitively expensive. Therefore,
we shall study the situation where a large backward error is allowed for the solution
of the linear system.
In order to give an indication of what we mean by large, a few words about
iterative linear system solvers are needed. A linear system is said to be solved
with a relative residual tolerance - when the solution, x, satisfies kb \Gamma Cxk -kbk for
any b: Krylov methods [15, 34] are typically used. GMRES [35], BiCGSTAB(') [37],
and QMR [16] are among those most widely used. See [3] for templates for all these
solvers. The performance of these solvers substantially improves when a suitable
preconditioner is employed. Hence what we mean by a large error is that 10 \Gamma8 -
By putting all the s j for together in Sm j
; we have
which we call an inexact rational Krylov sequence (I-RKS) relation. This relation
may be rewritten as
Km with
where ~
m is the generalized Moore-Penrose inverse. In other words,
we have computed an exact RKS for the pencil We caution the reader
not to confuse the Em 's with the unsubscripted E's of x 2.3.
Denote by oe
Lm ) the reciprocal of the minimum singular value of ~
fore, if
is large, then the Ritz pairs from x2.2.1 may not be those of a pencil near (A; B): This
situation implies that even if we use a direct method for the linear systems, a nearly
rank deficient ~
Lm might lead to inaccurate Ritz pairs. The matrix Em incorporates
the backward error of the linear solution and is the distance to the matrix pencil
We call the Ritz pairs for
We now define and discuss a few quantities that will prove helpful in the discussion
that follows.
ffl Cayley residual s C
this is the residual of the linear system (3.1) at step j of
the rational Krylov method.
ffl RKS residual f (j)
the RKS method computes a Ritz pair (' (j)
i and
and so the RKS residual
satisfies
~
INEXACT RATIONAL KRYLOV SEQUENCE METHOD 7
ffl Choose a starting vector v1 with
ffl For
1. Select a pole - j . If j ? 1 set the zero -
Otherwise, set
2. Compute the continuation vector r (j
3. Form by solving
ae
oe
for x j and set
4. See steps 4-6 of the IC-RKS method listed in Figure 3.2.
Fig. 3.1. Inverse iteration via the inexact RKS method.
ffl True residual r (j)
this is the residual defined by r (j)
(The sub- and superscript of the true residual is dropped with those of f (j)
These three residuals may be linked via the relationships
that follow from Equation (3.2) and the definition (3.3) of present numerical
evidence that demonstrates that although kf (j) k decreases in size for increasing j; r (j)
does not decrease when an inexact shift-invert transformation is employed. However,
when an inexact Cayley transformation is used instead, both kS j z (j) k and kf (j) k
decrease and the size of the true residual also decreases.
The continuation of this section is as follows. In x3.1, we present a relationship
with inverse iteration that includes a theorem that shows the convergence for inexact
inverse iteration. In x3.2, we fix the various parameters of the RKS method,
i.e. the poles, zeros and continuation vectors. This selection makes a link with the
generalized Davidson method [6, 25, 26]. In x3.3, an informal argument is given for
the convergence of the inexact Cayley rational Krylov sequence (IC-RKS) method,
described under x3.2, using the theoretical result from x3.1. We also illustrate this by
a numerical example.
3.1. Inverse Iteration. We first exploit a direct relationship with inverse iteration
that occurs with a special choice of the continuation vector when a Cayley
transformation is used. An example is then presented that compares this choice with
a shift-invert transformation. The subsection is concluded with a theorem that shows
the numerical behavior observed is not just a fortuitous event. Although the choice
of continuation vector does not exploit the entire space of vectors as in IC-RKS, the
theorem justifies the superior properties of combining an approximate linear solve via
a Cayley transformation.
From Equation (2.2) and the matrix identity (2.1), it follows that
Using (2.5) with ~ l
8 R. B. LEHOUCQ AND K. MEERBERGEN
and hence V j+1 ~ l j is the linear combination of the columns of V j+1 obtained by performing
one step of inverse iteration on the vector V j t j . An inductive argument easily
establishes the following property.
Lemma 3.1. If t
Y
is a scalar and v 1 is the starting vector of RKS.
Lemma 3.1 indicates how to compute an approximate eigenvalue. If we denote
~
Equation (2.4) gives the Rayleigh quotient
as an estimate of an eigenvalue without need to explicitly apply B
An algorithm for inverse iteration is given in Figure 3.1. The approximate eigenpair
on iteration j is (' so we can use the relationships (3.4)
with z Recall that we used - (j \Gamma1) and . The entries
' (0) and v 1 determine the initial estimates for the eigenpair. We now compare inexact
inverse iteration computed via the RKS method using the shift-invert and Cayley
transformations with an example.
Example 3.1. The Olmstead model [28] represents the flow of a layer of viscoelastic
fluid heated from below. The equations are
@t
with boundary conditions
the speed of the fluid and v is related to viscoelastic forces. The equation was discretized
with central differences with gridsize 1=(n=2). After the discretization,
the equation may be written as -
The size of the Jacobian matrix We consider the Jacobian
for the parameter values for the trivial steady state
Thus, the interest is in the eigenvalue of largest real part.
We ran the algorithm in Figure 3.1. The linear systems were solved by 20 iterations
of Gauss-Seidel starting with a zero initial vector. Since this solver is stationary,
the relative residual norm is almost constant. The initial guess for the eigenvalue was
0: The initial vector for RKS was
n: The poles - j were set
equal to 5 for all j: The residuals r (j) , f (j) and S j z (j) are shown in Table 3.1. All
three sequences decrease when the Cayley transform is used.
We redid the experiments using the shift-invert transformation. The results are
also shown in Table 3.1. Both kS j z (j) k and kr (j) k stagnate near the same value. Note,
however, that kf (j) k tends to zero.

Table

3.1 shows that the true residual decreases when the Cayley transformation is
used, but stagnates for the shift-invert transformation. The following result indicates
what occurs under some mild conditions when performing inexact inverse iteration
with either the shift-invert or the Cayley transformation.
INEXACT RATIONAL KRYLOV SEQUENCE METHOD 9

Table
Numerical results for inverse iteration on Example 3.1 using inexact Cayley and shift-invert
transformations. The table shows the norms of true residual r (j) , S j z (j) , and the RKS residual
The norm of ~ l j is also displayed for the Cayley transformation.
Cayley shift-invert
Theorem 3.2. Assume that there is an integer k - m and value fl ? 0 such that
is the relative residual tolerance used for the linear solves (see equations (3.12) and
(3.13)).
If a Cayley transformation is used, then for
kf
and when a shift-invert transformation is used,
Proof. With z
it follows that
ks j k=k ~ l j k
ks
For the Cayley transform, we prove (3.7) by induction on j. We clearly have that
which satisfies (3.7) for k. Suppose that (3.7) holds for some integer
From the hypothesis of the theorem, we have that
Combining this with equations (3.9) and (3.12) results in
ks (j \Gamma1) k:
Using our inductive hypothesis on kr (j \Gamma1) k gives
and (3.7) follows. For shift-invert, (3.8) follows from (3.9) and (3.13), which completes
the proof.
The theorem shows that if ae iteration computed via
the Cayley transformation will produce a Ritz pair with a small direct residual. Since
inexact inverse iteration can do no better than exact inverse iteration.
Although, the term kf (j) k will decrease when using the shift-invert transformation, the
size of the direct residual kr (j) k may stagnate. This occurs because the contribution
from solving the linear systems inexactly (s SI
j ) to the true residual is constant. When
a direct method is used for the linear system of equations, - is a multiple of machine
precision. Hence, whether a shift-invert or Cayley transformation is used, the true
residual kr (j) k decreases at a rate proportional to ae:
For the exact Cayley transformation, we have
and ~ l
1. Hence, we have
Thus, converges to zero, then 1 \Sigma
tends to one for increasing j: Computation reveals that quite often k ~ l j k - 1 after a
very small number of steps. This also holds for inexact inverse iteration, because it can
be seen as exact inverse iteration applied to B), as Table 3.1 demonstrates.
Hence, for large enough k, and the convergence rate of inverse iteration using
the Cayley transform is approximately ae As the method progresses, ae is easily
estimated and thus the largest relative residual tolerance that may be used is also
easily estimated.
3.2. Choosing a Pole, Zero and Continuation Vector. A robust and efficient
strategy for selecting the poles during the RKS method is a subject of research.
The present situation is further complicated because we employ approximate methods
for the linear solves. Since we are more concerned in showing that use of ~
K k and ~
for the computation of Ritz pairs, a fixed pole is used for the numerical experiments.
The choice of the zero of the Cayley transformation is crucial for computing a Ritz
pair with a small direct residual. This was demonstrated by the numerical examples
in [23]. We first formally analyze the choice of the zero and continuation vector and
then give an example.
Suppose that (' (j \Gamma1) ; y (j \Gamma1) ) is an (inexact) Ritz pair computed during the j-1st
step of an (inexact) RKS method. We select the zero - (j \Gamma1) and contination
vector
of interest.
For a Cayley transformation, this leads to
while a shift-invert transformation gives
(3.
INEXACT RATIONAL KRYLOV SEQUENCE METHOD 11
ffl Choose a starting vector v 1 with kv 1
ffl For
1. Select a pole - j . If j ? 1 set the zero -
(j \Gamma1) .
Otherwise, set -
2. Compute the continuation vector r (j
residual).
3. Form w /
4. See steps 4-7 of the RKS method listed in Figure 2.1.
5. Solve the eigenvalue problem ~
6. Check whether the approximate eigenpairs satisfy the convergence criterion

Fig. 3.2. Computing eigenvalues of the pencil (A; B) by the inexact Cayley rational Krylov
sequence (IC-RKS) method
as the linear systems to be solved. Although both transformations use the same continuation
vector, the Cayley transformation also uses the Ritz value for its zero. The
only difference in the two linear systems (3.10) and (3.11) is the righthand side. When
a preconditioner is used to solve the linear system (3.10), we have a generalization of
Davidson's method [6, 26] for computing eigenvalues of a matrix pencil.
Denote the computed solutions to (3.10) and (3.11) by x C
j and x SI
If an iterative method with relative residual tolerance - is used for the two linear
systems, then the residuals of the linear systems satisfy
ks C
ks SI
for the Cayley and shift-invert transformation, respectively. (We drop the superscripts
that denote whether a Cayley or shift-invert transformation when the context is clear.)
In view of the two bounds (3.12) and (3.13) on the computed solutions, a Cayley
transformation is preferred over a shift-invert transformation. It appears that use of
a Cayley transformation leads to better results with inexact linear solvers when the
zero and continuation vector are chosen as in (3.10). Our experimental results also
support this conclusion. The algorithm in Figure 3.2 lists an inexact Cayley RKS
method (IC-RKS). We now illustrate a few properties of this algorithm by means
of an example that demonstrates: (1) the inexact rational Krylov method is not a
Galerkin projection method; (2) the method can only compute one eigenvalue at a
time, just as in Davidson methods.
Example 3.2. Consider the matrices
(A; I) has eigenpairs (j; e j 5: The goal is to compute the smallest eigenvalue
1 and corresponding eigenvector e 1 with IC-RKS method using a fixed pole
0:7. The starting vector is set equal to v
5: The Cayley system
12 R. B. LEHOUCQ AND K. MEERBERGEN
is solved as x
Note that M simulates a stationary iterative solver with residual tolerance
which implies that f (5)
5: Thus, the computed eigenpairs are exact
eigenpairs of A+E 5 : We found that
1:0000 0:0120 \Gamma0:0697 0:3708 \Gamma0:4728
\Gamma0:0000 1:9987 \Gamma0:5981 4:4591 \Gamma5:6013
\Gamma0:0001 0:1003 0:4666 17:1897 \Gamma21:1757
\Gamma0:0002 0:0340 \Gamma4:4220 36:8172 \Gamma40:7251
and has eigenpairs
1:0000
0:0000
0:0000
0:0000
0:0000
0:9812
\Gamma0:1801
\Gamma0:0177
\Gamma0:1340
0:9229
0:3188
0:1681
0:0045 \Sigma 0:0068i
\Gamma0:0826 \Sigma 0:7214i
the true residual has the form r (5)
This example shows that E 5 is nearly rank deficient and that the desired eigenvector
of (A; I) is nearly its nullvector. Therefore, the desired eigenvalue, in this case,
can be computed with a small true residual. It should be noted that the
perturbation E 5 is small in the direction of only one eigenspace, hence IC-RKS is not
able to compute several eigenvalues simultaneously. This is not the situation when
the linear systems are solved more accurately with, for instance, a direct method.
In this example, IC-RKS computes the exact eigenpairs of A
steps. In general, however, r (5)
because the inexact Ritz pair is not computed
from a Galerkin projection with A: We also remark that ' (5)
4 and ' (5)
5 are non-real and
this would not be the case with a Galerkin projection because A is a real symmetric
matrix. This is in contrast with other iterative eigenvalue solvers, such as Arnoldi
and Jacobi-Davidson methods where Galerkin projections with A are employed.
3.3. Inexact Rational Krylov. We now informally discuss the algorithm listed
in

Figure

3.2 including a comparison with inexact inverse iteration of the previous
section.
From (3.5) with the Ritz vector y (i\Gamma1) computed as in x 2.2.1, it follows that
INEXACT RATIONAL KRYLOV SEQUENCE METHOD 13

Table
Numerical results for the Olmstead model of Example 3.3. The table shows the order of accuracy
for the residual norm of the rightmost Ritz pair, the norm of S j z (j) , and the first four components
of z (j) .
Numerical experiments reveal that the j-th component of z (j) is large relative the
components (see Table 3.2). This is because the best approximation of
the desired eigenvector among the columns of V j+1 ~
L j is given by V j+1 ~ l j -the improvement
of the previous Ritz vector by one step of inverse iteration. Thus, using
the continuation vector V i ~
should give better results because information
in the subspace R(V i+1 ~
used. Inexact inverse iteration only uses information in
the space spanned by the last column of V i+1 ~
The inexact Ritz pairs (' (i) ; y (i) ) lead to true residuals r (i) if the Cayley transform
is used. The Cayley residual on iteration i satisfies ks i k -kr (i\Gamma1) k. The true residual
on the jth iteration is decomposed as r
ks
gives an upper bound to kS j z (j) k. In the right-hand side, ks i k is independent of j and
can be quite large for small i. However, because je T
typically forms a decreasing
sequence for increasing j, we have a decreasing sequence kS j z (j) k:
Example 3.3. We now discuss an example for which e T
z (j) and S j z (j) tend
to zero in the IC-RKS method. The matrix arises from the same problem as in
Example 3.1, but now We ran Algorithm IC-RKS from Figure 3.2 with
fixed starting with vector
n: The linear systems were
solved by GMRES preconditioned by ILU. The number of iterations of GMRES was
determined by the relative error tolerance, which was selected as
shows the residual norm and the norm of the error term S j z (j) . Both kS j z (j) k and
tend to zero. For large j, kS j z (j) k - kr (j) k. This is the case because f (j)
converges more rapidly to zero than S j z (j) . Table 3.2 also illustrates the fact that
decreases for a fixed i and increasing j.
4. A Numerical Example. This example illustrates the use of inexact rational
Krylov methods for the solution of a generalized eigenvalue problem. We also make a
comparison between inexact inverse iteration with the Cayley transform and IC-RKS.
The simulation of flow of a viscous fluid with a free surface on a tilted plane,
leads, with a finite element approach, to an eigenvalue problem
singular matrix. The computation of the
eigenvalue nearest \Gamma10 is of interest. Since our theory is valid only for nonsingular
14 R. B. LEHOUCQ AND K. MEERBERGEN

Table
Numerical results for the tilted plane problem from x 4. The methods used are inexact rational
Krylov (IC-RKS) and inverse iteration with the Cayley transform. On iteration j, ' (j) is the inexact
Ritz value, s j the Cayley residual, and g
IC-RKS (Fig. 3.2) Inverse Iteration (Fig. 3.1)
ks
ks
9
B, we interchange the role of A and B by computing the eigenvalue
nearest
The fact that B is singular implies that is an eigenvalue. It has been shown
that the presence of this eigenvalue can disturb the calculation of a nonzero eigenvalue
when the shift-invert Arnoldi method [29, 24], or the rational Krylov method [36] are
used. One way to reduce the impact of is to start the IC-RKS method with an
initial vector v 1 that is poor in the eigenspace corresponding to This can
be achieved by selecting arbitrary.
The eigenvalue fl nearest \Gamma0:1 was computed by use of IC-RKS (Fig. 3.2) with
fixed pole - \Gamma0:1. The linear systems were solved by GMRES preconditioned with
. The initial vector v 1 was computed
from the system using the GMRES-ILUT
solver. The algorithm was stopped when kr (j)
The numerical results are shown in Table 4.1 for inexact rational Krylov (IC-
RKS) and inexact inverse iteration using the Cayley transform. First, note that
kf (j) k - kAkkg (j) k, so kg (j) k does not measure the RKS residual (see also (2.8)).
Also note that for both IC-RKS and inverse iteration, the sequences kr (j) k, ks j k and
decrease. Both methods converge to Finally, note that
IC-RKS is faster than inverse iteration.
5. A relation between inexact shift-invert and Cayley transforms. In
the previous section, we showed that the inexact rational Krylov method can be used
for the computation of eigenvalues of a matrix pencil. The example shows a substantial
difference in convergence behavior between the shift-invert and Cayley transforma-
tions. In this section, we show that an appropriate shift-invert transformation may
also be employed.
During each step of IC-RKS, the following relationship
results where s j is the residual of the linear system that is approximately solved.
Rearranging (5.1) and adding - j By (j) to both sides gives the equivalent shift-invert
system
INEXACT RATIONAL KRYLOV SEQUENCE METHOD 15
Hence, if the zero vector is used as the initial guess for the iterative method for linear
systems approximately solved via the Cayley transform, \Gammay (j) should be used for the
shift-invert transformation formulation.
Assume that - is a constant and that IC-RKS converges to some eigenpair. From
(3.13), it follows that when shift-invert is used, convergence to the same eigenpair is
attained for decreasing - (as j increases). In the context of inexact inverse iteration,
Lai, Lin and Lin [20] also observe that the approximate linear system solver requires
an increasingly tighter tolerance on the residual (of the linear system) as the number
of inverse iterations increases. In contrast, a Cayley transformation allows us to use
a fixed tolerance on the linear system residual.
6. A Connection with the Jacobi-Davidson Method. We now show a
connection between the Jacobi-Davidson [14, 39, 38] and RKS [32] methods.
Consider the linear system
(j \Gamma1) is a Ritz vector of interest. This amounts to selecting the
jth continuation vector t
(j \Gamma1) as in the Algorithm IC-RKS in Figure 3.2
with associated Ritz value
~
The right-hand side in (6.1) is then the residual of the eigenpair (~- j ; y j ) and is orthogonal
to y j . Since we are interested in expanding our search space (the span
of the columns of V j ), multiply both sides of Equation (6.1) by the projector I \Gamma
By
Using the fact that
results in
By
the component of w in the direction of y j does not play a role
when w is added to the subspace R(V j ). Thus, we are interested in finding only the
component of w orthogonal to y j and so the linear system
By
is solved instead. The Jacobi-Davidson method calls Equation (6.2) the correction
equation. Suppose that x j is a computed solution of Equation (6.2) with residual s j ,
given by
By
where s j is orthogonal to y j . Rewrite (6.3) with d j
The orthogonality of y j with d and with s j leads to
(6.
R. B. LEHOUCQ AND K. MEERBERGEN
Choosing the zero - between the Jacobi-Davidson
and RKS methods when Cayley transformations are used. When " j is computed, the
solution of the Jacobi-Davidson correction equation x can be inserted in the
RKS method. Note that, although, the Ritz vector y j is orthogonal to the right-hand
side of the Jacobi-Davidson correction equation (6.2), y j is not orthogonal to the
right-hand side of (6.4).
An advantage of the inexact rational Krylov method is that the matrices ~
and ~
K j do not require the explicit application of A and/or B as needed as in the
Jacobi-Davidson method. An efficient implementation of the Jacobi-Davidson method
requires dot products (the first elements in the last row of V H
We caution the reader to conclude that the Jacobi-Davidson method is an expensive
variant of IC-RKS because it fits an IC-RKS framework. A detailed numerical
comparison of the two methods requires examining the respective rates of convergence
and ability to obtain relative residual reductions during the linear solves. This is the
subject of future work.
7. Conclusions. This paper studied the use of approximate linear solves within
Ruhe's rational Krylov sequence method. The analysis of the convergence of inexact
inverse iteration showed the importance of using the Cayley transformation instead of
the usual shift-invert transformation, when the linear systems are solved with a given
relative residual tolerance.
A theoretical link between the inexact rational Krylov method that uses generalized
Cayley transformations and the Jacobi-Davidson methods was drawn resulting
in a connection between the correction equation and a Cayley transformation.
We called the eigenpairs computed by IC-RKS inexact Ritz pairs, because they
are Ritz pairs for a perturbed RKS method. The classical properties of Galerkin
projection are lost due to this inexactness. Since IC-RKS solves a perturbed problem,
the application of techniques developed for the RKS method (using approximate linear
solves) may be employed. These techniques include the use of complex poles and zeros
for real A and B [31], harmonic Ritz pairs, deflation and purging [32, 36], and the
implicit application of a rational filter [36].

Acknowledgments

The authors thank Dirk Roose for the financial support that
allowed the first author to visit second author. This visit initiated the collaboration
that lead to this article. The authors also thank Gorrik De Sambanx, Gerard Sleijpen
and the referees for helpful comments and suggestions that improved the quality of
the article. In particular, one of the referees provided numerous contructive criticisms
that improved the quality of the presentation.



--R


The principle of minimized iterations in the solution of the matrix eigenvalue problem.
Templates for the solution of linear systems: Building blocks for iterative methods.
A subspace preconditioning algorithm for eigenvector/eigenvalue computation.
Reorthogonalization and stable
The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real symmetric matrices.
An unsymmetric-pattern multifrontal method for sparse LU factor- ization
A combined unifrontal/multifrontal method for unsymmetric sparse matrices.
A supernodal approach to sparse partial pivoting.
ME28: A sparse unsymmetric linear equation solver for complex equations.
The design of MA48
The design of a new frontal code for solving sparse unsymmetric systems.
The spectral transformation Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems.

Iterative solution of linear systems.
QMRPACK: A package of QMR algorithms.
Matrix computations.
Convergence rate estimates for iterative methods for mesh symmetric eigenvalue problem.
The preconditioned gradient-type iterative methods in a subspace for partial generalized symmetric eigenvalue problem
An inexact inverse iteration for large sparse eigenvalue problems.
ARPACK Users' Guide: Solution of Large Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods.
Matrix transformations for computing rightmost eigenvalues of real nonsymmetric matrices.
The restarted arnoldi method applied to iterative linear system solvers for the computation of rightmost eigenvalues.
Implicitly restarted Arnoldi with purification for the shift-invert transformation
Davidson's method and preconditioning for generalized eigenvalue problems.
Generalizations of Davidson's method for computing eigenvalues of large non-symmetric matrices
How to implement the spectral transformation.
Bifurcation with memory.
Improving the spectral transformation block Arnoldi method.
Rational Krylov sequence methods for eigenvalue computation.
The Rational Krylov algorithm for nonsymmetric eigenvalue problems
Rational Krylov
tool kit for sparse matrix computations.
Iterative Methods for Sparse Linear Systems.
GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems.
The implicit application of a rational filter in the rks method.

A Jacobi-Davidson iteration method for linear eigenvalue problems

Implicit application of polynomial filters in a k-step Arnoldi method
Truncated QZ methods for large scale generalized eigenvalue problems.
Criteria for combining inverse and Rayleigh quotient iteration.
--TR

--CTR
Tsung-Min Hwang , Wen-Wei Lin , Wei-Cheng Wang , Weichung Wang, Numerical simulation of three dimensional pyramid quantum dot, Journal of Computational Physics, v.196 n.1, p.208-232, 1 May 2004
arnoldi method;rational Krylov sequence;eigenvalues;cayley transformation;generalized eigenvalue problem