• 沒有找到結果。

# 7.1 The Lanczos Algorithm

N/A
N/A
Protected

Academic year: 2022

Share "7.1 The Lanczos Algorithm"

Copied!
38
0
0

(1)

## Lanczos Methods

In this chapter we develop the Lanczos method, a technique that is applicable to large sparse, symmetric eigenproblems. The method involves tridiagonalizing the given matrix A. However, unlike the Householder approach, no intermediate (an full) submatrices are generated. Equally important, information about A0s extremal eigenvalues tends to emerge long before the tridiagonalization is complete. This makes the Lanczos algorithm particularly useful in situations where a few of A0s largest or smallest eigenvalues are desired.

### 7.1 The Lanczos Algorithm

Suppose A ∈ Rn×n is large, sparse and symmetric. There exists an orthogonal matrix Q, which transforms A to a tridiagonal matrix T .

QTAQ = T ≡ tridiagonal. (7.1.1)

Remark 7.1.1 (a) Such Q can be generated by Householder transformations or Givens rotations.

(b) Almost for all A (i.e. all eigenvalues are distinct) and almost for any q1 ∈ Rn with kq1k2 = 1, there exists an orthogonal matrix Q with first column q1 satisfying (7.1.1). q1 determines T uniquely up to the sign of the columns (that is, we can multiply each column with -1).

Let (x ∈ Rn)

K[x, A, m] = [x, Ax, A2x, · · · , Am−1x] ∈ Rn×m. (7.1.2) K[x, A, m] is called a Krylov-matrix. Let

K(x, A, m) = Range(K[x, A, m]) = Span(x, Ax, · · · , Am−1x). (7.1.3) K(x, A, m) is called the Krylov-subspace generated by K[x, A, m].

Remark 7.1.2 For each H ∈ Cn×m or Rn×m (m ≤ n) with rank(H) = m, there exists an Q ∈ Cn×m or Rn×m and an upper triangular R ∈ Cm×m or Rm×m with QQ = Im

such that

H = QR. (7.1.4)

(2)

Q is uniquely determined, if we require all rii> 0.

Theorem 7.1.1 Let A be symmetric (Hermitian), 1 ≤ m ≤ n be given and dimK(x, A, m) = m then

(a) If

K[x, A, m] = QmRm (7.1.5)

is an QR factorization, then QmAQm = Tm is an m × m tridiagonal matrix and satisfies

AQm = QmTm+ rmeTm, Qmrm = 0. (7.1.6) (b) Let kxk2 = 1. If Qm ∈ Cn×m with the first column x and QmQm = Im and satisfies

AQm = QmTm+ rmeTm, where Tm is tridiagonal, then

K[x, A, m] = [x, Ax, · · · , Am−1x] = Qm[e1, Tme1, · · · , Tmm−1e1] (7.1.7) is an QR factorization of K[x, A, m].

Proof: (a) Since

AK(x, A, j) ⊂ K(x, A, j + 1), j < m. (7.1.8) From (7.1.5), we have

Span(q1, · · · , qi) = K(x, A, i), i ≤ m. (7.1.9) So we have

qi+1 ⊥ K(x, A, i)(7.1.8)⊃ AK(x, A, i − 1) = A(span(q1, · · · , qi−1)).

This implies

qi+1 Aqj = 0, j = 1, · · · , i − 1, i + 1 ≤ m.

That is

(QmAQm)ij = (Tm)ij = qiAqj = 0 for i > j + 1.

So Tm is upper Hessenberg and then tridiagonal (since Tm is Hermitian).

It remains to show (7.1.6). Since

[x, Ax, · · · , Am−1x] = QmRm and

AK[x, A, m] = K[x, A, m]





0 0

1 . ..

. .. ...

0 1 0



+ AmxeTm,

(3)

we have

AQmRm = QmRm





0 0

1 . ..

. .. ...

0 1 0



+ QmQmAmxeTm+ (I − QmQm)AmxeTm.

Then

AQm = Qm[Rm





0 0

1 . ..

. .. ...

0 1 0



+ QmAmxeTm]R−1m + (I − QmQm)AmxeTmRm−1

= Qm[Rm





0 0

1 . ..

. .. ...

0 1 0



R−1m + γQmAmxeTm] + γ(I − QmQm)Amx

| {z }

rm

eTm

= QmHm+ rmeTm with Qmrm = 0,

