Numerical Algorithms 18 (1998) 31–50

31

Nested Lanczos: implicitly restarting an unsymmetric Lanczos algorithm Gorik De Samblanx and Adhemar Bultheel K.U. Leuven, Department of Computer Science, Celestijnenlaan 200A, B-3001 Heverlee, Belgium

Received 12 September 1997; revised 24 June 1998 Communicated by C. Brezinski

In this text, we present a generalization of the idea of the Implicitly Restarted Arnoldi method to the unsymmetric Lanczos algorithm, using the two-sided Gram–Schmidt process or using a full Lanczos tridiagonalization. The resulting implicitly restarted Lanczos method is called Nested Lanczos. Nested Lanczos can be combined with an implicit filter. It can also be used in case of breakdown and offers an alternative for look-ahead. Keywords: unsymmetric Lanczos, Implicitly Restarted Arnoldi, two-sided Gram–Schmidt AMS subject classification: 65F15

1.

Introduction

In the last few years, the concept of restarting an iterative eigenvalue solver has been widely accepted. Restarting an algorithm – implicitly or explicitly – provides us with a solution for the specific problems that emerge with the use of iterative solvers. These problems can be: slow convergence along with growing computational work per iteration step, the occurrence of spurious solutions or an uncertain robustness. For Arnoldi’s method, the idea of explicitly restarting the factorization was introduced by Saad [15] and different strategies were formulated on how the restart must be applied in practice. The main problem was to find a set of new starting conditions that contains as much “useful” information as possible. A breakthrough was accomplished with the Implicitly Restarted Arnoldi method (IRA) of Sorensen [16], which led to the popular ARPACK eigenvalue package [10]. One can say that the advantage of implicitly restarting lies in the fact that it removes uninteresting information instead of trying to retain the interesting part. This approach appeared not only to be cheaper, but it gave a new insight into iterative methods. The implicitly restarting Arnoldi procedure is normally cheaper than an explicit one, since it does not require to redo the matrix–vector products and the orthogonalizations of the original method. Instead, it performs a QR factorization of the Hessenberg matrix, i.e., a special Arnoldi run. It can often be combined with an implicit filtering step. Restarting the Arnoldi relation can reduce the amount of required storage space  J.C. Baltzer AG, Science Publishers

32

G. De Samblanx, A. Bultheel / Nested Lanczos

for the algorithm, and it can remove lumber information when the method converges too slowly. The filtering can be used in order to remove spurious – but dominant – information that may mislead the algorithm. It also increases the robustness of the algorithm by implementing a validation strategy, or it can be used for acceleration. Many people studied the practical aspects of restarting [11,17]: how many vectors should be removed, which vectors contain the most information, etc. We show in this text how the idea of IRA can be generalized to the unsymmetric Lanczos algorithm. For the symmetric case, Arnoldi corresponds to Lanczos and the theory is much alike [5]. For the unsymmetric case, Grimme et al. [7] considered restarting a sign symmetric implementation of the Lanczos algorithm in order to compute a stable reduced order model of a large state space control system. Our approach is more general in that it can be applied to any Lanczos implementation. We view it here only in the context of finding solutions to eigenvalue problems. Restarting the Lanczos algorithm has an additional application that might be very important. If there is a breakdown in the Lanczos algorithm, then the restart can be used in order to remove this singularity. The restart algorithm provides us with an alternative for a look-ahead algorithm that does not introduce a block structure in the Lanczos matrix. Instead, it looks back on previous computations by applying a transformation on the basis vectors. The aim of this paper is to show how the unsymmetric Lanczos algorithm can be restarted implicitly in a most general way. However, it can be seen as a generalization of IRA, it contains new and unexpected aspects, e.g., the application in case of breakdown. In section 2, we briefly recall the unsymmetric Lanczos algorithm and the different types of breakdown. In section 3, we restart the Lanczos algorithm by using a biorthogonal factorization of the small, projected problem. This factorization can be computed by a Lanczos algorithm with special starting vectors or by the two-sided Gram–Schmidt algorithm (BioGS). Both approaches are mathematically equivalent. Section 4 focuses on the “classical” properties of an implicitly restarted method: the filtering property and the use of “exact shifts” in order to remove eigenvalues from the approximation. These properties are used to show that the computation of the implicit restart does not break down itself, unless the explicit restart with corresponding starting vectors would break down. In section 5, we show that implicitly restarted Lanczos may be used as an alternative for look-ahead, in case of breakdown. Look-ahead extends the approximating subspaces until it finds a set that is no longer singular. Implicitly restarted Lanczos chooses a number of vectors from this large set, which form a nonsingular pair of subspace bases. In section 6, we close the text with some conclusions. Notation. Roman characters denote matrices and vectors; the (i, j)th entry of a matrix Q is denoted by qi,j . The hermitian transpose is given by Q∗ . Greek characters denote scalars and α ¯ is the complex conjugate of α. The kth unit vector is ek , the identity matrix is Ik,l ∈ Ck×l and the column range of a matrix is denoted by R(Q). Two matrices Q and Z are called biorthogonal if Q∗ Z = I. If they are square, then Q∗ Z = I = QZ ∗ .

G. De Samblanx, A. Bultheel / Nested Lanczos

2.

33

The unsymmetric Lanczos algorithm

Before we show how to restart the Lanczos algorithm, we describe the algorithm itself [9]. The unsymmetric Lanczos algorithm reduces a matrix to tridiagonal form; however, the result cannot be guaranteed for an arbitrary matrix. There are many different implementations possible, but all are equivalent up to some scaling of the Lanczos vectors [13, theorem 2.2]. The restarting procedure that we present does not depend on the exact implementation of the algorithm, so we are free to choose one. It is well known that if the Lanczos algorithm converges, then the biorthogonalization can become erroneous [1]. We assume in the text that this problem is settled, e.g., by using some reorthogonalization steps. The Lanczos algorithm projects a complex, unsymmetric matrix A with an oblique projector. The projector is defined by a biorthogonal pair of matrices (Vk , Wk ), which represent the basis of the subspaces in which the left and right eigenvectors are approximated. The columns of the basis matrices are called Lanczos vectors. The range of the Lanczos vectors depends on the first columns of Vk = [v1 , . . . , vk ] and Wk = [w1 , . . . , wk ]. They correspond each to a Krylov subspace  R(Vk ) = Kk (A, v1 ) = v1 , Av1 , . . . , Ak−1 v1 and (1)

k−1 R(Wk ) = Kk A∗ , w1 = w1 , A∗ w1 , . . . , A∗ w1 . (2) The unsymmetric Lanczos algorithm that we will use is given by algorithm 1. Unless it breaks down, the algorithm computes a biorthogonal pair of basis matrices (Vk , Wk ) and an unreduced tridiagonal matrix Tk such that AVk = Vk Tk + βk+1 vk+1 e∗k ,

