An Implicitly Restarted Symplectic Lanczos Method for the Symplectic Eigenvalue Problem.
An implicitly restarted symplectic Lanczos method for the symplectic eigenvalue problem is presented. The Lanczos vectors are constructed to form a symplectic basis. The inherent numerical difficulties of the symplectic Lanczos method are addressed by inexpensive implicit restarts. The method is used to compute some eigenvalues and eigenvectors of large and sparse symplectic operators.
Introduction
. We consider the numerical solution of the real symplectic
eigenvalue problem
where M 2 IR 2n\Theta2n is large and possibly sparse. A matrix M is called symplectic iff
or equivalently, M T
and I n is the n \Theta n identity matrix. The symplectic matrices form a group under
multiplication. The eigenvalues of symplectic matrices occur in reciprocal pairs: If -
is an eigenvalue of M with right eigenvector x, then - \Gamma1 is an eigenvalue of M with
left eigenvector (Jx) T . The computation of eigenvalues and eigenvectors of such matrices
is an important task in applications like the discrete linear-quadratic regulator
problem, discrete Kalman filtering, or the solution of discrete-time algebraic Riccati
equations. See, e.g., [21, 22, 28] for applications and further references. Symplectic
matrices also occur when solving linear Hamiltonian difference systems [6].
In order to develop fast, efficient, and reliable methods, the symplectic structure
of the problem should be preserved and exploited. Then important properties of symplectic
matrices (e.g., eigenvalues occurring in reciprocal pairs) will be preserved and
not destroyed by rounding errors. Different structure-preserving methods for solving
have been proposed. In [25], Lin introduces the S
which can be used to compute the eigenvalues of a symplectic matrix by a structure-preserving
method similar to Van Loan's square-reduced method for the Hamiltonian
eigenvalue problem [38]. Flaschka, Mehrmann, and Zywietz show in [14] how to construct
structure-preserving methods based on the SR method [10, 11, 26]. Patel
[34, 33] and Mehrmann [27] developed structure-preserving algorithms for the symplectic
generalized eigenproblem
submitted in July 1998
y Universit?t Bremen, Fachbereich 3 - Mathematik und Informatik, Zentrum f?r Technomathe-
matik, 28357 Bremen, FRG. E-mail: benner@math.uni-bremen.de
z corresponding author, Universit?t Bremen, Fachbereich 3 - Mathematik und Informatik, Zentrum
f?r Technomathematik, 28357 Bremen, FRG. E-mail: heike@math.uni-bremen.de
Benner and Fa-bender
Recently, Banse and Bunse-Gerstner [2, 3] presented a new condensed form for
symplectic matrices. The 2n \Theta 2n condensed matrix is symplectic, contains
nonzero entries, and is determined by 4n \Gamma 1 parameters. This condensed form, called
symplectic butterfly form, can be depicted as a symplectic matrix of the following
@ @@ @7 5
Once the reduction of a symplectic matrix to butterfly form is achieved, the SR
algorithm [10, 11, 26] is a suitable tool for computing the eigenvalues/eigenvectors of a
symplectic matrix. The SR algorithm preserves the butterfly form in its iterations and
can be rewritten in a parameterized form that works with the 4n\Gamma1 parameters instead
of the (2n) 2 matrix elements in each iteration. Hence, the symplectic structure, which
will be destroyed in the numerical process due to roundoff errors, can be restored in
each iteration for this condensed form. An analysis of the butterfly SR algorithm can
be found in [2, 4, 5].
In [2, 3] an elimination process for computing the butterfly form of a symplectic
matrix is given which uses elementary unitary symplectic transformations as well as
non-unitary symplectic transformations. Unfortunately, this approach is not suitable
when dealing with large and sparse symplectic matrices as an elimination process can
not make full use of the sparsity. Hence, symplectic Lanczos methods which create
the symplectic butterfly form if no breakdown occurs are derived in [2, 4]. Given
and a symplectic matrix M 2 IR 2n\Theta2n , these Lanczos algorithms produce
a matrix S which satisfies a recursion
of the form
MS
is a butterfly matrix of order 2k \Theta 2k, and the columns of S 2n;2k are
orthogonal with respect to the indefinite inner product defined by J (1.3). The latter
property will be called J-orthogonality throughout this paper. The residual r k+1
depends on v k+1 and w Such a symplectic Lanczos
method will suffer from the well-known numerical difficulties inherent to any Lanczos
method for unsymmetric matrices. In [2], a symplectic look-ahead Lanczos algorithm
is presented which overcomes breakdown by giving up the strict butterfly form. Un-
fortunately, so far there do not exist eigenvalue methods that can make use of that
special reduced form. Standard eigenvalue methods as QR or SR algorithms have to
be employed resulting in a full symplectic matrix after only a few iteration steps.
A different approach to deal with the numerical difficulties of the Lanczos process
is to modify the starting vectors by an implicitly restarted Lanczos process (see
the fundamental work in [9, 35]); for the unsymmetric eigenproblem the implicitly
restarted Arnoldi method has been implemented very successfully, see [24]). The
problems are addressed by fixing the number of steps in the Lanczos process at a prescribed
value k which depends upon the required number of approximate eigenvalues.
J-orthogonality of the k Lanczos vectors is secured by re-J-orthogonalizing these
vectors when necessary. The purpose of the implicit restart is to determine initial
vectors such that the associated residual vectors are tiny. Given (1.4), an implicit
Lanczos restart computes the Lanczos factorization
An implicitly restarted symplectic Lanczos method 3
which corresponds to the starting vector
(where p(M) 2 IR 2n\Theta2n is a polynomial) without having to explicitly restart the
Lanczos process with the vector - v 1 . Such an implicit restarting mechanism is derived
here analogous to the technique introduced in [4, 18, 35].
Section 2 reviews the symplectic butterfly form and some of its properties that will
be helpful for analyzing the symplectic Lanczos method which reduces a symplectic
matrix to butterfly form. This symplectic Lanczos method is presented in Section 3.
Further, that section is concerned with finding conditions for the symplectic Lanczos
method terminating prematurely such that an invariant subspace associated with
certain desired eigenvalues is obtained. We will also consider the important question
of determining stopping criteria. The implicitly restarted symplectic Lanczos method
itself is derived in Section 4. Numerical properties of the proposed algorithm are
discussed in Section 5. In Section 6, we present some preliminary numerical examples.
2. The Symplectic Butterfly Form. A symplectic matrix
is called a butterfly matrix if B 11 and B 21 are diagonal, and B 12 and B 22 are tridiag-
onal. Banse and Bunse-Gerstner [2, 3] showed that for every symplectic matrix M ,
there exist numerous symplectic matrices S such that MS is a symplectic
butterfly matrix. In [2], an elimination process for computing the butterfly form of a
symplectic matrix is presented (see also [4]).
In [4], an unreduced butterfly matrix is introduced in which the lower right tridiagonal
matrix is unreduced, that is, the subdiagonal elements of B 22 are nonzero. Using
the definition of a symplectic matrix, one easily verifies that if B is an unreduced
butterfly matrix, then B 21 is nonsingular. This allows the decomposition of B into
two simpler symplectic matrices:
I T
I @@ @
is tridiagonal and symmetric. Hence parameters that determine
the symplectic matrix can be read off directly. The unreduced butterfly matrices
play a role analogous to that of unreduced Hessenberg matrices in the standard QR
theory [2, 4, 5].
We will frequently make use of the decomposition (2.1) and will denote it by
a
a
(2.
4 Benner and Fa-bender
. d2
. dn
. b2 d2
.
. a2d2
.
an
Remark 2.1. (See [4].)
a) Any unreduced butterfly matrix is similar to an unreduced butterfly matrix
with
b) We will have deflation if d j. Then the eigenproblem can be
split into two smaller ones with unreduced symplectic butterfly matrices.
Eigenvalues and eigenvectors of symplectic butterfly matrices can be computed
efficiently by the SR algorithm [7], which is a QR like algorithm in which the QR decomposition
is replaced by the SR decomposition. Almost every matrix A 2 IR 2n\Theta2n
can be decomposed into a product A = SR where S is symplectic and R is J -
triangular, that is
R
where all submatrices R ij 2 IR n\Thetan are upper triangular, and R 21 is strictly upper
triangular [12]. In the following a matrix D 2 IR 2n\Theta2n will be called trivial if it is
both symplectic and J-triangular. D is trivial if and only if it has the form
where C and F are diagonal matrices, C nonsingular.
If the SR decomposition A = SR exists, then other SR decompositions of A can
be built from it by passing trivial factors back and forth between S and R. That
is, if D is a trivial matrix, ~
R is another SR
decomposition of A. If A is nonsingular, then this is the only way to create other SR
decompositions. In other words, the SR decomposition is unique up to trivial factors.
The SR algorithm is an iterative algorithm that performs an SR decomposition
at each iteration. If B is the current iterate, then a spectral transformation function
An implicitly restarted symplectic Lanczos method 5
q is chosen (such that q(B) 2 IR 2n\Theta2n ) and the SR decomposition of q(B) is formed,
if possible:
Then the symplectic factor S is used to perform a similarity transformation on B to
yield the next iterate, which we will call b
B:
If is an unreduced symplectic butterfly matrix, then so is
B in (2.5) [2, 3]. If rank is an unreduced symplectic
butterfly matrix, then b
in (2.5) is of the form (see [4])
@ @@ @
-z-z-z-z-
where
is a symplectic butterfly matrix and
ffl the eigenvalues of
are just the - shifts that are eigenvalues of B.
An algorithm for explicitly computing S and R is presented in [8]. As with explicit
QR steps, the expense of explicit SR steps comes from the fact that q(B) has to be
computed explicitly. A preferred alternative is the implicit SR step, an analogue to
the Francis QR step [15, 17, 20]. As the implicit SR step is analogous to the implicit
QR step, this technique will not be discussed here (see [4, 5] for details).
A natural way to choose the spectral transformation function q is to choose a
polynomial
these choices make use of the symmetries
of the spectrum of symplectic matrices. But, as explained in [5], a better choice is a
Laurent polynomial to drive the SR step. For example, instead of p 4 (-) we will use
2:
This reduces the size of the bulges that are introduced, thereby decreasing he number
of computations required per iteration. Moreover, the use of Laurent polynomials improves
the convergence and stability properties of the algorithm by effectively treating
each reciprocal pair of eigenvalues as a unit. Using a generalized Rayleigh-quotient
strategy, the butterfly SR algorithm is typically cubic convergent [5].
The right eigenvectors of unreduced butterfly matrices have the following property
which will be helpful when analyzing the symplectic Lanczos method introduced in
the next section.
6 Benner and Fa-bender
Lemma 2.2. Suppose that B 2 IR 2n\Theta2n is an unreduced butterfly matrix as in
(2.4). If
In order to proof this lemma we need the following definition. Let Pn be the
permutation matrix
If the dimension of Pn is clear from the context, we leave off the subscript.
Proof. The proof is by induction on the size of B. The entries of the eigenvector
x will be denoted by x
Suppose that 2. The second and fourth row of
a
Since B is unreduced, we know that a 2 6= 0 and d 2 6= 0. If x then from (2.9) we
obtain
while (2.8) gives b 2 Using (2.10) we obtain x
The third row of
a
Using since B is unreduced, we obtain x
which contradicts the assumption x 6= 0.
Assume that the lemma is true for matrices of order 2(n\Gamma1). Let B 2n;2n 2 IR 2n\Theta2n
be an unreduced butterfly matrix. For simplicity we will consider the permuted
equation B 2n;2n
Partition
2n\Gamma2 an an c n5 ;
~
x
~
is an unreduced butterfly matrix and y 2
This implies
since an 6= 0 as B 2n;2n is unreduced. Further we have
Hence, using (2.11) we get ~ x
-y. Using
~
further obtain from (2.11) y This is a contradiction,
because by induction hypothesis e T
Remark 2.3. If y be the right eigenvector
of B to - y, then (Jy)
that e T
2n y 6= 0, hence the nth component of the left eigenvector of B corresponding to
- is 6= 0.
An implicitly restarted symplectic Lanczos method 7
3. A Symplectic Lanczos Method for Symplectic Matrices. In this sec-
tion, we review the symplectic Lanczos method to compute the butterfly form (2.4)
for a symplectic matrix M derived in [4]. The usual unsymmetric Lanczos algorithm
generates two sequences of vectors. Due to the symplectic structure of M it is easily
seen that one of the two sequences can be eliminated here and thus work and storage
can essentially be halved. (This property is valid for a broader class of matrices, see
[16].) Further, this section is concerned with finding conditions for the symplectic
Lanczos method terminating prematurely such that an invariant subspace associated
with certain desired eigenvalues is obtained. Finally we will consider the important
question of determining stopping criteria.
In order to simplify the notation we use in the following permuted versions of M
and B as in the previous section. Let
with the permutation matrix P as in (2.7).
3.1. The symplectic Lanczos factorization. We want to compute a symplectic
matrix S such that S transforms the symplectic matrix M to a symplectic
butterfly matrix B; in the permuted version MS = SB yields
Equivalently, as
, we can consider
where
a
a
. dn 0
.
The structure preserving Lanczos method generates a sequence of permuted symplectic
matrices
satisfying
8 Benner and Fa-bender
T is a permuted 2k \Theta 2k symplectic butterfly matrix.
The vector r k+1 := d k+1 (b k+1 is the residual vector and is JP -
orthogonal to the columns of S 2n;2k
P , the Lanczos vectors. The matrix B 2k;2k
P is the
-orthogonal projection of MP onto the range of S 2n;2k
Here J 2k;2k
P denotes a permuted 2k \Theta 2k matrix J of the form (1.3). Equation (3.5)
defines a length 2k Lanczos factorization of MP . If the residual vector r k+1 is the zero
vector, then equation (3.5) is called a truncated Lanczos factorization when k ! n.
Note that r n+1 must vanish since (S 2n;2n
and the columns of S 2n;2n
form a JP -orthogonal basis for IR 2n . In this case the symplectic Lanczos method
computes a reduction to permuted butterfly form.
The symplectic Lanczos factorization is, up to multiplication by a trivial matrix,
specified by the starting vector v 1 (see [4, Theorem 4.1]).
wn ]. For a given v 1 , a Lanczos method constructs
the matrix SP columnwise from the equations
From this we obtain the algorithm given in Table 3.1 (for a more detailed discussion
see [4]).
method
Choose an initial vector e
do
(update of wm )
set
e
am
e
wm
(computation of c m )
(update of v m+1 )
e

Table
Symplectic Lanczos Method
Remark 3.1. Using the derived formulae for w k+1 , the residual term r
can be expressed as
An implicitly restarted symplectic Lanczos method 9
There is still some freedom in the choice of the parameters that occur in this
algorithm. Essentially, the parameters b m can be chosen freely. Here we set b
Likewise a different choice of the parameters am ; dm is possible.
Note that M \Gamma1
since M is symplectic. Thus M \Gamma1
just a
matrix-vector-product with the transpose of MP . Hence, only one matrix-vector
product is required for each computed Lanczos vector wm or v m . Thus an efficient
implementation of this algorithm requires 6n nz is the
number of nonzero elements in MP and 2k is the number of Lanczos vectors computed
(that is, the loop is executed k times). The algorithm as given in Table 3.1 computes
an odd number of Lanczos vectors, for a practical implementation one has to omit
the computation of the last vector v k+1 (or one has to compute an additional vector
w
In the symplectic Lanczos method as given above we have to divide by parameters
that may be zero or close to zero. If such a case occurs for the normalization parameter
dm+1 , the corresponding vector e v m+1 is zero or close to the zero vector. In this case, a
(good approximation to a) JP -orthogonal invariant subspace of MP or equivalently, a
symplectic invariant subspace of M is detected. By redefining e v m+1 to be any vector
satisfying
m, the algorithm can be continued. The resulting butterfly matrix is no
longer unreduced; the eigenproblem decouples into two smaller subproblems. In case
e
wm is zero (or close to zero), an invariant subspace of MP with dimension
is found (or a good approximation to such a subspace). In this case the parameter
am will be zero (or close to zero). From Table 3.1 we further obtain that in this case
is a real eigenvalue of MP (and hence of M) with corresponding
Due to the symmetry of the spectrum of M , we also have
that 1=b m is an eigenvalue of M . Computing an eigenvector y of MP corresponding
to 1=b m , we can try to augment the (2m \Gamma 1)-dimensional invariant subspace to
an MP -invariant subspace of even dimension. If this is possible, the space can be
made JP -orthogonal by JP -orthogonalizing y against f and
normalizing such that y T JP
Thus if either v m+1 or wm+1 vanishes, the breakdown is benign. If v m+1 6= 0
and wm+1 6= 0 but then the breakdown is serious. No reduction of the
symplectic matrix to a symplectic butterfly matrix with v 1 as first column of the
transformation matrix exists.
A convergence analysis for the symplectic Lanczos algorithm analogous to the one
for the unsymmetric Lanczos algorithm presented by Ye [39] can be given. Moreover,
an error analysis of the symplectic Lanczos algorithm in finite-precision arithmetic
analogous to the analysis for the unsymmetric Lanczos algorithm presented by Bai
[1] can also be derived. These results will be presented in [13]. As to be expected, the
computed Lanczos vectors loose J(JP )-orthogonality when some Ritz values begin to
converge.
3.2. Truncated symplectic Lanczos factorizations. This section is concerned
with finding conditions for the symplectic Lanczos method terminating prema-
turely. This is a welcome event since in this case we have found an invariant symplectic
(Following [17], we define each floating point arithmetic operation together with the associated
integer indexing as a flop.)
Benner and Fa-bender
subspace S 2n;2k and the eigenvalues of B 2k;2k are a subset of those of M . We will first
discuss the conditions under which the residual vector of the symplectic Lanczos factorization
will vanish at some step k. Then we will show how the residual vector and
the starting vector are related. Finally a result indicating when a particular starting
vector generates an exact truncated factorization is given.
First the conditions under which the residual vector of the symplectic Lanczos
factorization will vanish at some step k will be discussed. From the derivation of the
algorithm it is immediately clear that if no breakdown occurs, then
where K(X; v; vg. Further it is
easy to see that
If dim K(MP
Hence, there exist real scalars ff such that
Using the definition of a k+1 as given in Table 3.1 and the above expression we obtain
because of J-orthogonality,
a
0:
As e
This implies that an invariant subspace of MP with dimension 2k
If dim K(MP
g. Hence
a
for properly chosen ff and from the algorithm in Table 3.1
An implicitly restarted symplectic Lanczos method 11
Therefore e v This implies that the residual vector of the
symplectic Lanczos factorization will vanish at the first step k such that the dimension
of K(M; is equal to 2k and hence is guaranteed to vanish for some k - n.
Next we will discuss the relation between the residual term and the starting vector.
If dim K(M;
and Cn is a generalized companion matrix of the form
. 1
(see [2, proof of Satz 3.6]). Thus,
Define the residual in (3.7) by
Note that
where
We will now show that f k+1 is up to scaling the residual of the length 2k symplectic
Lanczos iteration with starting vector v 1 . Together with (3.9) this reveals the relation
between residual and starting vectors. Since det (C
J-orthogonal columns
(that is, (S 2n;2k ) T JnS is a J-triangular matrix. Then
. The diagonal elements of R are nonzero if and only if the columns of
are linear independent. Choosing
12 Benner and Fa-bender
assures that (\GammaJ k (S 2n;2k multiplying (3.7) from the right by
is an unreduced butterfly matrix (see [2, proof of Satz 3.6])
with the same characteristic polynomial as C k . Equation (3.10) is a valid symplectic
Lanczos recursion with starting vector v residual vector f k+1 =r 2k;2k .
By (3.9) and due to the essential uniqueness of the symplectic Lanczos recursion any
symplectic Lanczos recursion with starting vector v 1 yields a residual vector that can
be expressed as a polynomial in M times the starting vector v 1 .
Remark 3.2. From (3.8) it follows that if dim K(M; then we
can choose c 1 ; :::; c 2k such that f This shows that if the Krylov subspace
forms an 2k-dimensional M-invariant subspace, the residual of the
symplectic Lanczos recursion will be zero after k Lanczos steps such that the columns
of S 2n;2k span a symplectic basis for the subspace K(M; 1).
The final result of this section will give necessary and sufficient conditions for a
particular starting vector to generate an exact truncated factorization in a similar
way as stated for the Arnoldi method in [35]. This is desirable since then the columns
of S 2n;2k form a basis for an invariant symplectic subspace of M and the eigenvalues
of B 2k;2k are a subset of those of M . Here, -
will denote the Lanczos vectors
after permuting them back, i.e., - v
Theorem 3.3. Let MS 2n;2k
be the symplectic Lanczos factorization after k steps, with B 2k;2k unreduced. Then
Jordan matrix of order 2k.
Proof. If d
XJ be the Jordan canonical form of B 2k;2k
and put
X. Then
Suppose now that
it follows that
Hence by (3.6) dim K(M;
unreduced, dim K(M; k. Hence dim K(M;
and therefore, d
A similar result may be formulated in terms of Schur vectors or symplectic Schur
vectors (see, e.g., [28, 29] for the real symplectic Schur decomposition of a symplectic
matrix). These theorems provide the motivation for the implicit restart developed
in the next section. Theorem 3.3 suggests that one might find an invariant subspace
by iteratively replacing the starting vector with a linear combination of approximate
eigenvectors corresponding to eigenvalues of interest. Such approximations are readily
available through the Lanczos factorization.
3.3. Stopping Criteria. Now assume that we have performed k steps of the
symplectic Lanczos method and thus obtained the identity (after permuting back)
MS
An implicitly restarted symplectic Lanczos method 13
If the norm of the residual vector is small, the 2k eigenvalues of B 2k;2k are approximations
to the eigenvalues of M . Numerical experiments indicate that the norm of
the residual rarely becomes small by itself. Nevertheless, some eigenvalues of B 2k;2k
may be good approximations to eigenvalues of M . Let - be an eigenvalue of B 2k;2k
with the corresponding eigenvector y. Then the vector
The vector x is referred to as Ritz vector and - as Ritz value of M . If the last
component of the eigenvector y is sufficiently small, the right-hand side of (3.11) is
small and the pair f-; xg is a good approximation to an eigenvalue-eigenvector pair
of M . Note that by Lemma 2.2 je T
is unreduced. The pair (-; x) is
exact for the nearby problem
A small jjEjj is not sufficient for the Ritz pair f-; xg being a good approximation
to an eigenvalue-eigenvector pair of M . The advantage of using the Ritz estimate
jd
w k+1 jj is to avoid the explicit formation of the residual
deciding about the numerical accuracy of an
approximate eigenpair.
It is well-known that for non-normal matrices the norm of the residual of an
approximate eigenvector is not by itself sufficient information to bound the error in
the approximate eigenvalue. It is sufficient however to give a bound on the distance
to the nearest matrix to which the given approximation is exact. In the following, we
will give a computable expression for the error. Assume that B 2k;2k is diagonalizable
Since MS
2k , it follows that
MS
or
2k Y: Thus
k. The last equation can be re-written as
Using Theorem 2' of [19] we obtain that (- is an eigen-triplet of
where
jjJx k+i jj g:
14 Benner and Fa-bender
Furthermore, when jjEjj is small enough, then
is an eigenvalue of M and
Consequently, the symplectic Lanczos algorithm should be continued until both jjEjj
is small and cond(- j )jjEjj is below a given threshold for accuracy.
4. An Implicitly Restarted Symplectic Lanczos Method. In the previous
sections we have briefly mentioned two algorithms for computing approximations
to the eigenvalues of a symplectic matrix M . The symplectic Lanczos algorithm is
appropriate when the matrix M is large and sparse. If only a small subset of the
eigenvalues is desired, the length k symplectic Lanczos factorization may suffice. The
analysis in the last chapter suggests that a strategy for finding 2k eigenvalues in a
length k factorization is to find an appropriate starting vector that forces the residual
r k+1 to vanish. The SR algorithm, on the other hand, computes approximations to all
eigenvalues and eigenvectors of M . From Theorem 4.1 in [4] (an implicit Q-theorem
for the SR case) we know that in exact arithmetic, when using the same starting
vector, the SR algorithm and the length n Lanczos factorization generate the same
symplectic butterfly matrices (up to multiplication by a trivial matrix). Forcing the
residual for the symplectic Lanczos algorithm to zero has the effect of deflating a subdiagonal
element during the SR algorithm: by Remark 3.1 r
from the symplectic Lanczos process we have d Hence a zero residual
implies a zero d k+1 such that deflation occurs for the corresponding butterfly matrix.
Our goal in this section will be to construct a starting vector that is a member of
the invariant subspace of interest. Our approach is to implicitly restart the symplectic
Lanczos factorization. This was first introduced by Sorensen [35] in the context of
unsymmetric matrices and the Arnoldi process. The scheme is called implicit because
the updating of the starting vector is accomplished with an implicit shifted SR mechanism
on This allows to update the starting vector by working with a
symplectic matrix in IR 2j \Theta2j rather than in IR 2n\Theta2n which is significantly cheaper.
The iteration starts by extending a length k symplectic Lanczos factorization by
steps. Next, 2p shifts are applied to B 2(k+p);2(k+p) using double or quadruple SR
steps. The last 2p columns of the factorization are discarded resulting in a length k
factorization. The iteration is defined by repeating this process until convergence.
For simplicity let us first assume that and that a 2n \Theta 2(k
P is known such that
as in (3.5). Let - be a real shift and
Then (using
will be a permuted butterfly matrix and SP is
an upper triangular matrix with two additional subdiagonals.
With this we can re-express (4.1) as
MP (S 2n;2k+2
An implicitly restarted symplectic Lanczos method 15
P SP this yields
The above equation fails to be a symplectic Lanczos factorization since the columns
of the matrix d k+2 (b k+2 v k+2
2k+2 SP are nonzero. Let
ij be the (i; j)th entry of SP . The residual term in (4.2) is
Rewriting (4.2) as
where Z is blocked as6 6 6 6 4
dk+1e T
dk+1e T
dk+2 bk+2s 2k+2;2k e T
dk+2ak+2s 2k+2;2k e T
we obtain as a new Lanczos identity
r
where
d
a
Here, - a k+1 , - b k+1 , -
d k+1 denote parameters of -
are
parameters of B 2k+2;2k+2
P . In addition, -
w k+1 are the last two column vectors
from -
are the two last column vectors of S 2n;2k+2
As the space spanned by the columns of S
orthogonal, and SP is a permuted symplectic matrix, the space spanned by the
columns of -
is J-orthogonal. Thus (4.3) is a valid symplectic
Lanczos factorization. The new starting vector is -
ae 2 IR. This can be seen as follows: first note that for unreduced butterfly matrices
B 2k+2;2k+2 we have q 2 (B 2k+2;2k+2
Hence, from q 2 (B 2k+2;2k+2
we obtain q 2 (B 2k+2;2k+2
is an upper triangular
matrix. As q 2 (B 2k+2;2k+2
Using (4.3) it follows that
ae S 2n;2k+2
ae S 2n;2k+2
=ae (MP S 2n;2k+2
=ae (MP S 2n;2k+2
Benner and Fa-bender
as r k+2 e T
using again (4.3) we get
\Gammaae
as e T
Note that in the symplectic Lanczos process the vectors v j of S 2n;2k
P satisfy the
condition and the parameters b j are chosen to be one. This is no longer
true for the odd numbered column vectors of SP generated by the SR decomposition
and the parameters - b j from -
P and thus for the new Lanczos factorization (4.3).
Both properties could be forced using trivial factors. Numerical tests indicate that
there is no obvious advantage in doing so.
Using standard polynomials as shift polynomials instead of Laurent polynomials
as above results in the following situation: In p 2 (B 2k+2;2k+2
is an upper triangular matrix with four (!) additional
subdiagonals. Therefore, the residual term in (4.2) has five nonzero entries.
Hence not the last two, but the last four columns of (4.2) have to be discarded in
order to obtain a new valid Lanczos factorization. That is, we would have to discard
wanted information which is avoided by using Laurent polynomials.
This technique can be extended to the quadruple shift case using Laurent polynomials
as the shift polynomials as discussed in Section 2. The implicit restart can
be summarized as given in Table 4.1. In the course of the iteration we have to choose
shifts in order to apply 2p shifts: choosing a real shift - k implies
that - \Gamma1
k is also a shift due to the symplectic structure of the problem. Hence, - \Gamma1
k is
not added to \Delta as the use of the Laurent polynomial q 2 guarantees that - \Gamma1
k is used
as a shift once - k 2 \Delta. In case of a complex shift - k , implies that - k is
also a shift not added to \Delta. For complex shifts - k ,
- k in \Delta.
Numerous choices are possible for the selection of the p shifts. One possibility is
the case of choosing p "exact" shifts with respect to B 2(k+p);2(k+p)
. That is, first the
eigenvalues of B 2(k+p);2(k+p)
are computed (by the SR algorithm), then p unwanted
eigenvalues are selected. One choice for this selection might be: sort the eigenvalues
by decreasing magnitude. There will be k eigenvalues with modulus greater than
or equal to 1
Select the 2p eigenvalues with modulus closest to 1 as shifts. If - k+1 is complex with
then we either have to choose 2p shifts or just 2p
as - k+1 belongs to a quadruple pair of eigenvalues of B 2(k+p);2(k+p)
P and in order to
preserve the symplectic structure either - k and - k+1 have to be chosen or none.
An implicitly restarted symplectic Lanczos method 17
restarted symplectic Lanczos method
perform k steps of the symplectic Lanczos algorithm to compute S 2n;2k
obtain the residual vector r k+1
while jjr k+1 jj ? tol
perform p additional steps of the symplectic Lanczos method
to compute S 2n;2(k+p)
select p shifts - i
compute -
via implicitly shifted SR steps
set S 2n;2k
obtain the new residual vector r k+1
while

Table
k-step restarted symplectic Lanczos method
A different possibility of choosing the shifts is to keep those eigenvalues that are
good approximations to eigenvalues of M . That is, eigenvalues for which (3.11) is
small. Again we have to make sure that our set of shifts is complete in the sense
described above.
Choosing eigenvalues of B 2(k+p);2(k+p)
P as shifts has an important consequence
for the next iterate. Assume for simplicity that B 2(k+p);2(k+p)
P is diagonalizable. Let
be a disjoint partition of the spectrum
of B 2(k+p);2(k+p)
. Selecting the exact shifts - in the implicit restart,
following the rules mentioned above yields a matrix
g. This follows from (2.6).
Moreover, the new starting vector has been implicitly replaced by the sum of 2k
approximate eigenvectors:
ae
ae
properly chosen. The last
equation follows since q(B 2(k+p);2(k+p)
)e 1 has no component along an eigenvector of
associated with Hence
It should be mentioned that the k-step restarted symplectic Lanczos method as
in

Table

4.1 with exact shifts builds a J-orthogonal basis for a number of generalized
Krylov subspaces simultaneously. The subspace of length 2(k +p) generated during a
restart using exact shifts contains all the Krylov subspaces of dimension 2k generated
from each of the desired Ritz vectors, for a detailed discussion see [13]. A similar
Benner and Fa-bender
observation for Sorensen's restarted Arnoldi method with exact shifts was made by
Morgan in [30]. For a discussion of this observation see [30] or [23]. Morgan infers
'the method works on approximations to all of the desired eigenpairs at the same time,
without favoring one over the other' [30, p. 1220,l. 7-8 from the bottom]. This remark
can also be applied to the method presented here.
In the above discussion we have assumed that the permuted SR decomposition
exists. Unfortunately, this is not always true. During the
bulge-chase in the implicit SR step, it may happen that a diagonal element a j of B 1
(2.2) is zero (or almost zero). In that case no reduction to symplectic butterfly form
with the corresponding first column -
does exist. In the next section we will prove
that a serious breakdown in the symplectic Lanczos algorithm is equivalent to such
a breakdown of the SR decomposition. Moreover, it may happen that a subdiagonal
element d j of the (2; 2)-block of B
2 (2.3) is zero (or almost zero) such that
The matrix -
P is split, an invariant subspace of dimension j is found. If
shifts have been applied, then the iteration is halted. Otherwise we
continue similar to the procedure described by Sorensen in [35, Remark 3].
As the iteration progresses, some of the Ritz values may converge to eigenvalues of
long before the entire set of wanted eigenvalues have. These converged Ritz values
may be part of the wanted or unwanted portion of the spectrum. In either case it
is desirable to deflate the converged Ritz values and corresponding Ritz vectors from
the unconverged portion of the factorization. If the converged Ritz value is wanted
then it is necessary to keep it in the subsequent factorizations; if it is unwanted then
it must be removed from the current and the subsequent factorizations. Lehoucq and
Sorensen develop in [23, 36] locking and purging techniques to accomplish this in the
context of unsymmetric matrices and the restarted Arnoldi method. These ideas can
be carried over to the situation here.
5. Numerical Properties of the Implicitly Restarted Symplectic Lanczos
Method.
5.1. Stability Issues. It is well known that for general Lanczos-like methods
the stability of the overall process is improved when the norm of the Lanczos vectors is
chosen to be equal to 1 [32, 37]. Thus, Banse proposes in [2] to modify the prerequisite
our symplectic Lanczos method to
\Gammaoe
and
An implicitly restarted symplectic Lanczos method 19
For the resulting algorithm and a discussion of it we refer to [2]. It is easy to see that
BP SP is no longer a permuted symplectic matrix, but it still has the desired form
of a butterfly matrix. Unfortunately, an SR step does not preserve the structure of
and thus, this modified version of the symplectic Lanczos method can not
be used in connection with our restart approaches.
some form of reorthogonalization any Lanczos algorithm is numerically
unstable. Hence we re-J P -orthogonalize each Lanczos vector as soon as it is computed
against the previous ones via
wm
where for
defines the indefinite inner product implied
by J n
This re-J P -orthogonalization is costly, it requires 16n(m \Gamma 1) flops for the vector
wm and 16nm flops for v m+1 . Thus, if 2k Lanczos vectors are
computed, the re-J P -orthogonalization adds a computational cost of the order of
flops to the overall cost of the symplectic Lanczos method.
For standard Lanczos algorithms, different reorthogonalization techniques have
been studied (for references see, e.g., [17]). Those ideas can be used to design analogous
re-J P -orthogonalizations for the symplectic Lanczos method. It should be noted
that if k is small, the cost for re-J P -orthogonalization is not too expensive.
Another important issue is the numerical stability of the SR step employed in
the restart. During the SR step on the 2k \Theta 2k symplectic butterfly matrix, all but
are orthogonal. These are known to be numerically stable. For
the nonorthogonal symplectic transformations that have to be used, we choose
among all possible transformations the ones with optimal (smallest possible) condition
number (see [8]).
5.2. Breakdowns in the SR Factorization. If there is a starting vector -
aeq(M)v 1 for which the explicitly restarted symplectic Lanczos method breaks down,
then it is impossible to reduce the symplectic matrix M to symplectic butterfly form
with a transformation matrix whose first column is - v 1 . Thus, in this situation the SR
decomposition of q(B) can not exist.
As will be shown in this section, this is the only way that breakdowns in the
SR decomposition can occur. In the SR step, most of the transformations used are
orthogonal symplectic transformations; their computation can not break down. The
only source of breakdown can be one of the symplectic Gaussian eliminations L j .
For simplicity, we will discuss the double shift case. Only the following elementary
elimination matrices are used in the implicit SR step: elementary symplectic Givens
matrices [31]
where
20 Benner and Fa-bender
elementary symplectic Householder transformation
and elementary symplectic Gaussian elimination matrices [8]
where
Assume that k steps of the symplectic Lanczos algorithm are performed, then
from (3.5)
Now an implicit restart is to be performed using an implicit double shift SR step. In
the first step of the implicit SR step, a symplectic Householder matrix H 1 is computed
such that
H 1 is applied to B 2k;2k
introducing a small bulge in the butterfly form: additional elements are found in the
positions (2; 1), (1; 2), (n
1). The remaining implicit transformations perform a bulge-chasing
sweep down the subdiagonal to restore the butterfly form. An algorithm for this is
given in [2] or [4]; it can be summarized for the situation here as in Table 5.1, where
~
G j and G j both denote symplectic Givens transformation matrices acting in the same
planes but with different rotation angles.
compute G '+1 such that (G '+1 B 2k;2k )
compute L '+1 such that (L '+1 B 2k;2k )
compute ~
G '+1 such that (B 2k;2k ~
G '+1
compute H '+1 such that (B 2k;2k H '+1

Table
Reduction to butterfly form - double shift case.
An implicitly restarted symplectic Lanczos method 21
Suppose that the first exist and that we
have computed
~
In order to simplify the notation, we switch to the permuted version and rewrite the
permuted symplectic matrix b
SP as
I 2n\Gamma2j \Gamma2
making use of the fact that the accumulated transformations
affect only the rows 1 to j and j. The leading (2j
principal submatrix of
is given by
e
x
x
x
x
x
where the hatted quantities denote unspecified entries that would change if the SR
update could be continued. Next, the (2j should be annihilated
by a permuted symplectic Gaussian elimination. This elimination will fail to exist if
the SR decomposition of q(B 2k;2k ) does not exist.
As will be needed later, - a implies that - y This follows as e
P is
From
e
we obtain
- a j
x
x
If - a
(otherwise the last Gaussian transformation
did not exist).
Next we show that this breakdown in the SR decomposition implies a breakdown
in the Lanczos process started with the starting vector -
22 Benner and Fa-bender
For this we have to consider (5.1) multiplied from the right by b
SP . From the
derivations in the last section we know that the starting vector of that recursion is
given by -
As the trailing (2n
submatrix of b
SP is the identity, we can just as well consider
multiplied from the right by SP
P SP corresponds to the matrix in (5.2) (no butterfly
w j+1 ]. The
columns of -
are JP -orthogonal
The starting vector of the recursion (5.3) is given by - Deleting the
last four columns of -
P in the same way as in the implicit restart we obtain a
valid symplectic Lanczos factorization of length 2.
In order to show that a breakdown in the SR decomposition of q(B) implies a
breakdown in the above symplectic Lanczos recursion, we need to show
From (5.2) and (5.3) we obtain
and
Further we do know from the symplectic Lanczos algorithm
all of these quantities are already known. Now consider
x3
Obviously, x Using (5.6) we obtain
2. Hence x Using (5.5) end (5.4) will see that x
z3
An implicitly restarted symplectic Lanczos method 23
As - a
From (5.3) we obtain
Hence using (5.4) yields
Similar, it follows that z
This argumentation has shown that an SR breakdown implies a serious Lanczos
breakdown. The opposite implication follows from the uniqueness of the Lanczos
factorization. The result is summarized in the following theorem.
Theorem 5.1. Suppose the symplectic butterfly matrix B 2k;2k corresponding to
(3.5) is unreduced and let - 2 IR. Let L j be the jth symplectic Gauss transformation
required in the SR step on (B If the first
symplectic Gauss transformations of this SR step exist, then L j fails to exist if and
only if - v T
j as in (4.3).
6. Numerical Experiments. Some examples to demonstrate the properties of
the (implicitly restarted) symplectic Lanczos method are presented. The computational
results are quite promising but certainly preliminary. All computations were
done using Matlab Version 5.1 on a Sun Ultra 1 with IEEE double-precision arithmetic
and machine precision
Our code implements exactly the algorithm as given in Table 4.1. In order to
detect convergence in the restart process, the rather crude criterion
was used. This ad hoc stopping rule allowed the iteration to halt quite early. Usually,
the eigenvalues largest in modulus (and their reciprocals) of the wanted part of the
spectrum are much better approximated than the ones of smaller modulus. In a black-box
implementation of the algorithm this stopping criterion has to be replaced with
a more rigorous one to ensure that all eigenvalues are approximated to the desired
accuracy (see the discussion in Section 3.3). Benign breakdown in the symplectic
Lanczos process was detected by the criterion
while a serious breakdown was detected by
Our implementation intends to compute the k eigenvalues of M largest in modulus
and their reciprocals. In the implicit restart, we used exact shifts where we chose the
shifts to be the 2p eigenvalues of B 2k+p;2k+p closest to the unit circle.
Our observations have been the following.
Benner and Fa-bender
ffl Re-J-orthogonalization is necessary; otherwise J-orthogonality of the computed
Lanczos vectors is lost after a few steps, and ghost eigenvalues (see,
e.g., [17]) appear. That is, multiple eigenvalues of B 2k;2k correspond to simple
eigenvalues of M .
ffl The implicit restart is more accurate than the explicit one.
ffl The leading end of the 'wanted' Ritz values (that is, the eigenvalues largest
in modulus and their reciprocals) converge faster than the tail end (closest to
cut off of the sort). The same behavior was observed in [35] for the implicitly
restarted Arnoldi method. In order to obtain faster convergence, it seems
advisable (similar to the implementation of Sorensen's implicitly restarted
Arnoldi method in Matlab's eigs) to increase the dimension of the computed
Lanczos factorization. That is, instead of computing S 2n;2k
as a basis for the restart, one should compute a slightly larger factorization,
e.g. dimension 2(k instead of dimension 2k. When 2' eigenvalues have
converged, a subspace of dimension 2(k computed as a basis for
the restart, followed by p additional Lanczos steps to obtain a factorization
of length k Using implicit SR steps this factorization is reduced
to one of length k If the symplectic Lanczos method would be implemented
following this approach, the convergence check could be done using
only the k Ritz values of largest modulus (and their reciprocals) or those that
yield the smallest Ritz residual
jd
where the y j are the eigenvectors of B 2k;2k .
ffl It is fairly difficult to find a good choice for k and p. Not for every possible
choice of k, there exists an invariant subspace of dimension 2k associated to
the k eigenvalues - i largest in modulus and their reciprocals. If - k is complex
and - then we can not choose the 2p eigenvalues with modulus
closest to the unit circle as shifts as this would tear a quadruple of eigenvalues
apart resulting in a shift polynomial q such that q(B 2(k+p);2(k+p)
we can do is to choose the 2p \Gamma 2 eigenvalues with modulus closest to 1
as shifts. In order to get a full set of 2p shifts we add as the last shift
the real eigenvalue pair with largest Ritz residual. Depending on how good
that real eigenvalue approximates an eigenvalue of M , this strategy worked,
but the resulting subspace is no longer the subspace corresponding to the k
eigenvalues largest in modulus and their reciprocals. If the real eigenvalue
has converged to an eigenvalue of M , it is unlikely to remove that eigenvalue
just by restarting, it will keep coming back. Only a purging technique like the
one discussed by Lehoucq and Sorensen [23, 36] will be able to remove this
eigenvalue. Moreover, there is no guarantee that there is a real eigenvalue of
P that can be used here. Hence, in a black-box implementation
one should either try to compute an invariant subspace of dimension
or of dimension 2(k 1). As this is not known a priori, the algorithm should
adapt k during the iteration process appropriately. This is no problem, if as
suggested above, one always computes a slightly larger Lanczos factorization
than requested.
Example 6.1. The first test performed concerned the loss of J-orthogonality of
the computed Lanczos vectors during the symplectic Lanczos method and the ghost
An implicitly restarted symplectic Lanczos method 25
eigenvalue problem (see, e.g. [17]). To demonstrate the effects of re-J-orthogonali-
zation, a 100 \Theta 100 symplectic matrix with eigenvalues
200; 100; 50;
was used. A symplectic block-diagonal matrix with these eigenvalues on the block-diagonal
was constructed and a similarity transformation with a randomly generated
orthogonal symplectic matrix was performed to obtain a symplectic matrix M .
As expected, when using a random starting vector M 's eigenvalues largest in
modulus (and the corresponding reciprocals) tend to emerge right from the start,
e.g., the eigenvalues of B 10;10 are
199:99997; 100:06771; 48:71752; 26:85083; 8:32399
and their reciprocals. Without any form of re-J-orthogonalization, the J-orthogo-
nality of the Lanczos vectors is lost after a few iterations as indicated in Figure 6.1.
number of Lanczos steps
100,2k
JS
100,2k
||Fig. 6.1. loss of J-orthogonality after k symplectic Lanczos steps
The loss of J-orthogonality in the Lanczos vectors results, as in the standard
Lanczos algorithm, in ghost eigenvalues. That is, multiple eigenvalues of B 2k;2k correspond
to simple eigenvalues of M . For example, using no re-J-orthogonalization,
after 17 iterations the 6 eigenvalues largest in modulus of B 34;34 are
Using complete re-J-orthogonalization, this effect is avoided:
200; 100; 49:99992; 47:02461; 45:93018; 42:31199:
The second test performed concerned the question whether an implicit restart
is more accurate than an explicit one. After nine steps of the symplectic Lanczos
method (with a random starting vector) the resulting butterfly
had the eigenvalues (using the Matlab function eig)
200:000000000000 99:999999841718
13:344815062428 3:679215125563 \Sigma5:750883779240i
26 Benner and Fa-bender
and their reciprocals. Removing the 4 complex eigenvalues from B 18;18 using an
implicit restart as described in Section 4, we obtain a symplectic butterfly matrix
impl whose eigenvalues are
200:000000000000 99:999999841719
13:344815062428
and their reciprocals. From (2.6) it follows that these have to be the 14 real eigenvalues
of B 18;18 which have not been removed. As can be seen, we lost one digit during
the implicit restart (indicated by the 'underbar' under the 'lost' digits in the above
table). Performing an explicit restart with the explicitly computed new starting vector
butterfly
expl whose eigenvalues are
200:000000000000 99:999999841793
and their reciprocals. This time we lost up to nine digits.
The last set of tests performed on this matrix concerned the k-step restarted
symplectic Lanczos method as given in Table 4.1. As M has only one quadruple
of complex eigenvalues, and these eigenvalues are smallest in magnitude there is no
problem in choosing k - n. For every such choice there exists an invariant symplectic
subspace corresponding to the k eigenvalues largest in magnitude and their reciprocals.
In the tests reported here, a random starting vector was used. Figure 6.2 shows a plot
of jjr k+1 jj versus the number of iterations performed. Iteration Step 1 refers to the
norm of the residual after the first k Lanczos steps, no restart is performed. The three
lines in Figure 6.2 present three different choice for k and p:
Convergence was achieved for all three examples (and many more,
not shown here). Obviously, the choice results in faster convergence
than the choice 8. Convergence is by no means monotonic, during the major
part of the iteration the norm of the residual is changing quite dramatically. But once
a certain stage is achieved, the norm of the residual converges. Although convergence
quite fast, this does not imply that convergence is
as fast for other choices of k and p. The third line in Figure 6.2 demonstrates that
the convergence for does need twice as many iteration steps as for
Example 6.2. Symplectic matrix pencils that appear in discrete-time linear-quadratic
optimal control problems are typically of the form
- I \GammaBB T
(Note: For F 6= I , L and N are not symplectic, but L \Gamma -N is a symplectic matrix
pencil.) Assuming that L and N are nonsingular (that is, F is nonsingular), solving
this generalized eigenproblem is equivalent to solving the eigenproblem for the
symplectic matrix
- I \GammaBB T
If one is interested in computing a few of the eigenvalues of L \Gamma -N , one can use the
An implicitly restarted symplectic Lanczos method 27
number of iterations
norm(r
Fig. 6.2. k-step restarted symplectic Lanczos method, different choices of k and p
restarted symplectic Lanczos algorithm on In each step of the symplectic
Lanczos algorithm, one has to compute matrix-vector products of the form Mx and
Making use of the special form of L and N this can be done without explicitly
inverting us consider the computation of First compute
Next one has to solve the
linear system analogous to x and z, then from Ny = z
we obtain
In order to solve y we compute the LU decomposition of F and solve the
linear system F T y using backward and forward substitution. Hence, the explicit
inversion of N or F is avoided. In case F is a sparse matrix, sparse solvers can be
employed. In particular, if the control system comes from some sort of discretization
scheme, F is often banded which can be used here by computing an initial band LU
factorization of F in order to minimize the cost for the computation of y 2 . Note that
in most applications, such that the computational cost for C T Cx 1 and
significantly cheaper than a matrix-vector product with an n \Theta n matrix.
In case of single-input the corresponding operations
come down to two dot products of length n each.
Using Matlab's sparse matrix routine sprandn sparse normally distributed random
matrices F; B; C (here, n) of different dimensions and with different
densities of the nonzero entries were generated. Here an example of dimension
presented, where the density of the different matrices was chosen to be
matrix - nonzero entries
28 Benner and Fa-bender
Matlab computed the norm of the corresponding matrix to be - 5:3 \Theta
In the first set of tests k was chosen to be 5, and we tested
As can be seen in Figure 6.3, for the first 3 iterations, the norm of the residual
decreases for both choice of p, but then increases quite a bit. During the first step,
the eigenvalues of B 10;10 are approximating the 5 eigenvalues of L \Gamma -N largest in
modulus and their reciprocals. In step 4, a 'wrong' choice of the shifts is done in
both cases. The extended matrices B 20;20 and B 30;30 both still approximate the 5
eigenvalues of L \Gamma -N largest in modulus, but there is a new real eigenvalue coming
in, which is not a good approximation to an eigenvalue of L \Gamma -N . But, due to the
way the shifts are chosen here, this new eigenvalue is kept, while an already good
approximated eigenvalue - a little smaller in magnitude - is shifted away, resulting
in a dramatic increase of jjr k+1 jj. Modifying the choice of the shifts such that the
good approximation is kept, while the new real eigenvalue is shifted away, the problem
is resolved, the 'good' eigenvalues are kept and convergence occurs in a few steps (the
'o'-line in Figure 6.3).
Using a slightly larger Lanczos factorization as a basis for the restart, e.g., a
factorization of length k + 3 instead of length k and using a locking technique to
decouple converged approximate eigenvalues and associated invariant subspaces from
the active part of the iteration, this problem is avoided.
number of iterations
norm(r
modified
Fig. 6.3. k-step restarted symplectic Lanczos method, different choices of the shifts

Figure

6.4 displays the behavior of the k-step restarted symplectic Lanczos method
for different choices of k and p, where k is quite small. Convergence is achieved in
any case.
So far, in the tests presented, k was always chosen such that there exists a deflating
subspace of L \Gamma -N corresponding to the k eigenvalues largest in modulus and their
reciprocals. For there is no such deflating subspace (there is one for
and one for Figure 6.5 for a convergence plot. The eigenvalues of
B 2(k+p);2(k+p) in the first iteration steps approximate the k eigenvalues of largest
modulus and their reciprocals (where 5 - j - p) quite well. Our choice of shifts is
to select the 2p eigenvalues with modulus closest to 1, but as - k+1 is complex with
1, we can only choose shifts that way. The last shift is chosen
An implicitly restarted symplectic Lanczos method 29
number of iterations
norm(r
Fig. 6.4. k-step restarted symplectic Lanczos method, different choices of k and p
according to the strategy explained above. This eigenvalue keeps coming back before
it is annihilated. A better idea to resolve the problem is to adapt k appropriately.
number of iterations
norm(r
Fig. 6.5. k-step restarted symplectic Lanczos method, different choices of k and p
7. Concluding Remarks. We have investigated a symplectic Lanczos method
for symplectic matrices. Employing the technique of implicitly restarting the method
using double or quadruple shifts as zeros of the driving Laurent polynomials, this
results in an efficient method to compute a few extremal eigenvalues of symplectic
matrices and the associated eigenvectors or invariant subspaces. The residual of the
Lanczos recursion can be made to zero by choosing proper shifts. It is an open problem
how these shifts should be chosen in an optimal way. The preliminary numerical tests
reported here show that for exact shifts, good performance is already achieved.
Before implementing the symplectic Lanczos process in a black-box algorithm,
some more details need consideration: in particular, techniques for locking of con-
Benner and Fa-bender
verged Ritz values as well as purging of converged, but unwanted Ritz values, needs
to be derived in a similar way as it has been done for the implicitly restarted Arnoldi
method.



--R

analysis of the Lanczos algorithm for the nonsymmetric eigenvalue problem
Symplektische Eigenwertverfahren zur L-osung zeitdiskreter optimaler Steuerungs- probleme
A condensed form for the solution of the symplectic eigenvalue problem
The symplectic eigenvalue problem
SR and SZ algorithms for the symplectic (butterfly) eigenproblem
Linear Hamiltonian difference systems: Disconjugacy and Jacobi-type conditions
Matrix factorization for symplectic QR-like methods
A symplectic QR-like algorithm for the solution of the real algebraic Riccati equation
An implicitly restarted Lanczos method for large symmetric eigenvalue problems
Sur quelques Algorithmes de recherche de valeurs propres
Numerical linear algorithms and group theory
On some algebraic problems in connection with general eigenvalue algorithms
Symplectic Methods for Symplectic Eigenproblems
An analysis of structure preserving methods for symplectic eigenvalue problems
The QR transformation Part I and Part II

Matrix Computations
Model reduction of state space systems via an implicitly restarted Lanczos method
Residual bounds on approximate eigensystems of nonnormal matrices
On some algorithms for the solution of the complete eigenvalue problem
The Algebraic Riccati Equation
Invariant subspace methods for the numerical solution of Riccati equations
Deflation techniques for an implicitly restarted Arnoldi itera- tion
Solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods
A new method for computing the closed loop eigenvalues of a discrete-time algebraic Riccati equation



Canonical forms for Hamiltonian and symplectic matrices and pencils
On restarting the Arnoldi method for large nonsymmetric eigenvalue problems
A Schur decomposition for Hamiltonian matrices

Computation of the stable deflating subspace of a symplectic pencil using structure preserving orthogonal transformations

Implicit application of polynomial filters in a k-step Arnoldi method

Analysis of the look ahead Lanczos algorithm
A symplectic method for approximating all the eigenvalues of a Hamiltonian matrix
A convergence analysis for nonsymmetric Lanczos algorithms
--TR
implicit restarting;symplectic Lanczos method;symplectic matrix;eigenvalues