where Hm is an upper Hessenberg matrix. But QmAQm = Hm is Hermitian, so Hm = Tm is tridiagonal.

(b) We check (7.1.7):

x = Qme1 coincides the first column. Suppose that i-th columns are equal, i.e.

Ai−1x = QmTmi−1e1 Aix = AQmTmi−1e1

= (QmTm+ rmeTm)Tmi−1e1

= QmTmi e1+ rmeTmTmi−1e1.

But eTmTmi−1e1 = 0 for i < m. Therefore, Aix = QmTie1 the (i + 1)-th columns are equal.

It is clearly that (e1, Tme1, · · · , Tmm−1e1) is an upper triangular matrix.

Theorem 7.1.1 If x = q1 with kq1k2 = 1 satisfies rank(K[x, A, n]) = n

(that is {x, Ax, · · · , An−1x} are linearly independent), then there exists an unitary matrix Q with first column q1 such that QAQ = T is tridiagonal.

Proof: From Theorem 7.1.1(a) m = n, we have Qm = Q unitary and AQ = QT . Uniqueness: Let QAQ = T , ˜QA ˜Q = ˜T and Q1e1 = ˜Qe1

⇒ K[q1, A, n] = QR = ˜Q ˜R

Q = ˜QD, R = D ˜R.

Substitute Q by QD, where D = diag(²1, · · · , ²n) with |²i| = 1. Then (QD)A(QD) = DQAQD = DT D = tridiagonal.

(4)

So Q is unique up to multiplying the columns of Q by a factor ² with |²| = 1.

In the following paragraph we will investigate the Lanczos algorithm for the real case, i.e., A ∈ Rn×n.

How to find an orthogonal matrix Q = (q1, · · · , qn) with QTQ = Insuch that QTAQ = T = tridiagonal and Q is almost uniquely determined. Let

AQ = QT, (7.1.10)

where

Q = [q1, · · · , qn] and T =





α1 β1 0

β1 . .. ...

. .. ... βn−1

0 βn−1 αn



.

It implies that the j-th column of (7.1.10) forms:

Aqj = βj−1qj−1+ αjqj + βjqj+1, (7.1.11) for j = 1, · · · , n with β0 = βn= 0. By multiplying (7.1.11) by qjT we obtain

qjTAqj = αj. (7.1.12)

Define rj = (A − αjI)qj − βj−1qj−1. Then rj = βjqj+1 with

βj = ±krjk2 (7.1.13)

and if βj 6= 0 then

qj+1 = rjj. (7.1.14)

So we can determine the unknown αj, βj, qj in the following order:

Given q1, α1, r1, β1, q2, α2, r2β2, q3, · · · . The above formula define the Lanczos iterations:

j = 0, r0 = q1 , β0 = 1 , q0 = 0 Do while (βj 6= 0)

qj+1 = rjj , j := j + 1 αj = qjTAqj ,

rj = (A − αjI)qj− βj−1qj−1, βj = krjk2.

(7.1.15)

There is no loss of generality in choosing the βj to be positive. The qj are called Lanczos vectors. With careful overwriting and use of the formula αj = qjT(Aqj − βj−1qj−1), the whole process can be implemented with only a pair of n-vectors.

(5)

Algorithm 7.1.1 (Lanczos Algorithm) Given a symmetric A ∈ Rn×n and w ∈ Rn having unit 2-norm. The following algorithm computes a j × j symmetric tridiagonal matrix Tj with the property that σ(Tj) ⊂ σ(A). The diagonal and subdiagonal elements of Tj are stored in α1, · · · , αj and β1, · · · , βj−1 respectively.

vi := 0 (i = 1, · · · , n) β0 := 1

j := 0

Do while (βj 6= 0) if (j 6= 0), then

for i = 1, · · · , n,

t := wi, wi := vij, vi := −βjt.

end for end if

v := Aw + v, j := j + 1, αj := wTv, v := v − αjw, βj := kvk2.

Remark 7.1.3 (a) If the sparity is exploited and only kn flops are involved in each call (Aw) (k ¿ n), then each Lanczos step requires about (4+k)n flops to execute.

(b) The iteration stops before complete tridiagonalizaton if q1 is contained in a proper invariant subspace. From the iteration (7.1.15) we have

A(q1, · · · , qm) = (q1, · · · , qm)





α1 β1

β1 . .. ... βm−1 . .. ...

βm−1 αm