(3)

A∗ Wk = Wk Tk∗ + γ¯k+1 wk+1 e∗k

(4)

∗ Wk+1 Vk+1 = I

(5)

with

and

α

1

 β2  Tk =   

γ2 α2 .. .

 γ3 .. .

..

.

  .  

(6)

αk−1 γk βk αk A tridiagonal matrix is called unreduced if all its off-diagonal elements are nonzero. Equations (3) and (4) are called the Lanczos relations. In an eigenvalue context, they can be seen as approximate Schur decompositions of A and A∗ , from which approximate eigenvalues and eigenvectors can be computed as follows. A pair (θi , yi ) is an approximate (right) eigenpair of A if yi ≡ Vk ui

and

βk−1

Tk ui − θi ui = 0,

i = 1, . . . , k.

(7)

34

G. De Samblanx, A. Bultheel / Nested Lanczos

We assume that there exist k linear independent vectors yi . The corresponding (right) residual norm can be computed as kAyi − θi yi k2 = |ui,k βk+1 | kvk+1 k2 .

(8)

The algorithm has converged if |ui,k βk+1 |kvk+1 k2 ' 0. Analogously, a left eigenpair (θ¯i , yi(∗) ) can be defined by yi(∗) ≡ Wk u(∗) i such that

and

¯ (∗) Tk∗ u(∗) i − θi ui = 0,

i = 1, . . . , k,



∗ (∗)

A y − θ¯i y (∗) = u(∗) γ¯k+1 kwk+1 k2 . i i i,k 2

(9)

(10)

The residuals of the approximate eigenpairs fulfill the following Galerkin conditions Ayi − θi yi ⊥ R(Wk )

and

A∗ yi(∗) − θ¯i yi(∗) ⊥ R(Vk ),

(11)

which can be used as a definition for (θi , yi ). Other definitions of approximate eigenvalues and eigenvectors are possible, e.g., analogous to the concept of Harmonic Ritz values [12]. Algorithm 1 (Simple Lanczos). In: A ∈ Cn×n , v1 , w1 ∈ Cn , with kv1 k2 = kw1 k2 = 1, k ∈ N. Out: Vk , Wk ∈ Cn×k , αi , βi , γi for i = 1, . . . , k. 1. Set V1 = [v1 ], W1 = [w1 ], β1 = 0, γ1 = 0. 2. For j = 1, 2, . . . , k − 1 do 2.1. Set v+ = Avj and w+ = A∗ wj . 2.2. Compute αj = wj∗ v+ . 2.3. Biorthogonalize v+ = v+ − αj vj − βj vj−1 , ¯ j wj − γ¯j wj−1. w+ = w+ − α ∗ v 6= 0 then normalize 2.4. If w+ +

vj+1 = v+ /βj+1 , γj+1 , wj+1 = w+ /¯ ∗v . with βj+1 γj+1 = w+ +

Else Breakdown and Stop.

G. De Samblanx, A. Bultheel / Nested Lanczos

35

It is well known that the Lanczos algorithm can suffer from breakdown. Break∗ v = 0 (step 2.4 in algporithm 1), such that no β down occurs whenever w+ + k+1 and γk+1 can be found that normalize the new basis vectors vk+1 and wk+1 . Breakdown can have different causes, which can be classified as follows. Lucky breakdown. If v+ = 0 or w+ = 0, then an invariant subspace is found for A or A∗ . E.g., if v+ = 0, then the approximate eigenvectors in R(Vk ) are true eigenvectors for the problem – supposing that (3) holds exactly. The approximate eigenvalues θi are also equal to eigenvalues of A. If the wanted eigenvalues are among these eigenvalues, then the algorithm may stop. Otherwise, we can set βk+1 ≡ 0 and set vk+1 equal to any vector that fulfills the biorthogonality condition (5) and proceed with the algorithm. At this point, the matrix Tk+1 is no longer unreduced since it contains a zero subdiagonal element. ∗ v = 0, then the Lanczos algoSerious breakdown. If v+ 6= 0 and w+ 6= 0 but w+ + rithm cannot be continued any more. This event will occur if the moment matrix Mk is singular for some k.

Definition 2. Given a matrix A ∈ Cn×n and two vectors v1 , w1 ∈ Cn , the moment matrix Mk (A, v1 , w1 ) is defined as 
k−1 ∗   w1 v1 , Av1 , . . . , Ak−1 v1 . (12) Mk (A, v1 , w1 ) = w1 , A∗ w1 , . . . , A∗ A submatrix of a moment matrix is denoted by its size k and a “shift” parameter l as Mkl (A, v1 , w1 ) = Mk (A, Ai v1 , (A∗ )j w1 ), with i + j = l. If the parameters in this definition are obvious, then we omit them: Mk = Mk (A, v1 , w1 ) and Mkl = Mkl (A, v1 , w1 ). The matrix Mk (A, v1 , w1 ) is a Hankel matrix, since its (i, j)th element w1∗ Ai+j−2 v1 is equal to its (i−m, j +m)th element. If serious breakdown occurs, then the algorithm must be continued in a different way. One can restart the method explicitly by choosing new starting vectors v1 and w1 . But this is expensive and is no guarantee to avoid a new breakdown. A better solution is to give up the tridiagonal form of Tk and perform a look-ahead strategy, e.g., [2,6,14]. Roughly speaking, a look-ahead strategy tries to find a p > 0 such that Mk+p has full rank. After the look-ahead step, the algorithm continues with a block tridiagonal matrix Tk+p . We will show that an implicitly restarted Lanczos algorithm offers an alternative for look-ahead. Incurable breakdown. If none of the Mk , Mk+1 , . . . , Mn has full rank, then the breakdown is called incurable, since it cannot be cured with look-ahead. The only solution then is to restart the algorithm with a new pair of starting vectors. One can prove that if incurable breakdown occurs, then the eigenvalues of Tk are correct eigenvalues of A [18], as with lucky breakdown. However, the approximate eigenvectors are no true eigenvectors.

36

G. De Samblanx, A. Bultheel / Nested Lanczos

Summarizing, we can say that if breakdown occurs, then the algorithm must be restarted or the matrix Tk can no longer be an unreduced tridiagonal matrix. 3.

Implicitly Restarted Lanczos

Let us now show how the Lanczos algorithm can be restarted implicitly. With the term “implicitly restarting”, we denote the reorganization of equations (3) and (4) such that the new relations are correct Lanczos relations. This means that in exact arithmetic, they could have been generated by the Lanczos algorithm if we restarted it with the proper starting vectors. Since the new matrices V + and W + must be biorthogonal, the reorganization must consist of a biorthogonal transformation. Two algorithms can be used in order to compute this transformation: the Lanczos algorithm itself or the two-sided Gram–Schmidt algorithm [14, p. 107]. Theorem 3. Given the Lanczos relations (3) and (4) and suppose that there exist a pair of biorthogonal upper Hessenberg matrices Q, Z ∈ Ck×k and a tridiagonal matrix H ∈ Ck×k such that, for some µ ∈ C, (Tk − µI)Q = QH

