| Home > Research > Ifsdocs > ENSEMBLE > |
IFS documentation front page
Table of contents
Chapter 1. Theory
Chapter 2. Computational
details
Previous Section
APPENDIX A Eigenvalue algorithms
A.1 The Lanczos algorithm
Algorithms based on Lanczos theory are very useful to solve an eigenvalue problem when only a few of the extreme eigenvectors are needed. It can be applied to large and sparse problems. The algorithm does not access directly the matrix elements of the operator that defines the problem, but it gives an estimate of the eigenvectors through successive application of the operator.
Consider the eigenvalue problem
|
(A.1) |
where the matrix
is 
dimensional, and symmetric. Without loss of generality, we can also
suppose that it is real.If
is a real, symmetric matrix, then there exists an orthogonal
real matrix 
such that
 ,
|
(A.2) |
where
is a diagonal matrix, and 
denotes the transpose of 
(Schur decomposition theorem).The Lanczos algorithm does not directly compute the diagonal matrix
, but it first computes a partial transformation of the matrix
using a tridiagonal matrix 
 ,
|
(A.3) |
with
 ,
|
(A.4) |
and with
|
(A.5) |
where the vectors
are column vectors, and where the number of
iterations 
is much smaller than the dimensionality
of the problem, 
. Then, the Lanczos algorithm finds the diagonal decomposition of
 .
|
(A.6) |
The elements of the diagonal matrix
are an estimate of the eigenvalues of 
, and an estimate of the eigenvectors are given by 
, with
 ,
|
(A.7) |
The actual computation is performed by writing Eq. (A.3) as
 .
|
(A.8) |
Equating columns of Eq. (A.8), it follows that
|
(A.9) |
for
. The orthogonality of the vectors 
implies that
 .
|
(A.10) |
Moreover, if
|
(A.11) |
is non-zero, then
 ,
|
(A.12) |
where
. An iterative application of these equations, with a
randomly chosen starting vector 
, defines the Lanczos iterative procedure. The total number of iterations
determines the accuracy of the computation.
As this number increases, more eigenvalues/eigenvectors can be separated
from the others, independently from the choice of the starting vector 
. This separation starts from the boundaries of the eigenvalue spectrum.
The accuracy of the eigenvectors is less than the accuracy of the singular
values, say to order 
when the precision of the singular values is of order 
.The reader is referred to Golub and van Loan (1983) for a theoretical description of the Lanczos algorithm. The Lanczos code is available in NAG issue 17.
A.2 The Jacobi-Davidson algorithm
The generalized eigenproblem
|
(A.13) |
is usually handled by bringing it back to a standard eigenproblem
|
(A.14) |
The matrix
is in general nonsymmetric, even if both 
and 
are symmetric. However, if 
is symmetric and positive definite, the 
-inner product is well defined. The matrix 
is symmetric in this inner product if 
is symmetric:
 .
|
(A.15) |
The proposed method to solve (A.14) constructs a set of basis vectors
of a search space 
, c.f. the Lanczos method. The approximate eigenvectors are linear
combinations of the vectors 
. The classical and most natural choice for the search space 
, for instance utilized in the Lanczos method, is the
so-called Krylov subspace, the space spanned by the vectors
|
(A.16) |
This
-dimensional subspace is denoted by 
. The vector 
is a starting vector that has to be
chosen. The Krylov subspace is well suited for computing dominant eigenpairs
since the vector 
points more and more in the direction of 
the dominant eigenvector of for increasing 
.Given a search space
, the approximate eigenpair 
of (A.15) is a linear combination of the
basis vectors of 
:
 .
|
(A.17) |
A suitable criterion for finding an optimal pair is the Galerkin condition that the residual
|
(A.18) |
is
-orthogonal to the search space to 
. Hence
|
(A.19) |
and consequently, using (A.19),
|
(A.20) |
Note that the resulting eigenproblem is of the dimension of the search space, which is generally much smaller than of the original problem. The basis vectors are usually orthogonalized so that
Approximate eigenpairs that adhere to the Galerkin condition are
called Ritz pairs.It can be shown that the residuals
form a 
-orthogonal basis for 
when the approximate eigenpairs are computed according to (1.37).A.3 Solution method
As was stated before, the natural search space for the generalized eigenvalue problem is the Krylov subspace
. This basis can be generated by expanding the basis by new residual
vectors. The problem in the construction of a basis for this space is that
operations with 
are needed. Since in our application the inverse of 
is not known explicitly, its action is approximated by the
Conjugate Residual method (CR). To compute the vector
|
(A.21) |
one iteratively solves the system
|
(A.22) |
Iterative solution methods require, apart from vector operations, only multiplications with
. The vector 
can in principle be determined to high accuracy. This, however, may
require many multiplications with 
and hence may be very expensive. Therefore, the action of 
is approximated to low accuracy, by performing only a few steps
with an iterative solution method. The number of iterations is controlled
by NINNER. The subspace generated in this way is not a Krylov subspace and
the basis vectors are not the residuals (A.19)
but only approximations to it. As a consequence they are not perfectly 
-orthogonal. This has to be done explicitly.The complete algorithm can be summarized as follows.
• Choose a starting vector  |
• Compute  , B-normalize  |
| • Repeat the following 7 steps NITERL times |
1) compute  |
2) Solve small eigenproblem  |
3) Select Ritz value  and  |
4) Compute Ritz vector  and residual  |
5) Compute approximately  with the CR-method |
6) B-orthonormalize new  against  |
7) Expand  with the resulting vector |
The matrix
contains the basis for the search space, the vector
contains the new basis vector and 
is an approximate eigenpair. In step 3 the pair with the largest
is selected because the dominant part of the
spectrum is only of interest. However, if the eigenpair approximation reaches
a certain accuracy 
(XKAPA in namelist NAMLCZ), i.e., 
, a smaller 
is selected. In step 5 a few CR steps
are performed to approximate the action of 
. In step 6 the Modified Gram-Schmidt procedure is used for reasons
of numerical stability. For a more detailed descripton of the Jacobi-Davidson
algorithm the reader is referred to
Sleijpen and van der Vorst (1996)Next Section
Previous Section
 |
10.04.2002 |
 |
© ECMWF |