+ (0, · · · , 0,

rm

z }| { βmqm+1)

| {z }

rmeTm

βm = 0 if and only if rm = 0.

This implies

A(q1, · · · , qm) = (q1, · · · , qm)Tm. That is

Range(q1, · · · , qm) = Range(K[q1, A, m])

is the invariant subspace of A and the eigenvalues of Tm are the eigenvalues of A.

Theorem 7.1.2 Let A be symmetric and q1 be a given vector with kq1k2 = 1. The Lanc- zos iterations (7.1.15) runs until j = m where m = rank[q1, Aq1, · · · , An−1q1]. Moreover, for j = 1, · · · , m we have

AQj = QjTj + rjeTj (7.1.16) with

Tj =





α1 β1 β1 . .. ...

. .. ... βj−1 βj−1 αj



 and Qj = [q1, · · · , qj]

(6)

has orthonormal columns satisfying Range(Qj) = K(q1, A, j).

Proof: By induction on j. Suppose the iteration has produced Qj = [q1, · · · , qj] such that Range(Qj) = K(q1, A, j) and QTjQj = Ij. It is easy to see from (7.1.15) that (7.1.16) holds. Thus

QTjAQj = Tj + QTjrjeTj. Since αi = qiTAqi for i = 1, · · · , j and

qTi+1Aqi = qi+1T iqi+1+ αiqi+ βi−1qi−1) = qi+1T iqi+1) = βi

for i = 1, · · · , j − 1 we have QTjAQj = Tj. Consequently QTjrj = 0.

If rj 6= 0 then qj+1 = rj/krjk2 is orthogonal to q1, · · · , qj and qj+1 ∈ Span{Aqj, qj, qj−1} ⊂ K(q1, A, j + 1).

Thus QTj+1Qj+1 = Ij+1 and Range(Qj+1) = K(q1, A, j + 1).

On the other hand, if rj = 0, then AQj = QjTj. This says that Range(Qj) = K(q1, A, j) is invariant. From this we conclude that j = m = dim[K(q1, A, n)].

Encountering a zero βj in the Lanczos iteration is a welcome event in that it signals the computation of an exact invariant subspace. However an exactly zero or even small βj is rarely in practice. Consequently, other explanations for the convergence of Tj0s eigenvalues must be sought.

Theorem 7.1.3 Suppose that j steps of the Lanczos algorithm have been performed and that

SjTTjSj = diag(θ1, · · · , θj)

is the Schur decomposition of the tridiagonal matrix Tj, if Yj ∈ Rn×j is defined by Yj = [y1, · · · , yj] = QjSj

then for i = 1, · · · , j we have

kAyi− θiyik2 = |βj||sji| where Sj = (spq).

Proof: Post-multiplying (7.1.16) by Sj gives

AYj = Yjdiag(θ1, · · · , θj) + rjeTjSj, i.e.,

Ayi = θiyi + rj(eTjSjei) , i = 1, · · · , j.

The proof is complete by taking norms and recalling krjk2 = |βj|.

Remark 7.1.4 The theorem provides error bounds for Tj0s eigenvalues:

µ∈σ(A)min i− µ| ≤ |βj||sji| i = 1, · · · , j.

(7)

Note that in section 10 the (θi, yi) are Ritz pairs for the subspace R(Qj).

If we use the Lanczos method to compute AQj = QjTj + rjeTj and set E = τ wwT where τ = ±1 and w = aqj + brj, then it can be shown that

(A + E)Qj = Qj(Tj + τ a2ejeTj) + (1 + τ ab)rjeTj. If 0 = 1 + τ ab, then the eigenvalues of the tridiagonal matrix

T˜j = Tj + τ a2ejeTj

are also eigenvalues of A + E. We may then conclude from theorem 6.1.2 that the interval i(Tj), λi−1(Tj)] where i = 2, · · · , j, each contains an eigenvalue of A + E.

Suppose we have an approximate eigenvalue ˜λ of A. One possibility is to choose τ a2 so that

det( ˜Tj − ˜λIj) = (αj+ τ a2− ˜λ)pj−1λ) − βj−12 pj−2λ) = 0, where the polynomial pi(x) = det(Ti− xIi) can be evaluated at ˜λ using (5.3).

The following theorems are known as the Kaniel-Paige theory for the estimation of eigenvalues which obtained via the Lanczos algorithm.

