An Inverse Free Preconditioned Krylov Subspace Method for Symmetric Generalized Eigenvalue Problems.
In this paper, we present an inverse free Krylov subspace method for finding some extreme eigenvalues of the symmetric definite generalized eigenvalue problem x$. The basic method takes a form of inner-outer iterations and involves no inversion of B or any shift-and-invert matrix $A-\lambda_0 B$. A convergence analysis is presented that leads to a preconditioning scheme for accelerating convergence through some equivalent transformations of the eigenvalue problem. Numerical examples are given to illustrate the convergence properties and to demonstrate the competitiveness of the method.
Introduction
Iterative methods such as the Lanczos algorithm and the Arnoldi algorithm are widely used for
solving large matrix eigenvalue problems (see [21, 22]). Eective applications of these algorithms
typically use a shift-and-invert transformation, which is sometimes called preconditioning [22] and
requires solving a linear system of equations of the original size at each iteration of the process.
For truly large problems, solving the shift-and-invert equations by a direct method such as the LU
factorization is often infeasible or ine-cient. In those cases, one can employ an iterative method to
solve them approximately, resulting in two levels of iterations called inner-outer iterations. However,
methods like the Lanczos algorithm and the Arnoldi algorithm are very sensitive to perturbations
in the iterations and therefore require highly accurate solutions of these linear systems (see [8]).
Therefore, the inner-outer iterations may not oer an e-cient approach for these methods.
There has recently been great interest in other iterative methods that are also based on shift-
and-invert equations but tolerate low accuracy solutions. One simple example of such methods is the
inexact inverse iteration where the linear convergence property is preserved even when the inversion
is solved to very low accuracy (see [7, 13, 14, 27]). Several more sophisticated and competitive
methods have been developed that also possess such a property. They include the Jacobi-Davidson
Scientic Computing and Computational Mathematics Program, Department of Computer Science, Stanford
University, Stanford, CA 94305. E-mail : golub@sccm.stanford.edu. Research supported in part by National Science
Foundation Grant DMS-9403899.
y Department of Mathematics, University of Kentucky, Lexington, KY 40506-0027. E-mail: qye@ms.uky.edu. Part
of this Research was supported by NSERC of Canada while this author was with University of Manitoba.
method [5, 24, 25], truncated RQ iterations [26, 32] and others [14, 29, 30, 31]. One di-culty with
these methods is that it is not easy to determine to what accuracy the shift-and-invert equations
should be solved. On the other hand, there have been several works that aim at generalizing the
concept of preconditioning for linear systems to the eigenvalue problem [1, 2, 4, 11, 12, 17, 18, 15, 20,
28]. This is mostly done, however, by directly adapting a preconditioner used for inverting a certain
matrix into an eigenvalue iteration and in these situations, the role of preconditioners is usually
not clear, although some of them can be regarded as using inexact shift-and-invert [19]. Overall,
while all these new methods have been demonstrated to work successfully in some problems, there
is in general a lack of understanding of how and why they work. Furthermore, optimal eigenvalue
projection methods such as the Lanczos algorithm have mostly not been incorporated in these
developments.
In this paper, we shall present a variation of the Krylov subspace projection methods for computing
some extreme eigenvalues of the generalized eigenvalue problem
called the pencil problem for (A; B),
are symmetric matrices and B > 0. The method iteratively improves an approximate
eigenpair, each step of which uses either the Lanczos or the Arnoldi iteration to produce
a new approximation through the Rayleigh-Ritz projection on a Krylov subspace, resulting in a
form of inner-outer iterations. We shall present our theoretical and numerical ndings concerning
convergence properties of this method and derive bounds on asymptotic linear convergence rates.
Furthermore, we shall develop from the convergence analysis some equivalent transformations of
the eigenvalue problem to accelerate the convergence, which will be called preconditioning. In
particular, such transformations will be based on incomplete factorization and thus generalize the
preconditioning for linear systems. To the best of our knowledge, this is the rst preconditioning
scheme for the eigenvalue problem that is based on and can be justied by a convergence theory.
The paper is organized as follows. We present the basic algorithm in Section 2 and analyze
its convergence properties in section 3. We then give some numerical examples in Section 4 to
illustrate the convergence properties. We then present a preconditioned version of the algorithm
in Section 5, followed by some numerical examples on the preconditioning in Section 6. We nish
with some concluding remarks in Section 7.
Basic Inverse Free Krylov Subspace Method
In this section, we present our basic algorithm for nding the smallest eigenvalue and a corresponding
eigenvector (; x) of a pencil (A; B) where A; B are symmetric with B > 0. We note that
the method to be developed can be modied in a trivial way for nding the largest eigenvalue (or
simply considering ( A; B)).
Given an initial approximation we aim at improving it through the Rayleigh-Ritz
orthogonal projection on a certain subspace, i.e. by minimizing the Rayleigh quotient
on that subspace. Noting that the gradient of the Rayleigh quotient at x 0 is r
the well-known steepest descent method chooses a new approximate eigenvector x
by minimizing Clearly, this can be considered as the Rayleigh-Ritz projection method on
the subspace K 1  spanfx g. On the other hand, the inverse iteration constructs a
new approximation by x . If the inversion is solved inexactly by an iterative
solver (i.e. in an inexact inverse iteration [7]), then x 1 is indeed chosen from a Krylov subspace as
generated by A  0 B. Since x 1 is extracted from the Krylov subspace to solve the linear system,
it may not be a good choice for approximating the eigenvector.
We consider here a natural extension of these two approaches that nds a new approximate
eigenvector x 1 from the Krylov subspace
Km  spanfx
(for some xed m) by using the Rayleigh-Ritz projection method. The projection is carried out by
constructing a basis for Km and then forming and solving the projection problem for the pencil
B). Repeating the process, we arrive at the following iteration, which we call inverse free Krylov
method for (A; B).
Algorithm 1: Inverse Free Krylov Subspace Method.
Input m  1 and an initial approximation x 0 with
For
Construct a basis
Find the smallest eigenpair
End
In the algorithm, we apply the projection to the shifted pencil update the approximation
accordingly, which is theoretically equivalent to using the projection of (A; B) directly.
This formulation, however, may improve the stability while saving matrix-vector multiplications by
utilizing which need to be computed in the construction of the basis.
In constructing a basis z there are many possible choices and, theoretically,
they are all equivalent in that the new approximate eigenpair ( obtained will be the
same, which is dened by
(1)
However, numerically, we will consider a basis that is orthonormal on a certain inner product. Such
a basis of the Krylov subspace is typically constructed through an iterative method itself, which
will be called the inner iteration. The original iteration of Algorithm 1 will be called the outer
iteration. We shall discuss in the subsections later three methods for constructing the basis.
In our presentation of Algorithm 1, we have assumed for convenience that
that we construct a full basis z Generally, if
only. Then, the Rayleigh-Ritz projection is simply carried out by replacing Zm with
and (1) is still valid. Numerically, however, early termination at step p of the
inner iteration is not likely and a full basis is usually constructed even when p < m in theory,
but this causes no problem as the larger space spanned by more vectors would yield a better
approximation.
Given that the basic ingredients of Algorithm 1 are the projection and the Krylov subspaces,
it is not surprising that some similar methods have been considered before. In [10, 11], Knyazev
discussed and analyzed some very general theoretical methods and suggested several special cases,
among which is the use of Km and (1). Morgan and Scott's preconditioned Lanczos algorithm
[18] takes a similar iteration but uses the smallest Ritz value of the matrix A  k B rather than
that of the pencil to update the eigenvalue. With an m varied with each iteration,
it has a quadratic convergence property. We point out however that the quadratic convergence
is not a desirable property because it is achieved at the cost of increasingly larger m and it prevents
improvement of convergence by preconditioning (there is hardly any need to accelerate a
quadratic convergent algorithm). Our study will be somewhat of dierent nature in that we consider
accelerating convergence by changing certain conditions of the problem through equivalent
transformations (see section 5) as opposed to increasing m.
Also related to ours are methods based on inverting a shifted matrix A  k B or its projection,
which include the inverse iteration and the Jacobi-Davidson method [24]. When the inversion is
solved approximately by an iterative method, the solution is extracted from a Krylov subspace
generated by A  k B (or its projection). In these cases, it is chosen to satisfy the related linear
system. We note that the Jacobi-Davidson method also uses the Rayleigh-Ritz projection in the
outer iteration, the cost of which increases with the iteration. By xing the size of subspaces for
projection, the cost of Algorithm 1 is xed per outer iteration.
I, it is easy to see that Algorithm 1 is just the standard restarted Lanczos algorithm
for A. In this regard, our investigation is on the version with a xed m and on how m aects the
convergence. Furthermore, our development will lead to a preconditioning strategy that transforms
I) to the pencil problem (L 1 AL suitably chosen L), to which Algorithm
1 will be applied. This transformation to a more complicated problem may seem counter intuitive
but an important feature of Algorithm 1 is that the case I oers no advantage than a more
general B:
We now discuss in details the construction of a basis for Km in Algorithm 1.
2.1 Orthonormal basis by the Lanczos algorithm
One obvious choice of the basis of the Krylov subspace Km is the orthonormal one as constructed by
applying the Lanczos algorithm to C Simultaneously
with the Lanczos process, we produce a tridiagonal matrix
. The Lanczos process
requires m+1 matrix-vector multiplications by C Once the basis has been constructed,
we
multiplications by B. Note that
Here we state the Lanczos algorithm.
Algorithm 2: Orthonormal basis by Lanczos.
B, an approximate
For
End
With the orthonormal basis, Bm is in general a full matrix and we need to solve a generalized
eigenvalue problem for (Am ; Bm ). While in the exact arithmetic, it may not be valid in a
nite precision arithmetic for larger m when there could be severe loss of orthogonality among z i .
This can be corrected by either computing
explicitly or using reorthogonalization
[6] in the Lanczos process. We note that C k Zm has been computed in the Lanczos algorithm and
can be stored for forming Am .
2.2 B-orthonormal basis by the Arnoldi algorithm
We can also construct a B-orthonormal basis for Km by the modied Gram-Schmidt process in the
B-inner product, which is essentially the Arnoldi algorithm. The advantage of this approach is a
simpler projection problem with but it is at the cost of a longer recurrence. We also need
to compute
. We state the algorithm here.
Algorithm 3: B-orthonormal basis by Arnoldi.
Input
For
For
End
Each step of the Arnoldi algorithm requires 2 matrix-vector multiplications with one by C k and
one by B. In addition, we need to store Bz i from each iteration in order to save matrix-vector
multiplications, resulting in storage cost of m vectors. We note again that, for larger m, the B-
orthogonality among the columns of Zm may gradually be lost. This leads to deterioration of the
equation In that case, we need either reorthogonalization in the Arnoldi algorithm, or
explicit computations of
In comparing the two constructions, the computational costs associated with them are very
comparable. They both require 2(m multiplications. The Arnoldi recurrence
is more expensive in both
ops and storage than the Lanczos recurrence while it produces a more
compact projection matrix than the Lanczos algorithm. Clearly, these dierences are very minor
when m is not too large, which is the case of interest in practical implementations. In terms of
numerical stability of these two theoretically equivalent processes, our testing suggests that there
is very little dierence. However, for the preconditioned version of Algorithm 1 that we will discuss
in Section 5, the approach by the Arnoldi algorithm seems to have some advantage, see section 5.
2.3 C k -orthogonal basis by a variation of the Lanczos algorithm
It is also possible to construct Zm that is C k -orthogonal by a variation of the Lanczos algorithm
with a three term recurrence. Then the projection
will have a compact form,
leading to a computationally more eective approach than the previous two. However, it is less
stable owing to the indeniteness of C k . For the theoretical interest, we outline this variation of
the Lanczos algorithm for in the form of a full matrix tridiagonalization.
be the standard tridiagonalization of the Lanczos algorithm for C where T
is tridiagonal and Q is orthogonal with x k =kx k k as its rst column. For the sake of simplicity in
presentation, we assume here that  k is between the rst and the second eigenvalue, which implies
that C has exactly one negative eigenvalue. Noting that the (1; 1) entry of T is x T
be the block LDL T decomposition of T , where

I
and

Write
It is easy to check that
and
Now a Lanczos three term recurrence can be easily derived
from (2) to construct the columns of Z, which still form a basis for the Krylov subspace and is
essentially C-orthogonal. However, our tests show that this is numerically less stable. Therefore,
we shall not consider this further and omit a detailed algorithm here.
3 Convergence analysis.
In this section, we study convergence properties of Algorithm 1 that include a global convergence
result and a local one on the rate of linear convergence. In particular, we identify the factors
that aect the speed of the convergence so as to develop preconditioning strategy to improve the
convergence.
We rst prove that Algorithm 1 always converges to an eigenpair. For that, we need the
following proposition, the proof of which is straightforward.
Proposition 1 Let  1 be the smallest eigenvalue of (A; B) and ( k ; x k ) be the eigenpair approximation
obtained by Algorithm 1 at step k. Then
and
Theorem 1 Let ( k ; x k ) denote the eigenpair approximation obtained by Algorithm 1 at step k.
Then  k converges to some eigenvalue ^
of (A; B) and converges in
direction to a corresponding eigenvector).
Proof From Proposition 1, we obtain that  k is convergent. Since x k is bounded, there is a
convergent subsequence x n k . Let
x:
B)^x. Then it follows from (3) that,
Suppose now ^ r 6= 0. We consider the projection of (A; B) onto
rg by dening
r]:
Noting that f^x; ^ rg is orthogonal, we have ^