and

(Tk − µI)∗ Z = ZH ∗ ,

(13)

then the matrices + Vk−1 ≡ Vk QIk,k−1,

+ Wk−1 ≡ Wk ZIk,k−1

and

+ ∗ ≡ Ik,k−1 (H + µI)Ik,k−1 Tk−1

(14)

define a restarted Lanczos relation + + + AVk−1 = Vk−1 Tk−1 + βk+ vk+ e∗k−1 ,
∗ + + + = Wk−1 + γ¯k+ wk+ e∗k−1 , Tk−1 A∗ Wk−1

(16)

βk+ vk+ = hk,k−1 Vk Qek + βk+1 qk,k−1vk+1 , ¯ k−1,k Wk Zek + γ¯k+1 zk,k−1wk+1 . γ¯k+ wk+ = h

(17) (18)

(15)

with

Proof. We prove the theorem only for the first Lanczos relation (3). The proof for (4) is analogous. If we shift (3) with µ and multiply it on the right by Q, then we get, recalling that Q is Hessenberg, (A − µI)Vk Q = Vk (Tk − µI)Q + βk+1 vk+1e∗k Q
= Vk (Tk − µI)Q + βk+1 vk+1 qk,k−1e∗k−1 + qk,k e∗k .

(19) (20)

If we insert I = QZ ∗, then we find


(A − µI)(Vk Q) = (Vk Q) Z ∗ (Tk − µI)Q + [0, . . . , 0, βk+1 qk,k−1vk+1 , βk+1 qk,k vk+1 ],

(21)

G. De Samblanx, A. Bultheel / Nested Lanczos

37

where Z ∗ (Tk − µI)Q = H. If we now remove the last column of this equation and shift it back, then we find, using (14), + AVk−1 = (Vk Q)(H + µI)Ik,k−1 + βk+1 qk,k−1vk+1 e∗k−1 .

(22)

∗ Inserting Ik = Ik,k−1Ik,k−1 + ek e∗k , we finally derive + ∗ = (Vk Q)Ik,k−1Ik,k−1 (H + µI)Ik,k−1 AVk−1

+ (Vk Q)ek e∗k (H + µI)Ik,k−1 + + Tk−1 + (hk,k−1 Vk Qek + = Vk−1

+ βk+1 qk,k−1vk+1 e∗k−1 βk+1 qk,k−1vk+1 )e∗k−1 .

(23) (24) (25)

This corresponds to (15), where βk+ vk+ = hk,k−1Vk Qek + βk+1 qk,k−1vk+1. The new basis matrices are biorthogonal, since
∗ + + ∗ Wk−1 Vk−1 = Ik,k−1 Z ∗ Wk∗ Vk QIk,k−1 ∗ ∗ = Ik,k−1 Z ∗ QIk,k−1 = Ik,k−1 Ik,k−1 = Ik−1,k−1. + ∗ + + Also (Wk−1 ) vk = 0, unless βk+ vk+ = 0, because if we multiply (15) by wk−1 and + ∗ + + combine this with (Wk−1 ) AVk−1 = Tk−1 , then the result follows. Notice that since vk+1 ∈ / R(Vk ), we expect that βk+ vk+ 6= 0. 

The matrix pair (Q, Z) can be computed by a full run of the Lanczos algorithm applied to the small, tridiagonal matrix Tk . Indeed, if we perform Lanczos on a tridiagonal matrix, using starting vectors that are non-zero only in their first two entries, then the resulting matrices will be upper Hessenberg and biorthogonal. Unless the Lanczos process on the tridiagonal matrix Tk breaks down. Lemma 4. Given a tridiagonal matrix T ∈ Ck×k and two vectors q1 , z1 ∈ C such that q1∗ z1 = 1. If q1 and z1 can be written as q1 = q1,1 e1 + q1,2 e2

and

z1 = z1,1 e1 + z1,2 e2 ,

(26)

then the matrices Q and Z that are generated from q1 , z1 and T using a Lanczos algorithm will be upper Hessenberg. Moreover, if z1 ≡ z1,1 e1 , then Z will be upper triangular and Q lower bidiagonal. Proof. The proof can be easily done by induction. If qi ∈ R(Ii+1,i ) then T qi ∈ R(Ii+2,i+1 ), so each new column of Q will have an additional nonzero entry. Therefore, Q must be upper Hessenberg. The same holds for Z. If z1 = z1,1 e1 , then clearly R(ZIk,i ) = R(Ik,i ), so the strict upper triangular part of Q must be zero in order to fulfill Z ∗ Q = I. 

38

G. De Samblanx, A. Bultheel / Nested Lanczos

There is a second possibility to compute Q and Z by using the two-sided Gram– Schmidt algorithm (BioGS) [14, p. 107]. In general, BioGS computes two biorthogonal basis matrices Q and Z for the column range of a set of matrices F , G ∈ Ck×l , i.e., Z ∗ Q = Il ,

R(F ) ⊂ R(Q)

and

R(G) ⊂ R(Z).

(27)

The BioGS algorithm can be seen as an algorithm that computes a “biorthogonal QRfactorization”, since it also produces two upper triangular matrices RQ and RZ ∈ Cl×l such that F = QRQ

and

G = ZRZ .

(28)

If F and G have full rank, then RQ and RZ will be invertible. An implementation of the two-sided Gram–Schmidt algorithm is given by algorithm 5. Algorithm 5 (BioGS). In: F , G ∈ Ck×l (k > l). Out: Q, Z ∈ Ck×l , RQ , RZ ∈ Cl×l . 1. Set Q, Z, RQ , RZ to [ ]. 2. For j = 1, . . . , k do 2.1. Set fj = F ej and gj = Gej . 2.2. Compute rQ,j = Z ∗ fj , qZ,j = Q∗ gj . 2.3. Biorthogonalize fj = fj − QrQ,j , 2.4. If

gj∗ fj

gj = gj − ZrZ,j .

6= 0 then normalize

qj = fj /ρQ ,

zj = gj /ρZ with ρ¯Z ρQ = gj∗ fj .

Else Breakdown and Stop. 2.5. Set



RQ RQ = 0

 rQ,j , ρQ



RZ RZ = 0

 rZ,j . ρZ

The BioGS algorithm can break down, even if rank(G∗ F ) = l. We assume that for our application, this will not occur. We use the BioGS algorithm with F = Tk −µI ¯I, in order to apply theorem 3. and G = Tk∗ − µ Lemma 6. Suppose that T ∈ Ck×k is an unreduced tridiagonal matrix, and F = T −αI and G = (T −βI)∗ which are decomposed as in (28) with the BioGS algorithm. If Vk and Wk are the results of k Lanczos steps (without breakdown) applied to T with v1 = F e1 and w1 = Ge1 , then there exists a diagonal scaling matrix D such that Q = Vk D and Z ∗ = D−1 Wk∗ .