Theorem 7.1.4 Let A be n × n symmetric matrix with eigenvalues λ1 ≥ · · · ≥ λn and corresponding orthonormal eigenvectors z1, · · · , zn. If θ1 ≥ · · · ≥ θj are the eigenvalues of Tj obtained after j steps of the Lanczos iteration, then

λ1 ≥ θ1 ≥ λ1 −(λ1− λn) tan (φ1)2 [cj−1(1 + 2ρ1)]2 ,

where cos φ1 = |q1Tz1|, ρ1 = (λ1− λ2)/(λ2− λn) and cj−1 is the Chebychev polynomial of degree j − 1.

Proof: From Courant-Fischer theorem we have θ1 = max

y6=0

yTTjy

yTy = max

y6=0

(Qjy)TA(Qjy)

(Qjy)T(Qjy) = max

06=w∈K(q1,A,j)

wTAw wTw .

Since λ1 is the maximum of wTAw/wTw over all nonzero w, it follows that λ1 ≥ θ1. To obtain the lower bound for θ1, note that

θ1 = max

p∈Pj−1

q1Tp(A)Ap(A)q1 q1Tp(A)2q1 , where Pj−1 is the set of all j − 1 degree polynomials. If

q1 = Xn

i=1

dizi

then q1Tp(A)Ap(A)q1

qT1p(A)2q1 = Pn

i=1d2ip(λi)2λi Pn

i=1d2ip(λi)2

(8)

≥ λ1− (λ1− λn)

Pn

i=2d2ip(λi)2 d21p(λ1)2+Pn

i=2d2ip(λi)2.

We can make the lower bound tight by selecting a polynomial p(x) that is large at x = λ1

in comparison to its value at the remaining eigenvalues. Set

p(x) = cj−1[−1 + 2 x − λn

λ2− λn],

where cj−1(z) is the (j − 1)-th Chebychev polynomial generated by cj(z) = 2zcj−1(z) − cj−2(z), c0 = 1, c1 = z.

These polynomials are bounded by unity on [-1,1]. It follows that |p(λi)| is bounded by unity for i = 2, · · · , n while p(λ1) = cj−1(1 + 2ρ1). Thus,

θ1 ≥ λ1− (λ1− λn)(1 − d21) d21

1

c2j−1(1 + 2ρ1).

The desired lower bound is obtained by noting that tan (φ1)2 = (1 − d21)/d21.

Corollary 7.1.5 Using the same notation as Theorem 7.1.4

λn≤ θj ≤ λn+1− λn) tan2n) c2j−1(1 + 2ρn) , where ρn= (λn−1− λn)/(λ1− λn−1) and cos (φn) = |q1Tzn|.

Proof: Apply Theorem 7.1.4 with A replaced by −A.

Example 7.1.1

Lj−1 1

[Cj−1(2λλ1

2 − 1)]2 1

[Cj−1(1 + 2ρ1)]2 Rj−1 = (λ2

λ1)2(j−1) power method

λ12 j=5 j=25

1.5 1.1 × 10−4/3.9 × 10−2 1.4 × 10−27/3.5 × 10−9 Lj−1/Rj−1

1.01 5.6 × 10−1/9.2 × 10−1 2.8 × 10−4/6.2 × 10−1 Lj−1/Rj−1

Rounding errors greatly affect the behavior of algorithm 7.1.1, the Lanczos iteration.

The basic difficulty is caused by loss of orthogonality among the Lanczos vectors. To avoid these difficulties we can reorthogonalize the Lanczos vectors.

(9)

### 7.1.1 Reorthogonalization

Since

AQj = QjTj+ rjeTj, let

AQj − QjTj = rjeTj + Fj (7.1.17) I − QTjQj = CjT + ∆j+ Cj, (7.1.18) where Cj is strictly upper triangular and ∆jis diagonal. (For simplicity, suppose (Cj)i,i+1 = 0 and ∆i = 0.)

Definition 7.1.1 θi and yi ≡ Qjsi are called Ritz value and Ritz vector, respectively, if Tjsi = θisi.

Let Θj ≡ diag(θ1, · · · , θj) = SjTTjSj where Sj

s1 · · · sj

¤.

Theorem 7.1.6 (Paige Theorem) Assume that (a) Sj and Θj are exact ! (Since j ¿ n). (b) local orthogonality is maintained. ( i.e. qTi+1qi = 0, i = 1, . . . , j − 1, rTjqj = 0, and (Cj)i,i+1 = 0 ). Let