B)^r
is indenite. Thus the smallest eigenvalue of
B), denoted by ~ , is less than ^
, i.e.
~
Furthermore, at step k, dene r
Let ~  k+1 be the smallest eigenvalue of
B.
Hence by the continuity property of the eigenvalue, we have
On the other hand,  k+1 is the smallest eigenvalue of the projection of (A; B) on
which implies
Finally, combining the above together, we have obtained
which is a contradiction to (4). Therefore, ^
is an eigenvalue and
Now, to show
suppose there is a subsequence m k such that
> 0. From the subsequence m k , there is a subsequence n k for which x n k
is convergent. Hence by
virtue of the above proof, which is a contradiction. This completes the proof.
Next, we study the speed of convergence through a local analysis. In particular, we show that
k converges at least linearly.
be the smallest eigenvalue of (A; B), x be a corresponding unit eigenvector and
be the eigenpair approximation obtained by Algorithm 1 at step k. Let  1 be the smallest
eigenvalue of A  k B and u 1 be a corresponding unit eigenvector. Then
Asymptotically, if  k !  1 , we have
Proof First, from the denition, we have
Furthermore, A  1 I  k B  0 and A 0I is the smallest eigenpair
of is the smallest eigenpair of (A; B). Clearly A  k B is indenite and
hence  1  0: Now using Theorem 3 of Appendix, we have
which leads to the bound (5).
To prove the asymptotic expansion, let  1 (t) be the smallest eigenvalue of A tB. Then
. Using the analytic perturbation theory, we obtain  0
and hence
Choosing
from which the expansion follows.
We now present our main convergence result. We assume that  k is already between the rst
and the second smallest eigenvalues. Then by Theorem 1, it converges to the smallest eigenvalue.
Theorem 2 Let  1 <  2       n be the eigenvalues of (A; B) and ( k+1 ; x k+1 ) be the approximate
eigenpair obtained by Algorithm 1 from ( k ; x k ). Let  1   2       n be the eigenvalues
of A  k B and u 1 be a unit eigenvector corresponding to  1 . Assume  1 <  k <  2 . Then
kBk
where
and
with Pm denoting the set of all polynomials of degree not greater than m:
Proof First, write g. At step k of the algorithm, we
have
Let A  k be the eigenvalue decomposition of A  k B, where
orthogonal and g. Let q be the minimizing polynomial in  m with
and it follows from x T
and
Using Proposition 1, we have y T
and hence
On the other hand, we also have
and
where we note that q( 1
Thus
kBk
where we have used (8) and (9). Finally, combining (7), (10) and Lemma 1, we have
kBk
which leads to the theorem.
It is well known that  m in the theorem can be bounded in terms of  i as (see [16, Theorem
1.64] for example)


Then the speed of convergence depends on the distribution of the eigenvalues  i of A  k B but not
those of (A; B). This dierence is of fundamental importance as it allows acceleration of convergence
by equivalent transformations that change the eigenvalues of A  k B but leave those of (A; B)
unchanged (see the discussion on preconditioning in Section 5). On the other hand, the bound
shows accelerated convergence when m is increased. In this regard, our numerical tests suggests
that the convergence rate decreases very rapidly as m increases (see Section 4).
For the special case of just the steepest descent method for
A. It is easy to check in this case that
Using this in Theorem 2, we recover the classical convergence bound for the steepest descent
method [9, p.617]. We note that there is a stronger global convergence result in this case, i.e.  k is
guaranteed to converge to the smallest eigenvalue if the initial vector has a nontrivial component
in the smallest eigenvector (see [9, p.613]). There is no such result known for the case B 6= I.
Asymptotically we can also express the bound in terms of the eigenvalues of A  1 B instead
of  i which is dependent of k. We state it as the following corollary; but point out that the bound
of Theorem 2 is more informative.
n be the eigenvalues of A  1 B. Then, we have asymptot-
ically
p! 2m
The proof follows from combining (11) with  i
4 Numerical Examples - I
In this section, we present numerical examples to illustrate the convergence behavior of Algorithm
1. Here, we demonstrate the linear convergence property and the eect of m on the convergence
rate.
Example 1: Consider the Laplace eigenvalue problem with the Dirichlet boundary condition
on an L-shape region. A nite element discretization on a triangular mesh with 7585 interior
nodes (using PDE toolbox of MATLAB) leads to a pencil eigenvalue problem
apply Algorithm 1 to nd the smallest eigenvalue with a random initial vector and the stopping
criterion is set as kr k k=kr 0 We give the convergence history of
the residual kr k k for (from top down resp.) in Figure 1. We present in Figure
2 (a) the number of outer iterations required to achieve convergence for each m in the range of
correspondingly in Figure 2 (b), the total number of inner iterations.
We observe that the residual converges linearly with the rate decreased as m increases. Fur-
thermore, from Figure 2 (a), the number of outer iterations decreases very rapidly (quadratically
or even exponentially) as m increases and it almost reaches its stationery limit for m around 70.
Because of this peculiar property, we see from Figure 2 (b) that the total number of inner iterations
is near minimal for a large range of m (40 < m < 80 in this case).
Example 2: We consider a standard eigenvalue problem which A is a ve point
nite dierence discretization of the Laplace operator on the mesh of the unit
square. Again, we apply Algorithm 1 to nd the smallest eigenvalue with a random initial vector.

Figure

Example 1: convergence of r k for
outer iterations
2-norm
of
residuals
In this case, it is simply a restarted Lanczos algorithm and we shall consider its comparison with
the Lanczos algorithm without restart. In Figure 3, we present the convergence history of  k  1
where  k is the approximate eigenvalue obtained at each inner iteration. They are plotted in the
dot lines from top down for respectively. The corresponding plot for the Lanczos
algorithm (without restart) is given in the solid line.
We have also considered the number of outer iterations and the total number of inner iterations
as a function of m and observed the same behavior as in Example 1. We omit a similar gure
here. In particular, the nearly exponential decrease of the outer iteration count implies that the
convergence history with a moderate m (in this case even 16) will be very close to the
one with very large m (i.e. Lanczos without restart in the solid line ).
These examples conrm the linear convergence property of Algorithm 1. Furthermore, our
numerical testing has consistently shown that the number of outer iterations decreases nearly
exponentially as m increases. This implies that near optimal performance of the algorithm can
be achieved with a moderate m, which is very attractive in implementations. Unfortunately we
have not been able to explain this interesting behavior with our convergence results. Even for the
restarted Lanczos algorithm, it seems to be a phenomenon not observed before and can not be
explained by the convergence theory of the Lanczos algorithm either.
5 Preconditioning
In this section, we discuss how to accelerate the convergence of Algorithm 1 through some equivalent
transformations, which we call preconditioning, and we shall present the preconditioned version of
Algorithm 1.

Figure

2: Example 1: Outer and total inner iterations vs. m
outer
iterations
(a)
parameter m (inner iterations)
total
inner
iterations
(b)
From our convergence result (Theorem 2), the rate of convergence depends on the spectral
distribution of C i.e. the separation of  1 from the rest of eigenvalues  2 ;    ;  n of C k .
With an approximate eigenpair ( k ; x k ), we consider for some matrix L k the transformed pencil
which has the same eigenvalues as (A; B). Thus, applying one step of Algorithm 1 to
have the bound (6) of Theorem 2 with the rate of convergence  2
determined by ^
the eigenvalues of
We can now suitably choose L k to obtain a favorable distribution of ^  i and hence a smaller  m . We
shall call (13) a preconditioning transformation.
One preconditioning transformation can be constructed using the LDL T factorization of a
symmetric matrix [6]. For example, if
is the LDL T factorization with D k being a diagonal matrix of 1, choosing this L k results in ^
. Then, at the convergence stage with  1 <  k <  2 , we have ^
which implies
and thus by Theorem 2,

Figure

3: Example 2: Eigenvalue convergence history against each inner iteration for restarted
dot lines from top down and for Lanczos without restart in solid
line
total inner iterations
error
of
Ritz
value
We conclude that Algorithm 1, when applied to
k ) at step k using the exact LDL T factor-
ization, converges quadratically. This is even true with (i.e. the steepest descent method).
Similarly, in light of Corollary 1, if we use a constant L obtained from the LDL T factorization
(assuming  1 is known) with D k being a diagonal matrix of 0 and 1,
Algorithm 2 also converges quadratically.
What we have described above is the ideal situations of fast quadratic convergence property
achieved by using an exact LDL T factorization. In practice, we can use an incomplete LDL T
factorization A  k B  L k D k L T
(through incomplete LU factorization, see [23, Chapter 10]).
Then we will have a nonzero but small  m and hence fast linear convergence. Indeed, to be e-cient,
we can consider a constant L as obtained from an incomplete LDL T factorization of
where  0 is a su-ciently good approximation of  1 and apply Algorithm 1 to (13). Then, the
preconditioned algorithm converges linearly with the rate determined by the eigenvalues of
which has a better spectral distribution as long as ( 0  k )L 1 BL T is small relative to ^
We note that  0  k need not be very small if L 1 BL T is small (e.g. for the discretization of
dierential operators). It may work even when  0 <  1 , for which A  and the incomplete
LDL T factorization becomes incomplete Cholesky factorization. It is also possible to construct L
based on other factorization, such as approximate eigenvalue decomposition.
As in the preconditioned iterative methods for linear systems, the preconditioned iteration of
Algorithm 1 can be implemented implicitly, i.e. without explicitly forming the transformed problem
C k . We derive a preconditioned version of the algorithm in the rest of this section.
be the approximate eigenpair that has been obtained at step k for the pencil (A; B).
is the corresponding approximate eigenpair for the transformed
pencil (13). By applying one step of iteration to the transformed pencil, the new approximation is
obtained by constructing a basis ^ z
z m for the Krylov subspace
and form the projection problem
is the smallest eigenpair of the above projection problem,
then
v) is the new approximate eigenpair for (A; B).
z m ]:
Then, the new approximate eigenpair can be written as ( k+1 and the projection problem
is
Therefore, to complete the k-th iteration, we only need to construct
basis for the subspace L T
Km . The actual construction of z i depends on which method we use
and will be given in details in the subsections later.
Here, we summarize the preconditioned algorithm as follows.
Algorithm 4: Preconditioned Inverse Free Krylov Subspace Method.
Input m and an initial approximation x 0 with
For convergence
Construct a preconditioner L k ;
Construct a preconditioned basis
Km
Find the smallest eigenvalue  1 and a eigenvector v for (Am ; Bm );
End
Remark: As in the linear system case, the above algorithm takes the same form as the original
one except using a preconditioned search space L T
Km . In the following subsections, we discuss
the construction of a preconditioned basis by the Arnoldi algorithm and the Lanczos algorithm
corresponding to the construction in Sections 2.1 and 2.2. Our numerical testing suggests that the
Arnoldi algorithm might be more stable than the Lanczos algorithm in some cases.
5.1 Preconditioned basis by the Arnoldi method
In the Arnoldi method, we construct ^
B-orthonormal basis for
Km . Correspond-
ingly, z is a B-orthonormal basis for L T
Km . Starting from ^
B , the
recurrence for
z i is
where h
z T
z
and
with
is B-orthonormal, and h j;i above ensures this condition.
From this, we arrive at the following algorithm.
Algorithm 5: Preconditioned B-orthonormal basis by Arnoldi
Input and a preconditioner L k .
For
For
End
We see from the algorithm that only L T
k is needed in our construction. If we use  0 <  1
and is an incomplete Cholesky factor, i.e. A  0 B  LL T , then we can use any matrix
approximating
explicitly forming L k . For example, for dierential
operators, we can use the multigrid or domain decomposition preconditioners for A  0 B directly.
5.2 Preconditioned basis by the Lanczos method
In the Lanczos method, we construct ^ z
z m as an orthonormal basis for
Km . Then the corresponding
basis z
k . Starting from ^
the recurrence is
z
z
z T
z i and
. The resulting tridiagonal matrix T as constructed from 's
and 's satises
Zm . Thus, using ^
and
with
Clearly, the formulas for  i and  i+1 ensures
z j is M-orthonormal. Thus, we have the alternative formulas
From this, we can derive a recurrence to construct the basis. We note that this construction
normalizes z i in the M-norm. In practice, M could be nearly singular. Therefore, it is more
appropriate to normalize it in the 2-norm. The following algorithm is one of several possible
formulations here.
Algorithm Preconditioned basis by Lanczos
Input and a preconditioner L k .
For
End
6 Numerical Examples - II
In this section, we present some numerical examples to demonstrate the eectiveness and competitiveness
of the preconditioned inverse free Krylov subspace method (Algorithm 4).
Example 3: A and B are the same as in Example 1. We apply Algorithm 4 to nd the smallest
eigenvalue (closest to  We use a constant L k as obtained by the threshold incomplete
LDL T factorization of A  A with the drop tolerance 10 2 . We compare our algorithm
with the Jacobi-Davidson algorithm that uses the same number of inner iterations (m) and with
the same kind of preconditioner. We give in Figure 4 the convergence history of the residual kr k k
of Algorithm 4 in solid lines and that of the Jacobi-Davidson algorithm in dot lines from top down
respectively. In Figure 5, we also present the number of outer iterations and the
total number of inner iterations required to reduce kr k k=kr 0 k to 10 7 for each m in (+) mark for
Algorithm 4 and in (o) mark for the Jacobi-Davidson algorithm.
Comparing it with Example 1, the results clearly demonstrate the acceleration eect of pre-conditioning
by signicantly reducing the number of outer iterations (Fig. 4 and Fig. 5(a)).
Furthermore, the values of m at which the total number of inner iterations is near its minimum are
signicantly smaller with preconditioning (around Fig. 5). Although J-D algorithm
has smaller number of total inner iteration for very small m, the corresponding outer iteration
count is larger, which increases its cost.
We also considered for this example the ideal preconditioning with L k chosen as the exact LDL T
factorization of C k . In this case, we use an initial vector with kAx 0  so that  0
is su-ciently close to  1 . We present the residual convergence history in Figure 6 for
steepest descent method), 2 and 4. The result conrms the quadratic convergence property for all

Figure

4: Example 3 residual convergence history for lines - Algorithm 4;
dot lines - Jacobi-Davidson)
outer iterations
2-norm
of
residuals
m. We have also tested the case that uses L L from the exact factorization of A  1 B and in
this case it converges in just one iteration, conrming Corollary 1.
The next example is for the standard eigenvalue problem and the preconditioned
algorithm implicitly transforms it to a pencil problem.
Example 4: A is the same matrix as in Example 2 We use a constant L k as obtained
by the incomplete LDL T decompositions of A  A with no ll-in. We compare it with the
Jacobi-Davidson algorithm with the same kind of preconditioner. We also consider the shift-and-
spectral transformed Lanczos algorithm (i.e. applying the Lanczos to A 1 ).
We give in Figure 7 the convergence history of the residual kr k k of Algorithm 4 in solid lines
from top down for and that of the Jacobi-Davidson algorithm in dot lines for
with the corresponding marks. The residual for the spectral transformed Lanczos
is given in dash-dot (with +) line. Figure 8 is the number of outer iterations and the total number
of inner iterations vs. m.
Again, preconditioning signicantly accelerates convergence and our result compares very favorably
with the Jacobi-Davidson method. An interesting point here is that Algorithm 4 with
based on incomplete factorization outperforms the shift-and-invert Lanczos algorithm. Although
we do not suggest this is the case in general, it does underline the eectiveness of the preconditioned
algorithm.