G. De Samblanx, A. Bultheel / Nested Lanczos

39

Proof. It is clear that there exist a ξ and a ζ such that q1 = ξ(T − αI)e1 = ξv1 and z1 = ζ(T − βI)∗ e1 = ζw1 . Suppose that T Vk = Vk H and T ∗ Wk = Wk H ∗ . Since Vk , Wk ∈ Ck×k are unreduced upper Hessenberg (apply lemma 4 on a Lanczos run without breakdown), [ e1 , Vk−1 ] and [ e1 , Wk−1 ] are full rank upper triangular. If we combine this with the Lanczos relation T Vk−1 = Vk HIk,k−1, then     (T − αI) e1 , Vk−1 = Vk e1 , (H − αI)Ik,k−1   −1
⇒ F = (T − αI) = Vk e1 , (H − αI)Ik,k−1 e1 , Vk−1 = Vk RV . Analogously, G = (T − βI)∗ = Wk



e1 , (H − βI)Ik,k−1



e1 , Wk−1

−1


= Wk RW .

Up to scaling of the matrices Vk and Wk , this corresponds to (28).



Notice that applying BioGS on a differently shifted version of T and T ∗ , i.e., with α 6= β, still corresponds to the Lanczos algorithm applied on the matrix T . The values of α and β are included in the starting vectors of the inner Lanczos loop. The restarting procedure of theorem 3 requires k − 1 steps of the small Lanczos algorithm, i.e., the Lanczos algorithm applied to the small matrix Tk − µI. Only the first k − 1 columns of Q and Z are used, along with the last off-diagonal elements of H, which are also computed at step k − 1. The following lemma sets a condition for the algorithm that computes Q and Z under which it will not break down. Lemma 7. Consider theorem 3. If the matrices Q and Z are computed using the Lanczos (or BioGS) algorithm, then this algorithm will not break down unless for some j < k, Wj∗(A − µI)2 Vj does not have full rank. Proof. The algorithm breaks down at step j if Mj0 ≡ ((Tk − µI)∗ Ik,j )∗ (Tk − µI)Ik,j is singular. Since this matrix is equal to Mj0 = (Wk (Tk − µI)∗ Ik,j )∗ Vk (Tk − µI)Ik,j , we can insert (3) and (4), so
∗ Mj0 = (A − µI)∗ Wk Ik,j (A − µI)Vk Ik,j = Wj∗ (A − µI)2 Vj .  The result of lemma 7 can be interpreted in terms of the implicitly restarted process. Since a Krylov subspace is invariant when A is shifted, it holds that R(Vk ) = Kk (A, v1 ) = Kk (A − µI, v1 ) and a similar relation holds for R(Wk ). Lemma 8. Given A, v1 , w1 as in definition 2, then rank(Mk (A, v1 , w1 )) = rank(Mk (A − µI, v1 , w1 )). Proof.

The lemma follows from the fact that Kk (A − µI, v1 ) ⊂ Kk (A, v1 ), since (A − µI)m v1 =

m X i=0

m! µm−i Ai v1 ∈ Kk (A, v1 ), (m − i)!i!

40

G. De Samblanx, A. Bultheel / Nested Lanczos

and for the same reason


Kk (A, v1 ) ⊂ Kk (A − µI) + µI, v1 = Kk (A − µI, v1 ).

Analogously, Kk (A∗ , w1 ) = Kk ((A − µI)∗ , w1 ). rank(Mk (A − µI, v1 , w1 )).

Hence, rank(Mk (A, v1 , w1 )) = 

Therefore, Mk0 = Wk∗ (A − µI)2 Vk corresponds to the moment matrix Mk (A, (A − µI)v1 , (A − µI)w1 ) = Mk2 ((A − µI), v1 , w1 ) of the Lanczos process with starting vectors v1+ = (A − µI)v1 and w1+ = (A − µI)∗ w1 . Lemma 7 shows that the small Lanczos iteration will not break down if the “large” Lanczos process with (v1+ , w1+ ) does not break down. In the following section, we show in theorem 10 that this is exactly the implicit Lanczos run. So the implicit restart does not break down if the corresponding explicit restart does not break down, which is not very surprising. An analogous result was found by Grimme et al. [7, theorem 3] for the sign-symmetric case. Algorithm 9 gives an implementation for Implicitly Restarted Lanczos. Algorithm 9 (Implicitly Restarted Lanczos (IRL)). In: Vk+1 , Wk+1 , Tk , βk+1 , γk+1 , µQ , µZ . + , βk+ , γk+ . Out: Vk+ , Wk+ , Tk−1 1. Compute Q, Z and RQ , RZ ∈ Ck×k as Tk − µQ I = QRQ and (Tk − µZ I)∗ = ZRZ . + + + = Vk QIk,k−1, Wk−1 = Wk ZIk,k−1 and Tk−1 = Ik−1,k Z ∗ Tk QIk,k−1. 2. Set Vk−1

3. Compute v+ = hk,k−1Vk Qek + βk+1 qk,k−1vk+1 and w+ = ¯hk−1,k Wk Zek + γ¯k+1 zk,k−1wk+1 . ∗v . 4. Set vk+ = v+ /βk+ and wk+ = w+ /γk+ with βk+ γ¯k+ = w+ +

4.

Applying an implicit filtering step

The choice of the shift µ in theorem 3 or α and β in lemma 6 is crucial for an efficient restart of the Lanczos algorithm. Indeed, if µ would be chosen such that most of the “wanted” information is removed from the subspaces Vk and Wk , then the algorithm would slow down, since this information must be recovered. On the other hand, µ can be used in order to remove uninteresting information, e.g., eigenvectors that correspond to unwanted eigenvalues. As for IRA [16], we prove two properties for the restart. First, we show that an implicit restart that is performed with the BioGS algorithm implicitly applies a polynomial filter on Vk and Wk , as if it performed a step of subspace iteration. Then we show that the concept of exact shifts (a notion

G. De Samblanx, A. Bultheel / Nested Lanczos

41

that was first introduced for IRA) also holds for the Lanczos case: if µ is equal to an approximate eigenvalue, then this approximate eigenvalue will be removed from the spectrum of Tk . The other eigenvalues remain unaltered. Theorem 10 (Implicit filtering property). Suppose that in theorem 3, Q and Z are computed using the BioGS algorithm on F = Tk − µQ I and G = (Tk − µZ I)∗ , such that (28) holds. If RQ and RZ have full rank, then

R Vk+ = R (A − µQ I)Vk , (29)

+ ∗ R Wk = R (A − µZ I) Wk . (30) Proof.

Let us prove the first statement (29). From (3) and (28), we derive