FjTQj− QTjFj = Kj− KjT,

jTj − Tjj ≡ Nj− NjT,

Gj = SjT(Kj + Nj)Sj ≡ (rik).

Then

(a) yTi qj+1 = riiji, where yi = Qjsi, βji = βjsji. (b) For i 6= k,

i− θk)yiTyk = rii(sjk

sji) − rkk(sji

sjk) − (rik− rki). (7.1.19) Proof: Multiplied (7.1.17) from left by QTj, we get

QTjAQj − QTjQjTj = QTjrjeTj + QTjFj, (7.1.20) which implies that

QTjATQj − TjQTjQj = ejrjTQj + FjTQj. (7.1.21) Subtracted (7.1.20) from (7.1.21), we have

(QTjγj)eTj − ej(QTjγj)T

= (CjTTj − TjCjT) + (CjTj − TjCj) + (∆jTj− Tjj) + FjTQj − QjFjT

= (CjTTj − TjCjT) + (CjTj − TjCj) + (Nj − NjT) + (Kj− KjT).

(10)

This implies that

(QTjrj)eTj = CjTj− TjCj+ Nj+ Kj. Thus,

yTi qj+1βji = sTi (QTjrj)eTjsi = sTi (CjTj− TjCj)si+ sTi (Nj + Kj)si

= (sTi Cjsii− θi(sTi Cjsi) + rii, which implies that

yTi qj+1 = rii

βji.

Similarly, (7.1.19) can be obtained by multiplying (7.1.20) from left and right by sTi and si, respectively.

Remark 7.1.5 Since yTi qj+1 = rii

βji =

½ O(esp), if |βji| = O(1), (not converge!)

O(1), if |βji| = O(esp), (converge for (θj, yj))

we have that qTj+1yi = O(1) when the algorithm converges, i.e., qj+1 is not orthogonal to

< Qj > where Qjsi = yi.

(i) Full Reorthogonalization by MGS:

Orthogonalize qj+1 to all q1, · · · , qj by

qj+1 := qj+1 Xj

i=1

(qj+1T qi)qi.

If we incorporate the Householder computations into the Lanczos process, we can produce Lanczos vectors that are orthogonal to working accuracy:

r0 := q1 (given unit vector)

Determine P0 = I − 2v0v0T/v0Tv0 so that P0r0 = e1; α1 := q1TAq1;

Do j = 1, · · · , n − 1,

rj := (A − αj)qj− βj−1qj−1 0q0 ≡ 0), w := (Pj−1· · · P0)rj,

Determine Pj = I − 2vjvTj/vTjvj such that Pjw = (w1, · · · , wj, βj, 0, · · · , 0)T, qj+1 := (P0· · · Pj)ej+1,

αj+1 := qTj+1Aqj+1.

This is the complete reorthogonalization Lanczos scheme.

(11)

(ii) Selective Reorthogonalization by MGS If |βji| = O(√

eps), (θj, yj) “good” Ritz pair Do qj+1⊥q1, . . . , qj

Else not to do Reorthogonalization (iii) Restart after m-steps

(Do full Reorthogonalization) (iv) Partial Reorthogonalization

Do reorthogonalization with previous (e.g. k = 5) Lanczos vectors {q1, . . . , qk} For details see the books:

Parlett: “Symmetric Eigenvalue problem” (1980) pp.257–

Golub & Van Loan: “Matrix computation” (1981) pp.332–

To (7.1.19): The duplicate pairs can occur!

i 6= k, (θi− θk) y|{z}iTyk = O(esp)

O(1), if yi = yk ⇒ Qi ≈ Qk

How to avoid the duplicate pairs ? The answer is using the implicit Restart Lanczos algorithm:

Let

AQj = QjTj + rjeTj be a Lanczos decomposition.

• In principle, we can keep expanding the Lanczos decomposition until the Ritz pairs have converged.

• Unfortunately, it is limited by the amount of memory to storage of Qj.

• Restarted the Lanczos process once j becomes so large that we cannot store Qj. – Implicitly restarting method

• Choose a new starting vector for the underlying Krylov sequence

• A natural choice would be a linear combination of Ritz vectors that we are interested in.

### 7.1.2 Filter polynomials

Assume A has a complete system of eigenpairs (λi, xi) and we are interested in the first k of these eigenpairs. Expand u1 in the form

u1 = Xk