Figure

5: Example 3 Outer and total inner iterations vs. m (+ - Algorithm 4;
outer
iterations
(a)
total
inner
iterations
(b)
7 Concluding Remarks
We have presented an inverse free Krylov subspace method that is based on the classical Krylov
subspace projection methods but incorporates preconditioning for e-cient implementations. A
convergence theory has been developed for the method and our preliminary tests of the preconditioned
version demonstrate its competitiveness. Comparing with the existing methods, it has a
relatively well understood theory and simple numerical behavior. We point out that the algorithm
has a xed cost per outer iteration, which makes it easy to implement.
For the future work, we will consider generalizations in three directions, namely, an e-cient
block version for computing several eigenvalues simultaneously, a strategy to compute interior
eigenvalues and an algorithm for the general nonsymmetric problem.
A

Appendix

. Perturbation Bounds for Generalized Eigenvalue
Problems
We present a perturbation theorem that is used in the proof of Lemma 1 but might be of general
interest as well. In the following, A; B; and E are all symmetric.
Theorem 3 Let  1 be the smallest eigenvalue of (A; B) with x a corresponding unit eigenvector
and let  1 be the smallest eigenvalue of with u a corresponding unit eigenvector, where
min (E)

Figure

convergence with ideal preconditioning
where  min (E) and  max (E) denote the smallest and the largest eigenvalues of E respectively.
Proof Using the minimax characterization, we have
z 6=0
z T
z T Bz
Similarly,
min (E)
which completes the proof.
We note that 1=x T Bx (or 1=u T Bu ) is the Wilkinson condition number for  1 (or  1 ). These
bounds therefore agree with the rst order analytic expansion and will be sharper than traditional
bounds based on kBk.

Figure

7: Example 4 residual convergence history for lines - Algorithm
4; dot lines - Jacobi-Davidson, dash-dot - shift-and-invert Lanczos)