(A − µQ I)Vk = Vk (Tk − µQ I) + βk+1 vk+1 e∗k = Vk QRQ + βk+1 vk+1 e∗k .

(31)

−1 Multiplying this on the right by RQ results in −1 = Vk Q + βk+1 /ρQ vk+1e∗k , (A − µQ I)Vk RQ

(32)

where ρQ is the (k, k)th element of RQ . Since the tridiagonal matrix H = Z ∗ Tk Q is factorized as H = RQ Q + µQ I (28), it holds that hk,k−1 = ρQ qk,k−1, such that by (17) −1 = [Vk QIk,k−1, Vk Qek + βk+1 qk,k−1/hk,k−1 ] (A − µQ I)Vk RQ  +  = Vk−1 , βk+ /hk,k−1vk+ ,

from which (29) follows. The proof of (30) is analogous.

(33) (34) 

Theorem 10 states that if the Implicitly Restarted Lanczos algorithm is repeated p times with shifts µQ,1 , . . . , µQ,p and µZ,1 , . . . , µZ,p using BioGS, then the subspace spanned by Vk−p+1 is multiplied implicitly by a polynomial function φ(A) = (A − µQ,pI) · · · (A − µQ,1 I). Similarly, Wk−p+1 is multiplied by ¯Z,p I) · · · (A∗ − µ ¯Z,1I). It should be noted that if this algorithm is φ(∗) (A∗ ) = (A∗ − µ applied on a symmetric matrix A and v1 = w1 , such that Wk = Vk , then the filtering with µi ≡ µQ,i = µZ,i is still different from an IRA filtering. With IRA, the filtering function for Vk is φ(A), so the same holds for Wk – since both matrices are equal. However, unless the µi are real, φ(A) 6= φ(∗) (A∗ ) = φ(∗) (A). It can be shown for IRA [16, lemma 3.10] that if the approximated spectrum is divided into two disjoint parts {θ1 , . . . , θp } ∪ {θp+1, . . . , θk }

(35)

and if the first set of p approximate eigenvalues is used as shifts µ1 ≡ θ1 , . . . , µp ≡ θp , then the approximate eigenvalues of the restarted relation will be given by {θp+1, . . . , θk }. Also the approximate eigenvectors (left and right) will be the same. The first p approximate eigenvectors are said to be deflated.

42

G. De Samblanx, A. Bultheel / Nested Lanczos

Before we derive a similar property for the IRL algorithm, we propose a simple condition for deflation. Suppose that the subspaces (Vk , Wk ) are restarted (or trans+ + , Wk−1 ). We say that the ith approximate eigenpair (θi , yi ) is not formed) into (Vk−1 deflated if the vector that is removed from Vk is orthogonal to yi . Lemma 11. Say that the right approximate eigenpairs of A are given by (θi , yi = + , v]) = Vk ui ), i = 1, . . . , k, where θi 6= 0. If w, v ∈ Cn are vectors such that R([Vk−1 + + + R(Vk ), R([Wk−1 , w]) = R(Wk ) and [Wk−1 , w]∗ [Vk−1 , v] = I, then w ⊥ yi ⇔ there k−1 such that T + u+ − θ u+ = 0 and y = V + u+ . exists a u+ i i i i ∈C k−1 i k−1 i In other words, w is orthogonal to yi iff (θi , yi ) is an approximate eigenpair of the restarted relation. The same holds for the left approximate eigenpairs. + , v]PV Proof. There exist transformation matrices PV and PW such that Vk = [Vk−1 + −1 ∗ ∗ and Wk = [Wk−1 , w]PW ; clearly PW PV = I. If we set u ˜i = PV ui = PW ui , then + yi = [Vk−1 , v]˜ ui and

 + ∗  +  ∗ ∗ Tk ui = θi PW ui ⇒ Wk−1 , w A Vk−1 ,v u ˜ i = θi u ˜i . PW

(36)

+ + + , w]∗ [Vk−1 , v] = I, we get w∗ [Vk−1 , v] = e∗k . Thus, by multiplying Since [Wk−1 with u ˜i ,   + ,v u ˜i = e∗k u ˜i . w∗ yi = w∗ Vk−1

If w∗ yi = 0, then this means that we can rewrite  u+ i , u ˜i = 0 

(37)

+ + + ∗ + and yi = Vk−1 u+ i . Because Tk−1 = (Wk−1 ) AVk−1 , it then follows by (36) that + + Tk−1 u+ ˜i fulfills (37), then it must follow that w∗ yi = 0, i − θi ui = 0. Inversely, if u  so w ⊥ yi .

We now show that if the Lanczos relations are restarted with BioGS and µ = θj , + ) is spanned by the remaining approximate eigenthen the restarted subspace R(Vk−1 vectors yi , i 6= j, because it fulfills the conditions of lemma 11. Lemma 12 (Exact shift property). Say that the approximate eigenpairs of A are given by (θi , yi = Vk ui ), i = 1, . . . , k. If theorem 3 is applied with BioGS and µ ≡ θj , + , then 1 6 j 6 k, then the following statement is true: If w ∈ R(Wk ) and w ⊥ Vk−1 + w ⊥ yi , for i such that θi 6= θj . Hence, the restarted subspace Vk−1 contains the same Ritz vectors as Vk , except for one vector that corresponds to θj .

G. De Samblanx, A. Bultheel / Nested Lanczos

43

Proof. If θj is an eigenvalue of Tk , then (Tk − θj I) = QRQ does not have full rank and the (k, k)th element of RQ is zero (if another diagonal element were zero, the ∗ RQ and thus BioGS would have had a breakdown). Hence, RQ = Ik,k−1Ik,k−1 + ∗ ∗ RQ = Vk−1 Ik,k−1 RQ . Vk QRQ = Vk QIk,k−1Ik,k−1

Following theorem 10, it holds that (see (32)) w∗ (A − θj I)Vk ui = w∗ Vk QRQ ui + βk+1 w∗ vk+1 e∗k ui + ∗ = w∗ Vk−1 Ik,k−1 RQ ui + βk+1 w∗ vk+1 e∗k ui = βk+1 w∗ vk+1 e∗k ui .

On the other hand, (31) implies w∗ (A − θi I)Vk ui = w∗ Vk (Tk − θi I)ui + βk+1 w∗ vk+1 e∗k ui = βk+1 w∗ vk+1 e∗k ui . (38) If we subtract both equations, then we find that (θj − θi )w∗ Vk ui = (θj − θi )w∗ yi = 0. Using lemma 11, we get the result.  A similar result can be formulated for Wk and for the left approximate eigenpairs, if µZ = θj . If a pair of Lanczos relations is restarted with BioGS and with a shift equal to the approximate eigenvalue, then the approximate eigenvector is filtered out of the subspace basis. By analogy to IRA, we call this choice of shifts the use of exact shifts. If the Lanczos relations are restarted with p exact shifts, then the corresponding eigenpairs are removed from the relations. Corollary 13. Suppose that the spectrum of Tk is divided into two disjoint parts {θ1 , . . . , θp } ∪ {θp+1, . . . , θk }

(39)

and suppose that the first set of p approximate eigenvalues is used as exact shifts for the implicitly restarting procedure. Then the right approximate eigenpairs of the resulting equations are given by (θi+ , yi+ ) = (θi+p , yi+p ), i = 1, . . . , k − p. The same holds for the left approximate eigenpairs. Proof.

The proof follows from lemmas 11 and 12, repeated p times.



Example 14. Let us illustrate these results with a small example. Consider a Toeplitz matrix A ∈ R100×100 that is defined by 1 1 for k > 0 and ai,i−k = − for k > 0. k−1 k−1 All eigenvalues λi of A have real part equal to one: real(λi ) = 1, i = 1, . . . , 100. We performed 10 steps of Lanczos with v1 = e1 = w1 and then restarted with BioGS and three different choices of µ = µQ = µZ : a true eigenvalue µ = λ1 , an approximate eigenvalue µ = θ1 and a random number µ = 1.5. The results are displayed in table 1. The table illustrates theorem 10 and lemma 12 to computer precision. If BioGS is ai,i+k =

44

G. De Samblanx, A. Bultheel / Nested Lanczos Table 1 Restarting the Lanczos algorithm using BioGS. For three different choices of µ (a true eigenvalue, an approximate eigenvalue and a “random” number), the error on the implicit filter is shown (PV⊥+ (A − µI)Vk is the projection of (A − µI)Vk on the nullspace of Vk+ ). Also the difference between the old and the new (second) approximate eigenvalue is shown, as well as the error on the Lanczos equation.

µ = λ1 µ = θ1 µ = 1.5

kPV⊥+ (A − µI)Vk k2

|θ2 − θ2+ |

+ + + kAVk−1 − Vk−1 Tk−1 − βk+ vk+ k2

8e−15 3e−15 2e−14

3e−3 9e−16 3e−1

3e−15 3e−15 2e−14

used, then the filtering property holds for all three cases. If an exact shift is used, then the corresponding approximate eigenvalue is removed from the spectrum of T + , without changing the other approximate eigenvalues. We also used an inner Lanczos loop to restart the algorithm. The results for the error on the Lanczos relation were the same, but neither the implicit filtering property nor the exact shift property was fulfilled, as expected.

5.

Implicitly restarting and breakdown

The Lanczos algorithm can break down. If this breakdown is not the result of the fact that some solution has been found, then it must be solved with look-ahead or it must be restarted completely. The main argument against an explicit restart is that it is expensive. It costs 2k matrix–vector multiplications and as many biorthogonalization steps. Also, it cannot guarantee that a breakdown will not reoccur in further iterations, although the same argument holds for look-ahead. In this section, we show how the implicit restart algorithm can be used to get around the serious breakdown problem. Theorem 15. Suppose that the Lanczos relations are restarted as in theorem 3, and + + + ∗ + ∗ v suppose that wk+1 k+1 = 0. Then βk γk (wk ) vk = hk,k−1 hk−1,k . Proof. Multiplying (17) by (18) results in
∗ βk+ γk+ wk+ vk+

¯ k−1,k Wk Zek ∗ (hk,k−1 Vk Qek ) + γ¯k+1zk,k−1wk+1 ∗ (βk+1 qk,k−1vk+1 ) = h ∗ vk+1 = hk−1,k hk,k−1e∗k Z ∗ Wk∗ Vk Qek + γk+1 βk+1 z¯k,k−1 qk,k−1wk+1

= hk−1,k hk,k−1 + 0.



Theorem 15 shows that if H is unreduced, then no breakdown occurs in the next Lanczos step of the restarted method. If hk,k−1 = 0, then the last column of Q turns out to be an eigenvector of Tk . In this case, a new decomposition (13) must be

G. De Samblanx, A. Bultheel / Nested Lanczos

45

found. Notice that this problem occurs with BioGS if the shift µ is chosen equal to an approximate eigenvalue θi . In that case, Tk − µI is singular and the last vector of Q will be the null-eigenvector (it must be the last vector, because otherwise the BioGS algorithm breaks down). The result of theorem 15 is at least against our intuition. The restarting procedure applies a polynomial filter on the subspaces Vk and Wk , but the new Krylov subspaces will be subsets of the old subspaces, even when the algorithm proceeds. One could argue that, since breakdown seems to be a property of these Krylov subspaces and since by restarting we remain in the same subspaces, breakdown cannot be avoided with an implicit restart. In general, this is not true. The new search spaces are indeed very similar to the Krylov subspaces that would have been found if the algorithm proceeded, but each new pair of vectors is biorthogonalized to one vector less (i.e., Vk Qek and Wk Zek ). It is this “missing” vector that prevents the new, biorthogonalized vectors to be orthogonal to each other. Indeed, following (17) and (18) the vectors vk+ and wk+ are combinations of vk+1 or wk+1 and of Vk Qek or Wk Zek . The term hk,k−1 hk−1,k corresponds to the parts of vk+ and wk+ that contain traces of Vk Qek and Wk Zek , which are not orthogonal. There is a relation between the moment matrix of restarted Lanczos relations and the moment matrix of the original relations. Suppose that in theorem 10, µ ≡ µQ = µZ , then the moment matrix of the restarted relation is Mk (A, (A − µI)v1 , (A − µI)∗ w1 ). The rank of this matrix is equal to


rank Mk A − µI, (A − µI)v1 , (A − µI)∗ w1 = rank Mk2 (A − µI, v1 , w1 ) . Indeed, the lower right k × k submatrix of Mk+1 (A − µI, v1 , w1 ) is equal to the lower left submatrix of Mk+2 (A − µI, v1 , w1 ), because of the Hankel structure (see figure 1). Therefore, implicitly restarting corresponds to shifting to a submatrix of the original moment matrix. In the following lemma, we show that if Mk+1 , . . . , Mk+p are singular, then it takes at least p restarts to avoid the singularity. We only consider restarts with µ = µQ = µZ = 0, because it is easier to understand. Since the rank of the moment matrix is equal for A and A − µI (lemma 8), it also holds for the general case. Lemma 16. Consider a matrix A and two starting vectors v1 , w1 . Suppose that all implicit restarts are computed with BioGS and with µ = 0. If rank(Mk ) = rank(Mk+1 ) = · · · = rank(Mk+p ) = k, then: (a) If l is such that rank(Mk+l ) = k + l, then l > p. If i 6 p, then rank(Mk+p+i ) 6 rank(Mk+p+i−1 ) + 2. 2j (b) If rank(Mk+p+1 ) > k + 1, then for j < p, we have rank(Mk+1 ) < k + 1 and 2j−1 2p rank(Mk+1 ) < k + 1. Moreover, rank(Mk+1 ) = k + 1, unless a unit vector ei ∈ R(Mk+2p+1 Ik+2p+1,2p ), for some 1 6 i 6 2p.

Proof.

Consider the rectangular Hankel matrix Mk,l ≡ Mk Ik,l , k > l.

46

G. De Samblanx, A. Bultheel / Nested Lanczos

1

1

M6 1 1 M1 4

1

1

1

*   1 1 1  M22 1     1 1 M24 1 @ R @

1 6

1

1

2

5

1

2

1

4

2

1

1

3

1

1

1

M2

1

1

2

rank(Mk )

6

c

c

2 c

1 1

2

1

c

1

1

c

c -k

1 1

2

3

4

5

6

Figure 1. A connection between IRL and look-ahead Lanczos. The look-ahead Lanczos algorithm jumps from the (full rank) moment matrices M1 to M6 . Implicitly restarted Lanczos results with each restart in a shifted moment matrix M2 → M22 → M24 , until the breakdown has disappeared. The circles in the right picture show the rank of the Mk matrices. The rank does not change from M1 to M3 and then grows with steps of one or two.

(a) Since Mk+1(A, v1 , w1 ) differs from Mk (A, v1 , w1 ) in one row and one column, the rank of Mk+1 can be at most rank(Mk ) + 2. If the nullspace of Mk+p has dimension p, then it will take at least p new rows and p new columns in order to end up with a full rank moment matrix. (b) In [8, corollary 5.1, p. 81], it is proven that  (40) rank(Mk+2l,k ) = min k, rank(Mk+l )   (41) rank(Mk+2l+1,k ) = min k, rank(Mk+l+1,k+l ) 6 min k, rank(Mk+l ) . 2i is the lower k + 1 × k + 1 submatrix of M Since Mk+1 k+2i+1,k+1 , we can derive that 
2i rank Mk+1 6 rank(Mk+2i+1,k+1 ) = min k + 1, rank(Mk+i+1 ) < k + 1. 2i−1 For the same reason rank(Mk+1 ) 6 min{k + 1, rank(Mk+i+1 )}. On the other hand,  2j−1
rank Mk+1 6 min k + 1, rank(Mk+p+1 ) = k + 1.

If ei ∈ / R(Mk+2p+1,k+1) for i = 1, . . . , 2p, then 2p
rank Mk+1 = rank(Mk+2p+1,k+1 ) = k + 1, since



[Mk+2p+1,k+1 then has full rank.

? Ik+2p+1,2p ] = 2p Mk+1

 I2p , O

G. De Samblanx, A. Bultheel / Nested Lanczos

47

Due to the correspondence of Mk2 to Mk (A, v1+ , w1+ ), applied recursively, the Lanczos algorithm will avoid the singularity by restarting p times.  Lemma 16 proves two properties of a singularity in the Lanczos process. The first part shows that a singularity in a sequence of Mk matrices consists of a part where the rank is constant and a second part where the rank in general grows with 2 (unless there is an overlap of singularities or if the size of the singularity is odd). The second part proves that if a look-ahead algorithm has to compute a look-ahead step of length 2p or 2p − 1, then the corresponding implicitly restarted Lanczos method has to be restarted at least p times. The connection of implicitly restarted Lanczos with look-ahead is illustrated in figure 1. In the figure, Mk+1 = M2 is singular. It is easy to see that M2 , . . . , M5 are singular too, but rank(M6 ) = 6. A look-ahead algorithm has to take a step of length 5 in order to find a new pair of biorthogonal matrices. If IRL is used, then one iteration 2 = M22 is after the first restart, a new breakdown will be encountered, since Mk+1 4 singular. If the algorithm is then restarted again, the moment matrix M2 is no longer singular, and the algorithm can proceed without breakdown. So, look-ahead computes M2 → M3 → M4 → M5 → M6 , whereas IRL computes M2 → M22 → M24 . If a Lanczos algorithm encounters an incurable breakdown, so that there exists no lookahead step that solves it, then IRL cannot solve the breakdown either. Example 17. We generate a set of matrices and starting vectors that lead to breakdown of the Lanczos algorithm as follows. Say that breakdown must occur at step m and that look-ahead must compute a step of length p. Consider then a random unreduced tridiagonal matrix B. If we perform Lanczos with v1 = e1 = w1 , then Vk = In,k = Wk . If we set bm+1,m = 0 and bm+p+1,m = 1, then it is easy to see that breakdown must occur: vm = em+p+1 ⊥ wm = em . The matrix can then be transformed, using a biorthogonal pair of matrices Q and Z ∈ Cn×n : A = QBZ ∗,

v1 = Qe1

and

w1 = Ze1 .

We generated four 100 × 100 matrices which suffered from breakdown at step m = 5, each with a growing singularity (p = 1, 2, 3, 4). Then we applied the following algorithms to these matrices: implicitly restarted Lanczos with BioGS √ (µ = µQ = µZ = 0), IRL restarted with Lanczos (µ = 0 and q1 = z1 = (e1 + e2 )/ 2) and look-ahead Lanczos. The results are shown in table 2. The table illustrates the results of lemma 16. The implicitly restarted Lanczos algorithm needs dp/2e restarts in order to handle the singularity. We also showed the error on the first Lanczos equation after the restart (i.e., after 10 iterations). The accuracy of all three applications is comparable, but it certainly depends on our implementation and on the parameters that are used for the restart. The optimal choice of the parameters, if it exists, should be the subject of further research.

48

G. De Samblanx, A. Bultheel / Nested Lanczos

Table 2 Comparison between IRL and look-ahead Lanczos for breakdown of different lengths (p = 1, 2, 3, 4) at step 5 of the algorithm. Shown is the number of required restarts (l, for look-ahead-Lanczos this indicates the length of the look-ahead jump) and the overall error on the first Lanczos equation when the size of Vk and Wk is 10. The results are generated with MATLAB4 on a DEC5000. p=1 l IRL-BioGS 1 IRL-Lanczos 1 LA-Lanczos 1

p=2

p=3

p=4

kAV − V T k

l

kAV − V T k

l

kAV − V T k

l

kAV − V T k

2e−11 5e−11 1e−11

1 1 2

2e−11 2e−10 2e−12

2 2 3

1e−10 2e−10 6e−10

2 2 4

9e−10 4e−10 1e−10

Nested Lanczos vs. look-ahead. There is one key disadvantage to the use of Nested Lanczos in order to avoid breakdown: transforming the bases Vk and Wk implies the availability of these matrices in the computer memory. Let us explain this briefly. The main motivations for the choice of an unsymmetric Lanczos algorithm for the computation of eigenvalues with a three-term recurrence relation are that the orthogonalization cost is low and the basis vectors need not be stored in fast memory. Also, left eigenvectors are approximated as well as right eigenvectors. On the other hand, there does not seem to exist a numerically stable scheme for computing eigenvalues of a large unsymmetric tridiagonal matrix [3,4]. The advantage of Nested Lanczos is that the matrix Tk is known to be tridiagonal, so the storage requirements for it can be determined in advance. Moreover, it has the favorable properties of an implicitly restarting scheme, such as implicit filtering. But implicitly restarting requires the biorthogonal transformation of the basis vectors, so the storage advantage is compromised. A look-ahead solution does not have this problem: if the eigenvectors are not needed, then it may even omit to store Vk and Wk , except for the last two vectors. Then again, if reorthogonalization is required, there are more than these vectors needed. Also, if more than one eigenvalue is needed, then restarting the eigenvalue solver seems inevitable. Summarizing, if it is a problem to restart, then clearly look-ahead must be used. If the cost of transforming the previous basis vectors is not dominating, then restarting as a cure for breakdown is the option to choose. 6.

Conclusions

This text generalizes the concept of implicitly restarting an iterative algorithm to the Lanczos method for eigenvalue problems. We showed that a full run of the Lanczos algorithm on the small, projected eigenvalue problem generates a biorthogonal transformation for the Lanczos bases. The small Lanczos run can also be written as a BioGS algorithm, and connects the restarted subspaces with implicit filtering of the original subspaces. The polynomials that define this filtering are governed by the starting vectors of the small Lanczos process or by the shifts µQ and µZ of the BioGS factorization.

G. De Samblanx, A. Bultheel / Nested Lanczos

49

The resulting IRL algorithm can be used as an alternative for look-ahead in case of serious breakdown. IRL and look-ahead are connected in that if look-ahead extends the subspaces Vk and Wk with size p, then IRL selects in these large subspaces a set of k-dimensional subspaces Vk+ and Wk+ that do not suffer from breakdown. At the cost of dp/2e extra iteration steps, IRL returns a tridiagonalization of A, whereas lookahead Lanczos results in a block tridiagonal matrix. Both algorithms fail if incurable breakdown occurs. The advantages of implicitly restarted Lanczos correspond to the advantages of implicitly restarted Arnoldi with respect to Arnoldi. If the Lanczos algorithm converges too slowly, then IRL can reduce the size of the subspaces Vk and Wk without losing too much information. If there are spurious eigenvalues (e.g., approximations of an infinite eigenvalue) or if a known, dominant eigenvalue must be avoided (e.g., a zero eigenvalue), then IRL can be used to filter away these eigenvalues. Finally, it is an alternative for look-ahead in case of breakdown. It has the additional advantage that even in case of near-breakdown, IRL can be employed and the tridiagonal structure of the Tk matrix is preserved. So the decision “is this a breakdown or not?” is less important when IRL is used, since the result is tridiagonal anyway – look-ahead makes the matrix block-tridiagonal. The disadvantage of implicitly restarted Lanczos is that it costs one iteration step per restart, since IRL reduces the size of the approximating subspaces with one. It is also possible that the restarting algorithm itself incorporates some of the disadvantages of the Lanczos algorithm, such as possible instabilities or inaccurate results due to a near-breakdown. Acknowledgements This research was supported by the National Fund for Scientific Research (FWO), project “Orthogonal systems and their applications”, grant #G.0278.97. This paper presents research results of the Belgian Programme on Interuniversity Poles of Attraction, initiated by the Belgian State, Prime Minister’s Office for Science, Technology and Culture. The scientific responsibility rests with the authors. We also thank the referees for the many useful comments. References [1] Z. Bai, Error analysis of the Lanczos algorithm for the nonsymmetric eigenvalue problem, Math. Comp. 62 (1994) 209–226. [2] C. Brezinski, M. Redivo-Zaglia and H. Sadok, A breakdown-free Lanczos-type algorithm for solving linear systems, Numer. Math. 63(1) (1993) 29–38. [3] A. Bunse-Gerstner, An analysis of the HR algorithm for computing the eigenvalues of a matrix, Linear Algebra Appl. 35 (1981) 155–173. [4] A. Bunse-Gerstner, On the similarity transformation to tridiagonal form, Linear and Multilinear Algebra 12 (1982) 51–56. [5] D. Calvetti, L. Reichel and D.C. Sorensen, An implicitly restarted Lanczos method for large symmetric eigenvalue problems, ETNA 2 (1994) 1–21.

50

G. De Samblanx, A. Bultheel / Nested Lanczos

[6] R.W. Freund, M.H. Gutknecht and N.M. Nachtigal, An implementation of the look-ahead Lanczos algorithm for non-Hermitian matrices, Part I, SIAM J. Sci. Statist. Comput. 14(1) (1993) 137–158. [7] E.J. Grimme, D.C. Sorensen and P. Van Dooren, Model reduction of state space systems via an implicitly restarted Lanczos method, Numer. Algorithms 12 (1996) 1–32. [8] G. Heinig and K. Rost, Algebraic Methods for Toeplitz-like Matrices and Operators, Mathematical Research, Vol. 19 (Akademie-Verlag, Berlin, 1984) p. 212. [9] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, J. Res. Nat. Bur. Stand. 45 (1950) 255–281. [10] R.B. Lehoucq, D.C. Sorensen and P. Vu, ARPACK: An implementation of the Implicitly Re-started Arnoldi Iteration that computes some of the eigenvalues and eigenvectors of a large sparse matrix, available from [email protected] under the directory scalapack (1995). [11] R.B. Morgan, On restarting the Arnoldi method for large nonsymmetric eigenvalue problems, Math. Comp. 65 (1996) 1213–1230. [12] C. Paige, B.N. Parlett and H.A. Van der Vorst, Approximate solutions and eigenvalue bounds from Krylov subspaces, Numer. Linear Algebra Appl. 2 (1995) 115–133. [13] B.N. Parlett, Reduction to tridiagonal form and minimal realizations, SIAM J. Matrix Anal. Appl. 13(2) (1992) 567–593. [14] B.N. Parlett, D.R. Taylor and Z.A. Liu, A look-ahead Lanczos algorithm for unsymmetric matrices, Math. Comp. 44 (1985) 105–124. [15] Y. Saad, Variations on Arnoldi’s method for computing eigenelements of large unsymmetric matrices, Linear Algebra Appl. 34 (1980) 269–295. [16] D.C. Sorensen, Implicit application of polynomial filters in a k-step Arnoldi method, SIAM J. Matrix Anal. Appl. 13 (1992) 357–385. [17] A. Stathopoulos, Y. Saad and K. Wu, Dynamic thick restarting of the Davidson, and the Implicitly Restarted Arnoldi methods, SIAM J. Sci. Comput., to appear (1997). [18] D.R. Taylor, Analysis of the look-ahead Lanczos algorithm, Ph.D. thesis, University of California, Berkeley (1982).