i=1

γixi+ Xn i=k+1

γixi. If p is any polynomial, we have

p(A)u1 = Xk

i=1

γip(λi)xi+ Xn i=k+1

γip(λi)xi.

(12)

• Choose p so that the values p(λi) (i = k +1, . . . , n) are small compared to the values p(λi) (i = 1, . . . , k).

• Then p(A)u1 is rich in the components of the xi that we want and deficient in the ones that we do not want.

• p is called a filter polynomial.

• Suppose we have Ritz values µ1, . . . , µm and µk+1, . . . , µm are not interesting. Then take

p(t) = (t − µk+1) · · · (t − µm).

### 7.1.3 Implicitly restarted algorithm

Let

AQm = QmTm+ βmqm+1eTm (7.1.22) be a Lanczos decomposition with order m. Choose a filter polynomial p of degree m − k and use the implicit restarting process to reduce the decomposition to a decomposition

A ˜Qk = ˜QkT˜k+ ˜βkq˜k+1eTk of order k with starting vector p(A)u1.

Let ν1, . . . , νm be eigenvalues of Tm and suppose that ν1, . . . , νm−k correspond to the part of the spectrum we are not interested in. Then take

p(t) = (t − ν1)(t − ν2) · · · (t − νm−k).

The starting vector p(A)u1 is equal to

p(A)u1 = (A − νm−kI) · · · (A − ν2I)(A − ν1I)u1

= (A − νm−kI) [· · · [(A − ν2I) [(A − ν1I)u1]]] .

In the first, we construct a Lanczos decomposition with starting vector (A−ν1I)u1. From (7.1.22), we have

(A − ν1I)Qm = Qm(Tm− ν1I) + βmqm+1eTm (7.1.23)

= QmU1R1+ βmqm+1eTm, where

Tm− ν1I = U1R1

is the QR factorization of Tm− κ1I. Postmultiplying by U1, we get (A − ν1I)(QmU1) = (QmU1)(R1U1) + βmqm+1(eTmU1).

It implies that

AQ(1)m = Q(1)m Tm(1)+ βmqm+1b(1)Tm+1, where

Q(1)m = QmU1, Tm(1) = R1U1+ ν1I, b(1)Tm+1 = eTmU1. (Q(1)m : one step of single shifted QR algorithm)

(13)

Remark 7.1.6

• Q(1)m is orthonormal.

• By the definition of Tm(1), we get

U1Tm(1)U1T = U1(R1U1+ ν1I)U1T = U1R1 + ν1I = Tm. (7.1.24) Therefore, ν1, ν2, . . . , νm are also eigenvalues of Tm(1).

• Since Tm is tridiagonal and U1 is the Q-factor of the QR factorization of Tm− ν1I, it implies that U1 and Tm(1) are upper Hessenberg. From (7.1.24), Tm(1) is symmetric.

Therefore, Tm(1) is also tridiagonal.

• The vector b(1)Tm+1 = eTmU1 has the form b(1)Tm+1 =

h

0 · · · 0 Um−1,m(1) Um,m(1)

i

; i.e., only the last two components of b(1)m+1 are nonzero.

• For on postmultiplying (7.1.23) by e1, we get

(A − ν1I)q1 = (A − ν1I)(Qme1) = Q(1)m R1e1 = r(1)11q1(1).

Since Tm is unreduced, r(1)11 is nonzero. Therefore, the first column of Q(1)m is a multiple of (A − κ1I)q1.

Repeating this process with ν2, . . . , νm−k, the result will be a Krylov decomposition AQ(m−k)m = Q(m−k)m Tm(m−k)+ βmqm+1b(m−k)Tm+1

with the following properties i. Q(m−k)m is orthonormal.

ii. Tm(m−k) is tridiagonal.

iii. The first k − 1 components of b(m−k)Tm+1 are zero.

iv. The first column of Q(m−k)m is a multiple of (A − ν1I) · · · (A − νm−kI)q1.

Corollary 7.1.1 Let ν1, . . . , νm be eigenvalues of Tm. If the implicitly restarted QR step is performed with shifts ν1, . . . , νm−k, then the matrix Tm(m−k) has the form

Tm(m−k) =

"

Tkk(m−k) Tk,m−k(m−k) 0 Tk+1,k+1(m−k)

# ,

where Tk+1,k+1(m−k) is an upper triangular matrix with Ritz value ν1, . . . , νm−k on its diagonal.