--R

A subspace preconditioning algorithm for eigenvec- tor/eigenvalue computation
The Davidson method
Applied Numerical Linear Algebra
Minimization of the computational labor in determining the

Matrix Computations
Inexact inverse iterations for the eigenvalue problems
Large sparse symmetric eigenvalue problems with homogeneous linear constraints: the Lanczos process with inner-outer iterations
Functional Analysis in Normed Spaces
Convergence rate estimates for iterative methods for a mesh symmetric eigenvalue problem.
Preconditioned Eigensolvers - An oxymoron? Elec
Toward the optimal preconditioned eigensolver: Locally optimal block preconditioned conjugate gradient method
An inexact inverse iteration for large sparse eigenvalue problems
The inexact rational Krylov sequence method
The restarted Arnoldi method applied to iterative linear solvers for the computation of rightmost eigenvalues
Computer Solution of Large Linear Systems
Generalizations of Davidson's method for computing eigenvalues of sparse symmetric matrices
Preconditioning the Lanczos algorithm for sparse symmetric eigenvalue problems
A geometric theory for preconditioned inverse iteration
a new method for the generalized eigenvalue problem and convergence estimate Preprint
The Symmetric Eigenvalue Problem
Numerical Methods for Large Eigenvalue Problems
Iterative Methods for Sparse Linear Systems
A Jacobi-Davidson iteration method for linear eigenvalue problems

A truncated RQ iteration for large scale eigenvalue calculations

Robust preconditioning of large sparse symmetric eigenvalue problems
Restarting techniques for the (Jacobi-)Davidson symmetric eigenvalue method Elec
Dynamic thick restarting of the Davidson and implicitly restarted Arnoldi methods
Inexact Newton preconditioning techniques for large symmetric eigenvalue problems Elec.
Convergence analysis of an inexact truncated RQ iterations Elec.
--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
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
preconditioning;eigenvalue problems;krylov subspace