Jacobi--Davidson Style QR and QZ Algorithms for the Reduction of Matrix Pencils.
Recently the Jacobi--Davidson subspace iteration method has been introduced as a new powerful technique for solving a variety of eigenproblems. In this paper we will further exploit this method and enhance it with several techniques so that practical and accurate algorithms are obtained. We will present two algorithms, JDQZ for the generalized eigenproblem and JDQR for the standard eigenproblem, that are based on the iterative construction of a (generalized) partial Schur form. The algorithms are suitable for the efficient computation of several (even multiple) eigenvalues and the corresponding eigenvectors near a user-specified target value in the complex plane. An attractive property of our algorithms is that explicit inversion of operators is avoided, which makes them potentially attractive for very large sparse matrix problems.We will show how effective restarts can be incorporated in the Jacobi--Davidson methods, very similar to the implicit restart procedure for the Arnoldi process. Then we will discuss the use of preconditioning, and, finally, we will illustrate the behavior of our algorithms by a number of well-chosen numerical experiments.
Introduction
. In this paper we expand on the usage of the Jacobi-Davidson
method [26], [24] for the computation of several solutions of the generalized eigenproblem(# A-
(1)
where A and B are large and sparse (n n)-matrices, which may be complex and/or
nonnormal. We will also discuss the standard eigenproblem
(2)
Of course, with I the generalized eigenproblem reduces to a standard eigen-
problem, and we could have restricted ourselves to the generalized eigenproblem case.
However, simplifications are possible when I that help reduce the memory requirements
and the computational complexity, and some phenomena are easier to
explain.
Our algorithms are based on the Jacobi-Davidson method described in [26] and
are adapted for generalized eigenproblems (and other polynomial eigenproblems)
in [24]. We have modified the Jacobi-Davidson approach so that partial (general-
ized) Schur forms are computed. The partial Schur forms have been chosen mainly
# Received by the editors March 4, 1996; accepted for publication (in revised form) March 5, 1997;
published electronically August 7, 1998.
http://www.siam.org/journals/sisc/20-1/30007.html
Department of Mathematics, Utrecht University, P.O. Box 80.010, NL-3508 Utrecht, The
Netherlands (sleijpen@math.ruu.nl, vorst@math.ruu.nl).
Current address: ISE Integrated Systems Engineering AG, Technopark Zurich, Technopark-
strasse 1, CH-8005 Zurich, Switzerland (fokkema@ise.ch).
1 The family A - #B is called a matrix pencil, and the generalized eigenvalues #, solutions
of (1), are also called eigenvalues of the matrix pencil (cf., e.g., [30]).
for numerical stability, since they involve orthogonal bases. These bases are also useful
for deflation, another ingredient of our algorithms.
In the Jacobi-Davidson approach a low-dimensional search subspace is generated
onto which the given eigenproblem is projected. This is the standard "Rayleigh-
Ritz" procedure that also underlies the Lanczos, Arnoldi, and Davidson methods.
The small projected eigenproblem is solved by standard techniques, and this leads to
approximations for the wanted eigenvectors and eigenvalues of the given large problem.
In the Davidson method [5], the solution of a simplified correction equation is
used for the expansion of the search subspace. Following an old idea of Jacobi [11],
we can also set up a correction equation, acting in the subspace orthogonal to the
current eigenvector approximation, which defines an optimal orthogonal expansion
of the search subspace. To be more precise, if the exact value for the eigenvalue is
known, then the correction equation defines the exact eigenvector. This modification
of Davidson's method is referred to as the Jacobi-Davidson method (note that this
has nothing to do with the diagonal preconditioning that is popular in combination
with Davidson's method).
The "Jacobi" correction equation may be solved by any method of choice, and
for large problems it is often more e#cient to solve this equation only approximately
by some iterative method. The speed of convergence of this iterative method may
be improved by preconditioning, and this is the only place where preconditioning is
exploited in the Jacobi-Davidson method. It should be noted that this preconditioning
does not a#ect the given eigenproblem. By including shifts in the Jacobi-Davidson
method, and by a proper selection of the approximate eigenpair for the correction
equation, the process can be tailored to find eigenpairs close to a given target value.
More details will be given in sections 2.2 and 3.
The small projected problem is reduced to (generalized) Schur form by the QZ
method [15] or by QR [8] when I. The construction of the subspace and the
projected system too, may be viewed as iterative inexact forms of QZ and QR. For
this reason we have named our new methods JDQZ and JDQR, respectively. JDQZ
produces a partial generalized Schur form for the generalized eigenproblem: a partial
"QZ-decomposition"; JDQR generates a partial Schur form for the standard eigen-
problem: a partial "QR-decomposition".
Restarts form an essential ingredient of almost any iterative method, and also
for the Jacobi-Davidson method, either for the computation of other eigenpairs, after
one eigenpair has converged, or because of limits on the dimension of the subspaces
(memory limitations). In any case, the usual restart procedure has the disadvantage
that a subspace that may contain very useful information is replaced by one single
vector, so that much valuable information may be lost. This problem has been solved
elegantly for the Arnoldi method [28] (see also [17]), and our approach (cf. section 2.3)
is related to this (see also [26, section 5.3]). In this approach the subspace is suitably
filtered to retain as much relevant information as possible. Expansion and filtering
are used in a repetitive way.
The generalized eigenproblem in section 2 forms the heart of the paper. It is
shown how (implicit) restart techniques and preconditioning can be used in order to
get inverse-free computationally e#cient algorithms. The resulting algorithm, JDQZ,
is enhanced with a deflation technique so that several solutions for the eigenproblem
can be computed. The computation of interior eigenvalues is a normally a risky
a#air, if we want to avoid shift-and-invert operations. We will discuss a rather robust
technique based on the idea of harmonic Ritz values [16], [19], [26].
D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
Section 3 focuses on the standard eigenproblem. Of course, the standard eigenproblem
can be viewed as a simplification of the generalized eigenproblem, and in
an obvious way JDQZ simplifies to the JDQR method for standard eigenproblems.
However, we have chosen to pay slightly more attention to the standard eigenproblem
since this simplification makes it easier to discuss some computational aspects of our
algorithms. In particular, we will consider the problem of preconditioning in more
detail, and we will ponder on the observed speed of convergence using well-known
arguments.
In section 4, we illustrate the convergence behavior of JDQZ and JDQR with
numerical experiments for a number of eigenproblems. Aspects that are investigated
concern, among others, the e#ect of approximation errors on the solution of the correction
equation (section 4.1), the e#ect of preconditioning (section 4.2), multiple eigen-values
(section 4.6), interior eigenvalues (sections 4.3 and 4.7), di#erent approaches
for the construction of the projected deflated problem (sections 4.3 and 4.7), and
implicit versus explicit deflation (section 4.5).
In section 5 we have collected some conclusions.
Remark 1. All computations can be done in complex arithmetic if necessary.
An alternative for real matrices would be to use quasi-Schur forms with 2  2 blocks
on the diagonal, which can be computed in real arithmetic for real matrices. It is
possible to derive a variant of Jacobi-Davidson based on this blocked form; we will
not discuss this alternative variant here.
Remark 2. With boldface letters we indicate that variables are associated with
the large n-dimensional space, and for low-dimensional spaces we use italic letters.
We use a tilde to indicate that a quantity approximates the corresponding quantity
without a tilde:
q approximates q, etc.
The algorithms are given in Matlab style. We use the Matlab conventions when
we refer to entries in matrices and vectors. In particular, where in the algorithms new
values overwrite old ones, the tildes are deleted.
2. The generalized eigenproblem.
2.1. Preliminaries.
Convention 1. We denote a generalized eigenvalue of the matrix pair (A, B) as
a pair #. This approach is preferred because underflow or overflow for #,
in finite precision arithmetic, may occur when # and/or # are zero or close to zero,
in which case the pair is still meaning and useful [15], [21], [30, Ch.VI].
Remark 3. Observe that, for each #= 0, the pairs # and # correspond
to the same generalized eigenvalue. Rather than scaling the coe#cients of # in
our algorithms (for instance, such that # [0, 1] and we follow the
advise in [15], and we show the results as produced by the QZ algorithm: the size of #
and # may give valuable information on the conditioning of the computed eigenpair.
However, in the construction of our algorithm, the choice of some parameters leads
to an implicit scaling.
For generalized eigenproblems, a partial generalized Schur form is defined as follows

Definition 1. A partial generalized Schur form of dimension k for a matrix pair
(A, B) is the decomposition
are orthogonal (nk)-matrices, and S k and T k are upper triangular
k)-matrices. A column q i of Q k is referred to as a generalized Schur vector, and
JACOBI-DAVIDSON STYLE QR 97
we refer to a pair (q i , # i , # i #) with # i , # i #S k (i, i), T k (i, i)# as a generalized Schur
The formulation in (3) is equivalent to
Furthermore, if (x, #) is a generalized eigenpair of (S k , T k ) then (Q k x, #) is
a generalized eigenpair of (A, B).
2.2. Jacobi-Davidson. We will briefly describe Jacobi-Davidson for the generalized
eigenproblem (1); for details we refer to [24].
Similar to subspace approaches for standard eigenproblems, in each step the approximate
eigenvector
# q is selected from a search subspace span{V}. The Galerkin
condition, with associated approximate generalized eigenvalue #
#, requires orthogonality
with respect to some test subspace span{W}:
# A
For the generalized case, it is, in view of (3) and (4), natural to take the test subspace
span{W} di#erent from the search subspace: the Petrov-Galerkin approach. Search
subspace and test subspace are of the same dimension, say, j. Equation (5) leads to
the projected eigenproblem
that can be solved by conventional techniques, and a solution (u, #
#) is selected
(note that (6) is a j-dimensional problem).
Then the Petrov vector
and the residual r #
# A
# q, associated
with the Petrov value #
#, are computed. The subspaces span{V} and span{W}
are expanded in each step of the iterative process. In the variant of the Jacobi-
Davidson method used in this paper, the search subspace is expanded by a vector v
that is orthogonal to
# q and that solves (approximately) the Jacobi correction equation
# z
The essential di#erence with the Davidson approach is in the inclusion of the left
and right orthogonal projections in this correction equation. It can be shown that
for nonzero
# z # span{A# q, B# q}, and if (7) is solved exactly, the convergence of the
generalized eigenvalue will be quadratic; see [24, Th. 3.2].
In the next iteration of the algorithm, span{V, v} defines the new search subspace,
and we will explain how to expand W appropriately. Since we prefer orthogonal
matrices V and W, similar to Z and Q in (3), the new columns of V and W are
orthonormalized by modified Gram-Schmidt (or some other stable variant of Gram-
Schmidt).
We use the QZ algorithm [15] to reduce (6) to a generalized Schur form. With j
the dimension of span{V}, this algorithm yields orthogonal (j  j)-matrices U R and
U L and upper triangular (j  j)-matrices S A and S B such that
The decomposition can be reordered such that the first column of U R and the (1, 1)-
entries of S A and S B represent the wanted Petrov solution of (6) [33], [34], [7]. For a
Matlab version of this reordering algorithm, see [7, Ch. 6.C].
D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
With the decomposition in (8), we construct an approximate partial generalized
Schur form (cf. (3)): VU R approximates a Q k , and WU L approximates the associated
Z k . Since span{Z k makes sense to choose
W such that, for some scalars # 0 ,  0 with, say, |# 0 | the space span{W}
coincides with span{# This choice is also in line with the restriction
on
# z for quadratic convergence.
In summary, the proposed method has the following main ingredients:
. Form the projected system (6) and reduce it to an ordered generalized Schur
form (8). Select as approximate generalized eigenpair
# (VU R (:, 1), #S A (1, 1), S B (1, 1)#).
. Form the Jacobi correction equation:
with
#B# q,
where # is a normalization constant; the choice for # 0 and  0 will be discussed later
(in section 2.4). Compute an approximate solution v #
q of (9). Note that
VUR (:, 1)) is normalized.
. Expand V with the orthonormal complement of v and W with the orthonormal
complement of w, where
Modified Gram-Schmidt is used for the computation of the orthonormal complements.
It can be shown that, with the above choices for
# z and W,
In this approach, the relation between the partial generalized Schur form for the large
problem and the complete generalized Schur form for the small problem (6) via right
vectors is similar to the relation via left vectors
The fact that
also convenient for restart purposes, as we will see
in section 2.3.
. After convergence, expand the partial Schur form with the converged Schur
vector, and repeat the algorithm with a deflated pencil for other eigenpairs. For
details on deflation, see section 2.5.
2.3. Practical selection and implicit restart. When we have reduced the
projected eigenproblem (6) to a generalized Schur form by the QZ algorithm [15],
then we can exploit the generalized Schur form for various purposes:
- selection of a Petrov pair (# q, #,
#),
- selection of the corresponding left vector
# z (cf. (12)),
- restriction of the dimension of the subspaces span{V} and span{W} if necessary
deflation after convergence of a Petrov pair.
We will explain the first and third points in more detail in this section.
JACOBI-DAVIDSON STYLE QR 99
Suppose that the generalized Schur form of the pencil
given by
is ordered with respect to # such that
where j is the dimension of span{V}. Then
# (VU R (:, 1), #S A (1, 1), S B (1, 1)#)
is the Petrov approximation corresponding to the projected system (6) with Petrov
value closest to the target # . The corresponding left vector is given by
Furthermore, VU R (:, spans the subspace that contains the i most
promising Petrov vectors. The corresponding test subspace is given by WU L
Therefore, when we want to reduce the dimension of the subspace ("implicit restart")
to j min , j min < j, then we simply discard the columns v j min+1 through v j and w j min+1
through w j , and continue the Jacobi-Davidson algorithm with
(cf. [26, section 5.3]).
Remark 4. Our restart strategy follows similar ideas as in the implicitly restarted
Arnoldi (IRA) [28]. However, in [28] implicit shifts are used to delete the unwanted
part, instead of explicitly selecting the wanted portion of the Krylov subspace as we
do. The situation for IRA is more complicated because the reduced search subspace
has to be a Krylov subspace. For further details, see [28], [13].
2.4. The values of # 0 and  0 . We will now explain how to select the scalars
# 0 and  0 in (10). The restriction |# 0 | scaling and avoids trivial
expansions. We discuss two approaches. The first one, in section 2.4.1, can be viewed
as a generalization of the approach by Ritz values for standard eigenproblems for
optimal expansion of the test subspace. The second one, in section 2.4.2, is related to
the approach by harmonic Ritz values and aims for optimal selection of Petrov pairs.
2.4.1. Fixed values for # 0 and  0 . If v is the expansion vector for the search
subspace, then in the general setting we have to expand the test subspace by # 0 Av+
0 Bv. Note that if
# q is the new approximate eigenvector then expanding the old
search subspace by v is equivalent to expanding it by
q, so that the new test subspace
can also be obtained by expanding with # 0 A# q I, the obvious
choice would be # I, the obvious choice would be
In this case, although B# q is in the direction of
# q is close to
some eigenvector q, multiplication by B may diminish the most important eigenvector
components of
# q if the eigenvalue of B associated with q is (very) small. Therefore,
expanding the test space by B# q may be (much) less optimal than expanding by
# q. In
the presence of rounding errors, this e#ect may be even more prominent.
100 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
The value of
as function of # and , is maximal if
and
where # denotes the complex conjugate of #. This approach can be seen as an attempt
to expand the test subspace span{W} optimally in the direction of z, where z is the
normalized vector Aq (or Bq).
Since we have to choose # 0 and  0 before we know the generalized eigenvalue
#, the best we can do, in particular in the initial phase of the process, is to select
and  0
where # is the target value.
For an approach by which the choice for the scalars is made adaptively, see [7,
Ch. 6, section 3.1]. In practice, we have not seen much advantage of that approach
compared with the approach to be discussed in section 2.4.2.
2.4.2. Values for # 0 and  0 based on harmonic Petrov values. In this
section, we will introduce harmonic Petrov values. We will see that the harmonic
Petrov values that are closest to a target can be considered as extremal Ritz values
for a specific test subspace, also if the target is in the interior of the spectrum. In
particular for the computation of interior eigenvalues, the harmonic Petrov values
appear to be attractive competitors for the standard Petrov values of the approaches
in section 2.4.1: for generalized eigenproblems the costs for the computation of the
harmonic Petrov values are the same as for standard Petrov values, and, because of
the extremality property, harmonic Petrov values closest to the target appear to be
the best choices, also in early stages of the process. This is in line with observations
for the standard eigenproblems made in [16], [26].
We first consider the computation of the eigenvalues of a standard eigenproblem
I) that are close to some target value # in the interior of (the convex hull of) the
spectrum. The transformation # 1/(#) maps these eigenvalues # to extremal
eigenvalues of and in that case the "correct" eigenpair approximations
can be obtained easily.
However, we want to avoid matrix inversion. With some formula manipulation,
it can be shown that this can be achieved by taking the search subspace and the
test subspace both equal to span{(A - # I)V} (cf. [26, section 5.1]): the resulting
eigenvalue approximations
# for A are then the solutions of
The solutions
# are called harmonic Ritz values of A, with respect to # (cf. [19], [26];
they have also been used in [16]); Vu is the associated harmonic Ritz vector. Since
W and #W, with # 1/ span the same space, the harmonic Ritz values
appear as Petrov values for the test subspace generated as in (11) with
and  0 # -
JACOBI-DAVIDSON STYLE QR 101
For generalized problems, with # 0 and  0 as in (17), the Petrov values closest to the
target value correspond to absolute largest Ritz values of the standard eigenproblem
with matrix Therefore, for this generalized case also a better
selection of appropriate eigenpair approximations may be expected. We refer to the
Petrov values associated with this choice of test subspace as harmonic Petrov values.
2.5. Expansion of Schur form and deflation. In this section, we focus on
the e#cient computation of a set of generalized eigenpairs. The idea is to use the
Jacobi-Davidson method for the computation of a partial generalized Schur form as
a major step in solving generalized eigenproblems.
Suppose that we have the partial generalized Schur form AQ
and BQ We want to expand this partial generalized Schur form with
a suitable q and z to
# and
# .
From this we deduce that the generalized Schur pair (q, #) satisfies
for
This leads to
Hence, (q, #) satisfies
and the generalized Schur pair (q, #) is therefore also an eigenpair of the deflated
matrix pair
In JDQZ we solve this eigenproblem again with the Jacobi-Davidson method.
In more detail, the procedure is as follows. Let V and W be orthogonal (n  j)-
matrices such that V # Q
and
and denote the generalized Schur form of the matrix pair (M A , M B ) by
2 Note that the Jacobi correction equation (7) has a similar structure and can be derived in a
similar way.
D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
If this generalized Schur form is ordered with respect to the target value # , then
# (VU R (:, 1), #S A (1, 1), S B (1, 1)#)
is a Petrov pair approximation for a solution of (20). The corresponding left vector
is given by
The Jacobi-Davidson method expands V with the orthogonal complement of v
that is an (approximate) solution of the generalized deflated Jacobi correction equation
q. Note that Q #
also have to expand W; we expand it with the complement of (I -Z k-1 Z #
# q, orthogonal with respect to W.
When the generalized Schur pair (# q, #
#) is su#ciently close to (q, #), then
we may continue for still another generalized Schur pair. In that case V and W are
replaced by VU R (:, in order to obtain a new search subspace
orthogonal to span{Q k-1 ,
q} and a new test subspace orthogonal to span{Z k-1 ,
respectively, and we continue the process.
2.6. Solution of (deflated) correction equation. In this section we discuss
how the generalized deflated Jacobi correction equation can be solved and how pre-conditioning
is involved.
The correction equation (20) involves an operator for which the domain and the
image space di#er. This means that Krylov subspace methods cannot be applied right
away. Fortunately, this can be fixed easily by incorporating preconditioning.
For preconditioning of the correction equation (20), we propose to use
for some preconditioner K # A- # B.
We introduce the following notation.
Notation 1.
# q ], the matrix Q k-1 expanded by
z ], the matrix Z k-1 expanded by
Z k , the expanded matrix of preconditioned vectors,
Y k , the projected preconditioner
Z k .
In this notation the left preconditioned correction equation for the generalized
correction equation (24) can be written as
where
section 3.2, we give more details for the simpler
standard eigenproblem).
JACOBI-DAVIDSON STYLE QR 103
Of course, right preconditioned generalized correction equations can be derived
as well. With
where
v.
Note that for the operators in the preconditioned correction equation (26), the
domain and the image space coincide, so that Krylov subspace methods can be used.
A pseudocode for the preconditioned Jacobi-Davidson QZ algorithm with harmonic
Petrov values (as discussed in section 2.4.2) is given in Algorithm 1.
In

Table

1 we have listed the main computational ingredients per iteration of
JDQZ.

Table
The computational costs of JDQZ per iteration. The integers j and k denote the dimensions
of span{V} and span{Q}, respectively.
Part dots axpys MVs K
The correction equation Variable costs
The projected problem 6j
The Petrov approximation 2k
a If Krylov subspace methods are used to solve the correction equation,
then the products Av and Bv are often already available, as side products.
Then no MVs are needed in this part.
b Instead of computing the residual r as
#B#q, r may also be
computed as
Algorithm 1); depending on the number of nonzeros in A, B, and the value
of j, this may be more e#cient.
3. Jacobi-Davidson for the standard eigenproblem. When I, the
JDQZ algorithm simplifies greatly, and this can be used to improve computational
e#ciency. We will also consider this situation in order to discuss specific aspects
such as preconditioning. It should be noted that these simplifications only lead to a
memory e#cient algorithm when we give up working with harmonic Ritz values since
these involve skew instead of orthogonal projections in the standard setting [26]. For
the computation of exterior eigenvalues this does not pose serious problems, but for
interior eigenvalues we prefer the more robust harmonic eigenvalue approximations
(section 2.4.2) and orthogonal projections, which means that we prefer JDQZ (with
instead of the algorithm discussed in this section.
With the natural choice # we have that so that the
projected eigenproblem reduces to
that is,
# (cf. (6)).
For this low-dimensional problem we select a solution (u,
#) by standard computational
techniques. The Ritz value
# and the Ritz vector
an approximate
104 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
Algorithm 1
Preconditioned JDQZ, using harmonic Petrov values.
function [ Q,Z,RA , RB
while k < k max ,
else
- the correction equation -
Y
Solve v (approximately) from:
Y
Y
- the projected problem -
while found,
"found and implicit restart part", see Algorithm 2
JDQZ returns a partial generalized Schur form (Q, Z, RA , RB ) of dimension kmax of the
matrix pair (A, B) with generalized eigenvalues near the target # . K is a preconditioner
for A - #B, v 0 is an initial guess, and # is the stopping tolerance. jmax and j min specify
the dimension of the search subspace before and after implicit restart, respectively. qz is
a Matlab function that computes a generalized Schur decomposition. The function mgs
performs modified Gram-Schmidt and qzsort sorts the generalized Schur form. Matlab
implementations of mgs and qzsort can be found in [7, Ch. 6.A-C].
JACOBI-DAVIDSON STYLE QR 105
Algorithm 2
"Found and implicit restart part" of preconditioned JDQZ with harmonic Petrov values.
if found,
eigenvalue and eigenvector, respectively, with residual r # (A-
# q (we will assume
that # q#
For the expansion of V, we take a vector v #
q that solves (approximately) the
correction equation
The expanded search subspace is span{V, v}. In exact arithmetic, V is an orthonormal
We have used modified Gram-Schmidt in our computations
for the construction of an orthonormal basis of the search subspace.
As mentioned in the introduction, if
# is replaced in the correction equation (29)
by an eigenvalue #, then the associated eigenvector is contained in the space spanned
by V and the exact solution of the Jacobi-correction equation. If the correction
equation (29), with
#, is solved exactly, then the speed of convergence for the selected
Ritz values is asymptotically quadratical (cf. [26], [24]).
We reduce the projected eigenproblem (28) to Schur form by the QR algorithm [8],
and then we can exploit the Schur form for the selection of a Ritz pair (# q,
#) and for
restriction of the dimension of the subspace span{V} in the same way as explained
in section 2.3 (note that in this case
A reordering algorithm for the Schur form can be found, for instance, in [29],
[9], [20]; a Fortran implementation is available from LAPACK [1]. A simple Matlab
implementation for reordering with respect to a target value # is given in [7, Ch. 6.B].
3.1. JDQR. For the standard eigenproblem we can use the JD algorithm for
the computation of a partial Schur form, which can be written as (cf. [22]):
106 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
is an orthogonal (n  k)-matrix, and R k is an upper triangular (k  k)-
matrix. A column q i of the matrix Q k is a Schur vector, and the pair (q i , # i ), with
i), is a Schur pair.
The diagonal entries of the matrix R k represent eigenvalues of A, and if (x, #) is
an eigenpair of R k then (Q k x, #) is an eigenpair of A.
The resulting simplification of JDQZ will be referred to as JDQR. Although most
of the simplifications are obvious, we will give the main expressions for JDQR for ease
of reference.
Suppose that k - 1 Schur pairs have detected; i.e., we already have the partial
Schur form AQ Then the new Schur pair (q, #) is an eigenpair of
the deflated matrix
We then solve the eigenproblem for the deflated matrix (30). More precisely, the
JD algorithm for the deflated matrix (30) constructs a subspace span{V} for finding
approximate eigenpairs, and V is an orthogonal matrix such that V # Q
the deflated interaction matrix M we have
The ordered Schur form gives an approximation (# q,
for a wanted eigenpair of the deflated matrix (30). Then, according to the Jacobi-
Davidson approach, the search subspace span{V} is expanded by the orthonormal
complement of v to V, where v is the (approximate) solution of the deflated Jacobi
correction equation
# q.
Note that the projections in (32) can be subdivided into two parts: the part
associated with Jacobi-Davidson and the deflation part
Observe also that Q # k-1
Remark 5. Two deflation techniques can be found in literature for subspace
methods like Arnoldi's method. They are referred to as explicit and implicit deflation
(cf., e.g., [22, Ch. VI, section 2.3]). In explicit deflation, the computation is continued
with a deflated matrix after detection of Schur vectors. For e#ciency reasons,
A - QRQ # is used (Schur-Wielandt deflation), rather than the more stable representation
In implicit deflation, each new vector for the search
subspace is generated with A itself and is then made orthogonal to the detected Schur
vectors before adding it to the search subspace. Our approach is a mixture of both
techniques. In the Jacobi correction equation we use the explicitly deflated matrix.
Since the solutions of the deflated correction equations are orthogonal to the detected
Schur vectors, there is no need to use the deflated matrix for computing the deflated
interaction matrix M ; we compute M as
Similar observations hold for the interaction matrices MA and MB in the generalized
case (cf. (22)).
Exclusively implicit deflation is possible as well: solve the correction equation
approximately with the nondeflated A and make the resulting solution orthogonal
JACOBI-DAVIDSON STYLE QR 107
to the detected Schur vectors. In this approach we avoid expensive matrix-vector
multiplications, but explicit deflation appears to improve the condition number of the
linear system, and that leads to a faster converging process for the Jacobi correction
equation (29). The decrease in the number of iteration steps, for the correction
equation, appears often to compensate for the more expensive multiplications (for a
numerical illustration of this, see section 4.5).
Moreover, the explicitly deflated correction equation (32) appears to lead to more
stable results. This can be understood as follows. Without deflation the resulting
solution of the correction equation may have a significant component in the space
spanned by the detected Schur vectors. By subtracting this component (as in implicit
cancellation may occur. If we work with an explicitly deflated matrix, such
cancellation is avoided.
Remark 6. As in implicitly deflated Arnoldi methods, the accuracy of an approximate
Schur pair in our method not only depends on the norm of the residual and
on the condition number of the pair but also on the approximation errors in the previously
detected Schur pairs (cf., e.g., [22, Ch. IV, section 2.5] and [13, section 6.4.1]):
in the derivation of the algorithms it is assumed that V #
which is true for exact Schur vectors. In practice, span{AQ} will not be contained in
span{Q}.
3.2. Preconditioning. In this section we will discuss preconditioning for the
correction equation. Preconditioning is not straightforward because of the
projections involved. We will derive explicit expressions for left and right preconditioned
correction equations.
In each iteration step we need to solve a deflated Jacobi correction equation (32)
for a given
# q and
(cf. (32)). For the approximate solution of this equation we
may use a Krylov subspace method, e.g., GMRES [23], or BiCGstab(#) [25]. The
rate of convergence and the e#ciency of Krylov subspace methods is often improved
by preconditioning. The identification of an e#ective preconditioner may be a prob-
lem. For instance, for interior eigenvalues the construction of an e#ective incomplete
LU-factorization [14], [10] for A-
#I may require much fill-in, 3 which makes the construction
expensive. As we will argue in section 3.3, it may be a good strategy to
compute a good (and possibly expensive) preconditioner K for A - #I for one fixed
value of # only and to use
as the preconditioner for various
# q and
#. Note that the projections on K are necessary
to let
K operate on the proper subspace (cf. [24]).
We will now give some more details on the derivation of the expressions for the
preconditioned correction equation. We use the same notation as in Notation 1 with
. The typical usage of the preconditioner in a Krylov subspace method would
look like
solve t, with
for s, with
3 These incomplete factorizations have not been designed for linear systems related to eigenprob-
lems. The solutions for which these factorizations are most e#ective are usually rather smooth, which
means that components of slowly varying eigenvectors are favored by the preconditioning.
108 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
The following lemma gives us an explicit expression for the solution t of (34) in
terms of easily computable matrix-vector products with
k and K -1 . The lemma
generalizes Proposition 7.5 in [24], and the proof runs along the same lines (for details,
see [7, Ch. 6, section 2.4]). Note that
H k is of small dimension, so it is cheaply inverted.
There is no need to invert K explicitly; instead can be computed by solving
v from
Lemma 1. If
H k is nonsingular, then the solution t of equation (34) is given by
Remark 7. If one stores the matrix Y k-1 # K of preconditioned
Schur vectors, one only has to compute the last column K -1
q of the matrix
q ] at each iteration step. Furthermore, when storing the projected preconditioner
only the last column and last row of
have to be computed in an iteration step.
Remark 8. If the preconditioner K is indefinite, then the matrix
may become
singular for an "unlucky" choice of approximate Ritz pair (# q,
#). This causes a
breakdown, but it never happened in our experiments. The breakdown may be cured
by selecting a di#erent nearby approximating Ritz pair (# q # ,
temporarily for the
current Jacobi-Davidson iteration.
preconditioning. From Lemma 1, it follows that the left preconditioned correction
equation is equivalent with
where
Note that the projection has to be applied explicitly to the residual. For the
unpreconditioned case there was no need for explicit projection, since there the fact
that the residual is associated with a deflated matrix and with a Ritz pair implied
orthogonality to
Observe that, for equation (36) is equivalent to the one in (32).
Of course, right preconditioned correction equations, similar to (27), can be derived
in a corresponding manner (for details, see [7, Ch. 6, section 2.4]).
Remark 9. If one uses Krylov subspace methods for solving the second equation
in (36), then one encounters matrix-vector products of the form
with t of the form
obviously,
and for the
approximate solution v we have that
provided that this is the case for the
initial guess as well. Moreover, the projection
k ) in front of t in (37) is
then redundant, and (37) reduces to
t.
JACOBI-DAVIDSON STYLE QR 109
Algorithm 3
Preconditioned JDQR.
function
while k < k max ,
else
- the correction equation -
Y
Solve v (approximately) from:
Y
Y
- the projected problem -
while found,
"found and implicit restart part", see Algorithm 4
JDQR returns a partial Schur form (Q, R) of the matrix A of dimension kmax with eigen-values
near the target # . K is a preconditioner for A - # I, v 0 is an initial guess, and #
is the stopping tolerance. jmax and j min specify the dimension of the subspaces V before
and after implicit restart, respectively. schur is a Matlab function that computes a Schur
decomposition. The function mgs performs modified Gram-Schmidt, while qrsort sorts the
Schur form. For Matlab implementations of mgs and qrsort, see [7, Ch. 6.A-B].
A pseudocode for the preconditioned Jacobi-Davidson QR algorithm is given in
Algorithm 3.
In

Table

2 we have listed the main computational ingredients per iteration of
JDQR.
D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
Algorithm 4
"Found and implicit restart part" of JDQR.
if found,

Table
The computational costs of JDQR per iteration. The integers j and k denote the dimensions
of span{V} and span{Q}, respectively.
Part dots axpys MVs K
The correction equation Variable costs
The projected problem 3j
a If Krylov subspace methods are used to solve the correction equation,
then the product Av is often already available as a side product. No MV
is needed in this part then.
b Instead of computing the residual r as A#q-#q, r may also be computed
as
#q, where V A # AV (cf. Algorithm 3); depending on the
number of nonzeros in A and the value j, this may be more e#cient.
3.3. The quality of the deflated preconditioner. Even when the preconditioner
K is constructed for a fixed # , then the correction equation still involves
projections that become more expensive after each Schur pair that has been detected,
but this does not necessarily lead to a more expensive computational process (com-
pared with explicit restart). When iterative solvers are used, they may converge faster
because the field of values of the projected operator
is contained in the field of values of A-
# I, and that may be smaller, especially after
exterior eigenvalues have been detected.
The projections may also have a positive e#ect on the preconditioner.
JACOBI-DAVIDSON STYLE QR 111
We see that, on the one hand, the preconditioning error is enlarged by a small shift
(#)I, but on the other hand, the projections diminish the error by filtering out
the detected Schur vectors. If the error R is large with respect to eigenvectors corresponding
to eigenvalues near # , then the projected error
will be significantly smaller, and the only penalty is a (small) shift due to #. It
seems plausible (cf. [32, Ch. IV]) that this will not lead to a significantly less e#ective
preconditioner, and it may help to explain the e#ectiveness of a fixed preconditioner
for JDQR in some of our experiments.
3.4. Notes on the speed of convergence. In this section we will make some
comments with respect to the convergence behavior of JDQR. We will use well-known
arguments.
The JDQR algorithm has nice properties with respect to the overall performance.
While adjusting for one Schur pair, the subspace span{V} also accumulates components
for other Schur pairs. As a result, after one Schur pair has been detected, other
Schur pairs may follow more quickly than after a complete restart. These components
will appear in a similar way as for the shift-and-invert Arnoldi [22] process, with a
shift
# for a (deflated) eigenproblem, as can be understood as follows.
For simplicity, suppose that A has a complete set of eigenpairs
and that we are trying to find an approximation (# q,
#) for
The exact solution of (29) is given by
with
(cf. [26, section 4.1]). Writing
# q as
that
We may assume, without loss of generality, that # i #= 0, because
# q is a Ritz vector
which means that #
#/2. The latter case is unlikely to
happen, due to rounding errors, and the first case indicates full convergence.
Hence, eigenvector components corresponding to eigenvalues closer to
# will be
amplified more in (A-
# q. The component orthogonal to
# q is used as an expansion
for V and thus as soon as
q has a large component in the direction of x 1 , say
that the angle is less than #/4, then necessarily components other than x 1 become
dominant. That is,
In Fig. 1 we have illustrated this phenomenon. The bullets represent the amplification
factors 1/|# i - #| for components in the direction of x
112 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
THETA3Fig. 1. Amplification factors of eigenvectors.
In the subsequent iterations similar amplifications will occur, and the closer # i is
to
# the more rapid the angle
This argument is repetitive: if the angle
very small, then the
corresponding # 2 will be very small and other components, due to orthogonalization,
will become more dominant.
Consequently, while the process converges to a Schur pair, the search subspace V
will provide good initial approximations for the nearby Schur pairs. Moreover, slow
convergence during one stage may be compensated for by faster convergence in the
next stage, because the subspace span{V} will be enriched with more components of
other Schur pairs due to repeated amplifications. This is observed in our numerical
experiments; see section 4.
4. Numerical experiments. In this section we present numerical results, obtained
with JDQZ and JDQR, for several generalized eigenproblems and standard
eigenproblems. The purpose of these experiments is to get an impression of the actual
behavior of these methods. We have not tried to find the most e#cient parameter
choices for each particular problem. We will illustrate the e#ect of more accurately
solving the correction equation and the e#ect of including appropriate precondition-
ing. We will show that the harmonic Petrov value choice for the test subspace may
lead to superior convergence behavior, not only for the generalized eigenproblem but
also for the standard eigenproblem. We will demonstrate that the projections in the
correction equation (32), involving detected Schur vectors, are essential components
of the algorithms. We will also consider eigenproblems where multiple eigenvalues are
involved.
The computations were done in double complex precision (# 15 digits) on a Sun
workstation. To facilitate comparison, we have selected for all cases j
(the dimension of the subspace before and after implicit restart, respectively), and
a fixed random real vector v 0 as an initial guess (cf. Algorithms 3 and 1).
As iterative solvers for the correction equation, we have considered full GMRES
[23] with a maximum of m steps, denoted by GMRESm , and BiCGstab(2) [25].
For BiCGstab(2) a maximum of 100 matrix multiplications was allowed. Of course,
this is comparing apples with pears; our main purpose was to mimic two realis-
JACOBI-DAVIDSON STYLE QR 113
tic scenarios, one with a fixed number of matrix-vector multiplications (GMRESm ),
and one with an iterative method to satisfy some stopping criterion (BiCGstab(2)).
As stopping criterion for the iterative methods for the correction equation, also for
GMRESm , we have used # r
is the initial residual,
is the
residual corresponding to the approximate solution produced by the inner method,
and j is the iteration number for the current eigenvalue approximation in the outer
iteration. Hence, as the outer iterations proceed, the inner iterations are solved more
accurately. This choice was inspired by the fact that the Jacobi-Davidson method
may be viewed as a Newton process [24], [27], and for Newton processes this stopping
criterion may lead to e#cient algorithms [6]. As the initial guess for the inner
iteration method we always took the null vector.
In the figures of the convergence behavior for JDQZ and JDQR, the performance
is plotted in terms of the actual amount of work, in millions of floating point operations
(flops), versus log 10 of the residual norm. The reason for this is that the computational
work in JDQZ and JDQR consists of two parts of a di#erent nature: one part is for the
inner iteration process, in which a correction equation is (approximately) solved; the
other part is for the outer iteration, in which an approximation for the (generalized)
Schur pair is constructed. If in the inner iteration the correction equation is solved
more accurately, then the number of outer iterations may decrease. Therefore, it
would be misleading to monitor the total number of matrix multiplications. It might
give a bad impression of the total costs, because most of the matrices are sparse and
therefore the dot products and vector updates in the outer and the inner iteration
represent substantial costs in JDQZ and JDQR.
Furthermore, we have plotted the entire convergence behavior. This means that
the convergence history of the residuals of all subsequentially selected approximate
eigenpairs is plotted. Whenever the residual norm curve drops below the acceptation
level, indicated by the dotted horizontal line, an eigenvalue is accepted and the search
process for the next one is continued. A large residual norm in the step immediately
after acceptance marks the start of a new search.
4.0.1. Construction of suitable initial subspaces. Specifically in the first
few steps of the process the Ritz or Petrov vectors are usually poor approximations
of the wanted eigenvectors, and the target value # may be relatively (much) closer to
the wanted eigenvalues than any of the approximate eigenvalues. In these cases, the
correction equations (26) and (36) lead to relatively poor expansions of the search
subspace. To see this, recall that the wanted eigenvector would be in the new search
subspace if this space would have been expanded by the exact solution for the correction
equation with the wanted eigenvalue instead of
# (cf. section 3). This observation
indicates how to improve the expansion in the first few steps: take in the correction
equation # instead of
#. To detect whether
# is close enough to replace # , we monitor
the norm of the residual: we take
#,
#) instead of #, 1#) in the correction equation
as soon as the first residual norm drops below a threshold value # tr . A similar
switch was proposed in [18].
Moreover, in all experiments we used GMRES 1 for the first j min iterations in order
to build up a search subspace span{V} in a relatively inexpensive way. Especially
when a preconditioner is involved, this approach can be justified with arguments
similar to those in the preceding paragraph (cf. [24, section 9.4]).
4.0.2. Stopping criterion. In our experiments, we considered an approximate
eigenpair
#) converged if the residual r is su#ciently small (#r#); then (# q,
114 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST

Table
Four eigenvalues of DW1024, computed by JDQR (cf. section 4.1).
9.6473e - 01
9.6551e - 01
9.7780e - 01
9.7880e - 01
log10
of
residual
norm
number of flops x 1e6
Fig. 2. Convergence history for DW1024, showing the e#ect of solving the correction equations
more accurately (cf. section 4.1).
is a "detected" eigenpair. In the algorithms, weighted residuals (e.g., #r#A#) or
more sophisticated stopping criteria can be employed as well (cf. [4]).
In our experiments we varied values for some of the parameters in JDQZ and
JDQR. For easy reference, we recall their meaning:
Parameter Description
# the target value
k max the number of wanted Schur pairs
# the stopping tolerance in the outer iteration
# tr threshold used for building initial subspaces
4.1. The influence of the correction equation. The purpose of this example
is to show the e#ect of a more accurate solution of the correction equation. We consider
the square dielectric waveguide standard eigenproblem DW1024 of order 2048 [2]. The
problem comes from an integrated circuit application. The rightmost eigenvalues and
their eigenvectors are wanted. We have used JDQR with standard Ritz values.
We took we have not used
preconditioning.
The computed eigenvalues are given in Table 3. The convergence history is plotted
in Fig. 2 for JDQR, GMRES 1 , and GMRES 10 . A summary of the number of iterations,
the number of matrix multiplications (MVs), and the number of flops is given in

Table

4.
JACOBI-DAVIDSON STYLE QR 115

Table
Summary of results for DW1024 (cf. section 4.1).
Method for the correction
equation
JDQR
iterations
MVs flops

Table
Five eigenvalues of BWM2000, computed by JDQR (cf. section 4.2).
2.4427e - 07 - 2.1395e
2.4427e
When solving the correction equation more accurately, the number of MVs is
increased, but the number of outer iterations is reduced significantly (see Table 4),
resulting in a much better overall performance.
With GMRES 1 the search subspace is the span of the residuals, and in that case
JD (with implicit restart) generates the same subspaces as IRA [28]. The eigenvalues
are not well separated in this case, and therefore Arnoldi converges only slowly. This
explains the poor convergence of JDQR with GMRES 1 .
Note that after an initial phase with two small bumps, JDQR converges quite
fast. For the next eigenvalues there is no such initial stagnation. Apparently, in the
iterations for the first eigenvalue, components for the next Schur vectors are already
collected in span{V} (cf. section 3.4).
4.2. The e#ect of preconditioning. When increasing the number of steps
in GMRES, the correction equation will be solved more accurately, and the number
of outer iterations may decrease as we have seen. But sometimes we need too many
inner iterations with GMRES for acceptable convergence of the outer iterations. With
appropriate preconditioning, we may see a dramatic improvement, as we will see in
this example.
We consider the standard eigenproblem BWM2000 of order 2000 for the Brusse-
lator wave model [2], [22]. The problem models the concentration waves for reaction
and transport interaction of chemical solutions in a tubular reactor. Our task is to
determine the eigenvalues with the largest real part in order to verify whether their
real parts are positive or negative (corresponding to stable or unstable modes). Again
we have used JDQR with standard Ritz values.
For this problem, we have selected
The computed eigenvalues are listed in Table 5. The convergence history is plotted
in Fig. 3 for JDQR with unpreconditioned GMRES 10 and with
preconditioning. A summary of the results is given in Table 6.
From Fig. 3 we see that JDQR with GMRES 10 does not converge (we checked
even up to GMRES 50 with little or no improvement), but with preconditioning JDQR
performs rather well. Again we see that the speed of convergence for the first eigenvalue
is somewhat slower than the speed of convergence for the other eigenvalues.
Note that, although the projections in the correction equation become more expen-
D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
log10
of
residual
norm
number of flops x 1e6
GMRES +ILU(0)Fig. 3. Convergence history for BWM2000, illustrating the e#ect of including preconditioning
in the solver of the correction equation (cf. section 4.2).

Table
Summary of results for BWM2000 (cf. section 4.2).
Method for the correction equation
JDQR
iterations
MVs flops
sive after each detected eigenvalue, the computational work for each eigenvalue is
roughly constant, except for the first eigenvalue.
It should be noted that popular software packages for solving eigenproblems have
been reported to fail for this problem BWM2000 [12].
The preconditioned correction equation (26) requires the computation and storage
of the term
Notation 1, with
If K is a good approximation
for A - #I, then
spans an
invariant subspace of A, then, in this nearly ideal case, we have that
and for this reason it is tempting to use (26) with
instead of
k . However,
this approach does not work well: it may even lead to slower convergence. In some
of our experiments (as for instance with QH882 as discussed below in section 4.3) we
lost convergence completely; in other experiments it had little or no e#ect.
4.3. Harmonic Ritz values. The JDQR algorithm computes a partial Schur
form for the standard eigenproblem with standard Ritz pairs for the Schur pairs.
However, with JDQZ for I, we can also compute a partial Schur form for the
standard eigenproblem with harmonic Ritz pairs. Here we give an example that
illustrates the improved convergence behavior with harmonic Ritz values.
We consider the Quebec Hydroelectric Power System problem QH882 of order
882 [2]. This matrix represents the Hydro-Quebec power system's small-signal model.
JACOBI-DAVIDSON STYLE QR 117

Table
Five eigenvalues of QH882, computed by JDQR (cf. section 4.3).
-4e+06
-2e+062e+06-8e+06 -6e+06 -4e+06 -2e+06 0
imaginary
axis
real axis
Fig. 4. Spectrum of QH882 (cf. section
imaginary
axis
real axis
eigenvalues
target
Fig. 5. Part of the spectrum of QH882.
The eigenvalues # of interest are the eigenvalues in the box -300 < Re(#) < 100,
in the complex plane.
For this problem, we have selected
The computed eigenvalues are given in Table 7. The convergence history is plotted
in Fig. 7 for JDQR with standard Ritz values and for JDQZ with the harmonic choice
(cf. section 2.4.2).
This problem is rather di#cult: the eigenvalues in the neighborhood of # are in
the interior of the spectrum; see Figs. 4 and 5. For all three methods, the correction
equation was solved with GMRES 20 and was preconditioned with the exact inverse
of A- # I. A summary of the results is given in Table 7.
Although the computational complexity of JDQR is less than the computational
complexity of JDQZ (cf. Tables 1 and 2), it is not the most e#cient method here. From
the irregular convergence behavior of JDQR in Fig. 6 we may conclude that JDQR
has problems in selecting the "correct" Ritz pairs and as a result the convergence is
delayed. Eventually JDQR loses track completely and stagnates. The peaks in the
convergence behavior show that sometimes the Ritz pair that is selected in the JDQR
process does not correspond to the close-by Schur pair. As a result the search subspace
is expanded in a poor direction. Clearly, for this example, this may lead to failure
of convergence. As anticipated (cf. section 2.2), JDQZ with the harmonic choice of
test subspace makes better selections, as is indicated by the smooth convergence, and
hence, its performance is much better.
4.4. Identification of suitable Ritz values (tracking). From a computational
point of view JDQR is less complex than JDQZ, even if the savings that
standard eigenproblems problem allow are exploited in JDQZ. Therefore, it may
be attractive to identify e#cient strategies for avoiding "incorrect" selection of Ritz
D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
of
residual
norm
JDQR
number of flops x 1e6
JDQZ Harmonic
Fig. 6. Convergence history for QH882. Although QH882 is a standard eigenproblem, for
computing interior eigenvalues it is more e#cient to use JDQZ with test subspaces that are di#erent
from the search subspaces (bottom picture). A better selection with harmonic Ritz values (bottom
picture) appears to compensate for a less optimal expansion of the test subspace (cf. section 4.3).
pairs. With correctly selected Ritz pairs convergence may be expected to be less irregular
and stagnation may be completely avoided. We will discuss such a strategy
now.
If the Ritz vector in the previous iteration is already a fair approximation, then
the norm of the residual gives information on the selected Ritz vector in the current
step: in case of a poor selection, the new residual can be much larger than the previous
one. It would then require additional computational work to find a Ritz pair with
small residual norm (and still be close enough to the target # ). A cheap alternative in
this case is to select a Ritz value that is close to a previously accepted one (and forget
about # ). In the experiment that we will discuss below, we have replaced in such
cases the target by the Ritz value that is selected and accepted in the previous step,
where we consider a Ritz value acceptable if the associated residual is smaller than
some specified threshold. In the example below, we took this threshold equal to # tr ,
the value used in the criterion for leaving the initialization stage (see section 4.0.1).
After convergence of the Ritz pair, the original target value is restored at the start of
the computation for the next eigenpair.
This tracking strategy does not require any additional computational costs per
step, while it appears to reduce the number of steps significantly, as we will see
below.
Because of the rather regular convergence of Ritz values to exterior eigenvalues,
improvement for these eigenvalues may not be expected with the above tracking strategy
for standard problems.
Here we illustrate (in Fig. 7) the e#ects one may see by including tracking for
JDQR. We applied JDQR to the example of section 4.3 with the same choice of
parameters.
Nevertheless, we see that the convergence behavior of JDQZ with the harmonic
choice of test subspace is still superior. See also Table 8.
4.5. The influence of Q k and Z k in the correction equation. In this
example we show that the projections with detected Schur vectors (cf. (24)) are very
JACOBI-DAVIDSON STYLE QR 119
of
residual
norm
JDQR
number of flops x 1e6
JDQZ Harmonic
Fig. 7. Convergence history for QH882 obtained with the tracking strategy (for all variants).
For interior eigenvalues, the tracking strategy improves JDQR significantly (compare the present
upper picture with the upper picture in Fig. 7), while there is no improvement for JDQZ (compare
the two bottom pictures of the present figure and Fig. 7) (cf. section 4.4).

Table
Summary of results for QH882 (cf. sections 4.3 and 4.4).
Method Iterations MVs flops
JDQR with tracking (cf. section 4.4) 99 1482 1.221e
JDQZ Harmonic 57
essential in the correction equation (cf. section 3.4), and we show what happens when
these projections are neglected. Note that we still take the Jacobi projections (with
q and
# z) into account.
We consider the bounded fineline dielectric waveguide generalized eigenproblem
BFW782 [2] of order 782. This problem stems from a finite element discretization
of the Maxwell equation for propagating modes and magnetic field profiles of a rectangular
waveguide filled with dielectric and PEC structures. The resulting matrix
A is nonsymmetric and the matrix B is positive definite. Of special interest are the
generalized eigenvalues # with positive real part (i.e., Re(# 0) and their
corresponding eigenvectors.
For this problem, the parameters were set to
The spectrum of this matrix pair is shown in Fig. 8. A magnification of the region
of interest is plotted in Fig. 9. The computed generalized eigenvalues, represented as
#, are given in Table 9. With JDQZ we discovered all four positive generalized
eigenvalues.
The convergence history, for the harmonic version of JDQZ with GMRES 10 , is
plotted in the upper picture of Fig. 10. A summary of the results is given in Table 10.
We see that JDQZ converges quite nicely.
In the bottom picture of Fig. 10 the convergence behavior of JDQZ is given for
the case where the correction equation (24) is solved without taking into account the
projections involving Q k and Z k . Of course, the correction equations that are used
120 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
-8000
-6000
-4000
-200020006000-3e+06 -2e+06 -1e+06 0e+00
imaginary
axis
real axis
Fig. 8. Spectrum of BFW782 (cf. section
4.5).
-0.4
imaginary
axis
real axis
eigenvalues target
Fig. 9. Part of the spectrum of BFW782.

Table
Five generalized eigenvalues of BFW782, computed by JDQZ (cf. section 4.5).
2.5233e
of
residual
norm
GMRES Harmonic-10
number of flops x 1e6
GMRES Harmonic without Q and ZFig. 10. Convergence history for BFW782 with (upper picture) and without (bottom picture)
deflating the matrices in the correction equations, with respect to the detected Schur vectors. (cf.
section 4.5).

Table
Summary of results for BFW782 (cf. section 4.5).
Method for the correction
equation
JDQZ
iterations
MVs flops
JACOBI-DAVIDSON STYLE QR 121

Table
generalized eigenvalues of AC1331, computed by JDQZ (cf. section 4.6).
include the rank-one projections involving
# q and
# z: these projections are essential for
Jacobi-Davidson. Furthermore, deflation in this case is realized by making the approximate
solution of the correction equation orthogonal to the detected Schur vectors
with modified Gram-Schmidt. By doing the latter twice, the overall performance is
improved significantly: in the results shown here (cf. Fig. 10) modified Gram-Schmidt
is applied twice.
However, as explained in section 3.4, we do not benefit from an improved operator
in the inner iteration. Although the resulting algorithm is computationally cheaper,
Fig. 10 shows that this does not lead to an overall better performance: the speed of
convergence becomes increasingly slower and even stagnates eventually.
4.6. Multiple eigenvalues. We consider the eigenproblem
with Neumann boundary conditions on the cube [0, 4] 3 . Finite element discretization
of this equation on an 11  11  11 regular grid, with tetrahedral elements and linear
interpolation functions, leads to a generalized eigenproblem of order 1331 (AC1331).
It has one positive generalized eigenvalue # relatively close to zero (i.e., # 0).
The other generalized eigenvalues are also positive and may be doublets or even
triplets.
For this problem, the parameters were set to
The computed 15 leftmost generalized eigenvalues represented as # are given
in

Table

11. The residual norm versus the number of flops is plotted in Fig. 11 for
the harmonic version of JDQZ with GMRES 10 and with BiCGstab(2), respectively.
A summary of the results is given in Table 12.
From the plots we see the e#ect that multiple generalized eigenvalues may have
on the convergence behavior. JDQZ converges initially quite fast until the point
that it "discovers" that the generalized eigenvalue is actually double or triple. The
convergence speed stagnates for a few iterations (two or three peaks in the plot with
GMRES and a plateau in the plot with BiCGstab(2)), after which the eigenvalues are
discovered quickly one after another. This behavior is in agreement with section 3.4:
during the stagnation phase components of other Schur vectors are amplified in the
inner iteration and collected in the search subspace, leading to faster convergence for
the next Schur pairs. The stagnation can be explained by the fact that with rank-one
Jacobi projections the correction equation may become (nearly) singular when
selecting Petrov approximations for multiple generalized eigenvalues. The iterative
methods used for solving the correction equation often su#er from this (see also [31]).
(Variable) block versions of the correction equation that take this multiplicity into
account may be preferable in such cases, but this falls outside the scope of this paper.
122 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
of
residual
norm GMRES Harmonic-10
log10
of
residual
norm
number of flops x 1e6
Harmonic
Fig. 11. Convergence history for AC1331: Stagnation followed by fast detection of triple generalized
eigenvalues (cf. section 4.6).

Table
Summary of results for AC1331 (cf. section 4.6).
Method for the correction
equation
JDQZ
iterations
MVs flops
4.7. Harmonic Petrov values for generalized problems. Our last example
shows again that for interior generalized eigenvalues the harmonic version JDQZ is
quite powerful.
We consider the MHD416 generalized eigenproblem of order 416 [2], [24], [3]. This
problem stems from a magnetohydrodynamics (MHD) model, where the interaction
of hot plasma and a magnetic field is studied. The matrix A is non-Hermitian and the
matrix B is Hermitian positive definite. Our goal is to compute interior generalized
eigenvalues corresponding to the so-called Alfven branch of the spectrum; see Figs. 12
and 13.
For this problem, the parameters were set to
and #
The computed generalized eigenvalues are plotted in Fig. 14. The convergence
history for the harmonic version of JDQZ, with GMRES 1 , are plotted in Fig. 15. The
exact inverse of A- # B (for # fixed) was used as a preconditioner for all eigenvalues.
For all generalized eigenvalues the rate of convergence is almost the same: in the
computation for one Schur pair, the search subspace apparently accumulates components
for the next Schur pairs as well.
JACOBI-DAVIDSON STYLE QR 123
imaginary
axis
real axis
Fig. 12. Spectrum of MHD416 (cf. section
4.7).0.20.61
imaginary
axis
real axis
Fig. 13. Alfven branch of MHD416 (cf.
section 4.7).0.50.60.7
imaginary
axis
real axis
eigenvalues target
Fig. 14. 20 generalized eigenvalues computed by JDQZ for MHD416 (cf. section 4.7).
5. Conclusions. We have proposed two algorithms, JDQZ and JDQR, for computing
several selected eigenpair approximations for generalized and standard eigen-
problems, respectively. The methods are based on the Jacobi-Davidson method and
compute iteratively a partial (generalized) Schur form with (generalized) eigenvalues
near a user-specified target value. For both methods, no exact inversion of any matrix
is strictly necessary, so that they are suitable for solving large eigenproblems.
Fast convergence is obtained with a projected correction equation that is solved
(approximately) by iterative methods with appropriate preconditioning. The convergence
of JDQZ and JDQR is asymptotically quadratical if this correction equation is
solved exactly. Furthermore, while converging to a particular Schur pair, the search
subspace accumulates components of other Schur pairs with (generalized) eigenvalues
near the target as well. This usually leads to faster convergence for the next
eigenpairs.
The dimension of the involved subspaces can be controlled by an e#cient implicit
restart technique in such a way that the most relevant part of the subspace is
maintained at restart.
124 D. R. FOKKEMA, G. L. G. SLEIJPEN, AND H. A. VAN DER VORST
-5
log10
of
residual
norm
number of flops x 1e6
GMRES HarmonicFig. 15. Convergence history of JDQZ for MHD416. Harmonic Petrov values allow a good
selection of Petrov pairs for computing interior generalized eigenvalues also for generalized problems
(cf. section 4.7).
The algorithms incorporate simple mechanisms for selecting the wanted eigenpair
approximations. Also multiple (generalized) eigenvalues can be detected.
Whereas in the Jacobi-Davidson method the test subspace can be chosen arbi-
trarily, in the JDQZ algorithm essentially two choices for the test subspace remain:
the standard Petrov value choice and the harmonic Petrov value choice. It is argued
and confirmed by our experiments that especially for interior eigenvalues the harmonic
approach is also superior for generalized eigenproblems.

Acknowledgments

. We acknowledge helpful discussions with Rich Lehoucq and
with Beresford Parlett (on naming conventions). We also appreciate the helpful comments
by the referees.



--R


A Test Matrix Collection for Non-Hermitian Eigenvalue Problems


The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real symmetric matrices

Subspace methods for linear
The QR transformation: A unitary analogue to the LR transformation
Matrix Computations
A class of first order factorizations methods
Ueber ein leichtes Verfahren
An Evaluation of Software for Computing Eigenvalues of Sparse Nonsymmetric Matrices
Deflation techniques within an implicitly restarted iter- ation
An iterative solution method for linear systems of which the coe
An algorithm for generalized matrix eigenvalue problems
Computing interior eigenvalues of large matrices

Generalizations of Davidson's method for computing eigen-values of sparse symmetric matrices
Approximate solutions and eigenvalue bounds from Krylov subspaces
Perturbation bounds for means of eigenvalues and invariant subspaces
Rational Krylov algorithms for nonsymmetric eigenvalue problems II.
Numerical Methods for Large Eigenvalue Problems
GMRES: A generalized minimum residual algorithm for solving nonsymmetric linear systems


A Jacobi-Davidson iteration method for linear eigenvalue problems
The Jacobi-Davidson method for eigenvalue problems and its relation to accelerated inexact Newton schemes
Implicit application of polynomial filters in a k-step Arnoldi method
Algorithm 406 HQR3 and EXCHNG: Fortran subroutines for calculating and ordering eigenvalues of a real upper Hessenberg matrix
Matrix Perturbation Theory
The rate of convergence of conjugate gradients
Preconditioning by Incomplete Decompositions
A generalized eigenvalue approach for solving Ricatti equations
Algorithm 590
--TR

--CTR
James H. Money , Qiang Ye, Algorithm 845: EIGIFP: a MATLAB program for solving large symmetric generalized eigenvalue problems, ACM Transactions on Mathematical Software (TOMS), v.31 n.2, p.270-279, June 2005
Ivo Bleylevens , Ralf Peeters , Bernard Hanzon, Efficiency improvement in an nD systems approach to polynomial optimization, Journal of Symbolic Computation, v.42 n.1-2, p.30-53, January, 2007
Lorenzo Valdettaro , Michel Rieutord , Thierry Braconnier , Valrie Frayss, Convergence and round-off errors in a two-dimensional eigenvalue problem using spectral methods and Arnoldi-Chebyshev algorithm, Journal of Computational and Applied Mathematics, v.205 n.1, p.382-393, August, 2007
Peter Arbenz , Martin Beka , Roman Geus , Ulrich Hetmaniuk , Tiziano Mengotti, On a parallel multilevel preconditioned Maxwell eigensolver, Parallel Computing, v.32 n.2, p.157-165, February 2006
Peter Arbenz , Roman Geus, Multilevel preconditioned iterative eigensolvers for Maxwell eigenvalue problems, Applied Numerical Mathematics, v.54 n.2, p.107-121, July 2005
Daniel Kressner, Block algorithms for reordering standard and generalized Schur forms, ACM Transactions on Mathematical Software (TOMS), v.32 n.4, p.521-532, December 2006
Richard Tran Mills , Andreas Stathopoulos , Evgenia Smirni, Algorithmic modifications to the Jacobi-Davidson parallel eigensolver to dynamically balance external CPU and memory load, Proceedings of the 15th international conference on Supercomputing, p.454-463, June 2001, Sorrento, Italy
P.-A. Absil , C. G. Baker , K. A. Gallivan, A truncated-CG style method for symmetric generalized eigenvalue problems, Journal of Computational and Applied Mathematics, v.189 n.1, p.274-285, 1 May 2006
Computing smallest singular triplets with implicitly restarted Lanczos bidiagonalization, Applied Numerical Mathematics, v.49 n.1, p.39-61, April 2004
Lawrence K. Saul , Sam T. Roweis, Think globally, fit locally: unsupervised learning of low dimensional manifolds, The Journal of Machine Learning Research, 4, p.119-155, 12/1/2003
James R. McCombs , Andreas Stathopoulos, Parallel, multigrain iterative solvers for hiding network latencies on MPPs and networks of clusters, Parallel Computing, v.29 n.9, p.1237-1259, September
linear eigenproblems;generalized eigenproblems;QZ-algorithm;jacobi-davidson;generalized Schur form;schur form;iterative methods;QR-algorithm