(14)

Therefore, the first k columns of the decomposition can be written in the form AQ(m−k)k = Q(m−k)k Tkk(m−k)+ tk+1,kqk+1(m−k)eTk + βkumkqm+1eTk,

where Q(m−k)k consists of the first k columns of Q(m−k)m , Tkk(m−k) is the leading principal submatrix of order k of Tm(m−k), and ukm is from the matrix U = U1· · · Um−k. Hence if we set

Q˜k = Q(m−k)k , T˜k = Tkk(m−k),

β˜k = ktk+1,kqk+1(m−k)+ βkumkqm+1k2,

˜

qk+1 = ˜βk−1(tk+1,kqk+1(m−k)+ βkumkqm+1), then

A ˜Qk = ˜QkT˜k+ ˜βkq˜k+1eTk

is a Lanczos decomposition whose starting vector is proportional to (A − ν1I) · · · (A − νm−kI)q1.

• Avoid any matrix-vector multiplications in forming the new starting vector.

• Get its Lanczos decomposition of order k for free.

• For large n the major cost will be in computing QU .

### 7.2 Approximation from a subspace

Assume that A is symmetric and {(αi, zi)}ni=1be eigenpairs of A with α1 ≤ α2 ≤ · · · ≤ αn. Define

ρ(x) = ρ(x, A) = xTAx xTx . Algorithm 7.2.1 (Rayleigh-Ritz-Quotient procedure)

Give a subspace S(m) = span{Q} with QTQ = Im; Set H := ρ(Q) = QTAQ;

Compute the p (≤ m) eigenpairs of H, which are of interest, say Hgi = θigi for i = 1, . . . , p;

Compute Ritz vectors yi = Qgi, for i = 1, . . . , p;

Check kAyi− θiyik2 ≤ T ol, for i = 1, . . . , p.

By the minimax characterization of eigenvalues, we have αj = λj(A) = min

Fj⊆Rnmax

f ∈Fjρ(f, A).

(15)

Define

βj = min

Gj⊆Smmax

g∈Gjρ(g, A), for j ≤ m.

Since Gj ⊆ Sm and S(m) = span{Q}, it implies that Gj = Q eGj for some eGj ⊆ Rm. Therefore,

βj = min

Gej⊆Rmmax

s∈ eGjρ(s, H) = λj(H) ≡ θj, for j = 1, . . . , m.

For any m by m matrix B, there is associated a residual matrix R(B) ≡ AQ − QB.

Theorem 7.2.1 For given orthonormal n by m matrix Q, kR(H)k ≤ kR(B)k

for all m by m matrix B.

Proof: Since

R(B)R(B) = QA2Q − B(QAQ) − (QAQ)B + BB

= QA2Q − H2+ (H − B)(H − B)

= R(H)R(H) + (H − B)(H − B)

and (H − B)(H − B) is positive semidefinite, it implies that kR(B)k2 ≥ kR(H)k2. Since

Hgi = θigi, for i = 1, . . . , m, we have that

QTAQgi = θigi, which implies that

QQTA(Qgi) = θi(Qgi).

Let yi = Qgi. Then QQTyi = Q(QTQ)gi = yi. Take PQ = QQT which is a projection on span{Q}. Then

(QQT)Ayi = θi(QQT)yi, which implies that

PQ(Ayi− θiyi) = 0, i.e., ri = Ayi − θiyi ⊥ Sm = span{Q}.

– The The readLine readLine method is the same method used to read method is the same method used to read  from the keyboard, but in this case it would read from a

• The start node representing the initial conﬁguration has zero in degree.... The Reachability

In this paper, we propose a practical numerical method based on the LSM and the truncated SVD to reconstruct the support of the inhomogeneity in the acoustic equation with

The hashCode method for a given class can be used to test for object equality and object inequality for that class. The hashCode method is used by the java.util.SortedSet

For the proposed algorithm, we establish its convergence properties, and also present a dual application to the SCLP, leading to an exponential multiplier method which is shown

/** Class invariant: A Person always has a date of birth, and if the Person has a date of death, then the date of death is equal to or later than the date of birth. To be

Biases in Pricing Continuously Monitored Options with Monte Carlo (continued).. • If all of the sampled prices are below the barrier, this sample path pays max(S(t n ) −

To convert a string containing floating-point digits to its floating-point value, use the static parseDouble method of the Double class..