• 沒有找到結果。

# Polynomial Jacobi Davidson Method for Large/Sparse Eigenvalue Problems

N/A
N/A
Protected

Share "Polynomial Jacobi Davidson Method for Large/Sparse Eigenvalue Problems"

Copied!
43
0
0

(1)

### Polynomial Jacobi Davidson Method for Large/Sparse Eigenvalue Problems

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University, Taiwan

April 28, 2011

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 1 / 42

(2)

### Outline

1 Jacobi’s orthogonal component correction

2 Davidson’s method

3 Jacobi-Davidson method

4 Polynomial Jacobi-Davidson method

5 Non-equivalence deflation

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 2 / 42

(3)

References

SIAM J. Matrix Anal. Vol. 17, No. 2, PP. 401-425, 1996, Van der Vorst etc.

BIT. Vol. 36, No. 3, PP. 595-633, 1996, Van der Vorst etc.

T.-M. Hwang, W.-W. Lin and V. Mehrmann, SIAM J. Sci. Comput.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 3 / 42

(4)

### Some basic theorems

Theorem 1

Let X be an eigenspace of A and let X be a basis for X . Then there is a unique matrix L such that

AX = XL.

The matrix L is given by

L = XIAX, where XI is a matrix satisfying XIX = I.

If (λ, x) is an eigenpair of A with x ∈ X , then (λ, XIx)is an eigenpair of L. Conversely, if (λ, u) is an eigenpair of L, then (λ, Xu) is an eigenpair of A.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 4 / 42

(5)

Proof: Let

X = [x1· · · xk] and Y = AX = [y1· · · yk] .

Since yi∈ X and X is a basis for X , there is a unique vector `i such that

yi = X`i. If we set L = [`1· · · `k], then AX = XL and

L = XIXL = XIAX.

Now let (λ, x) be an eigenpair of A with x ∈ X . Then there is a unique vector u such that x = Xu. However, u = XIx. Hence

λx = Ax = AXu = XLu ⇒ λu = λXIx = Lu.

Conversely, if Lu = λu, then

A(Xu) = (AX)u = (XL)u = X(Lu) = λ(Xu), so that (λ, Xu) is an eigenpair of A.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 5 / 42

(6)

Theorem 2 (Optimal residuals) Let [X X]be unitary. Let

R = AX − XL and SH = XHA − LXH. Then kRk and kSk are minimized when

L = XHAX, in which case

(a) kRk = kXHAXk, (b) kSk = kXHAXk, (c) XHR = 0.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 6 / 42

(7)

Proof: Set

 XH XH

 A

X X  =

 Lˆ H

G M

 . Then

 XH XH

 R =

 Lˆ H

G M

  XH XH

 X −

 XH XH

 XL =

 L − Lˆ G

 . It implies that

kRk =

 XH XH

 R

=

 L − Lˆ G

 , which is minimized when L = ˆL = XHAXand

min kRk = kGk = kXHAXk.

The proof for S is similar. If L = XHAX, then

XHR = XHAX − XHXL = XHAX − L = 0.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 7 / 42

(8)

Definition 3

Let X be of full column rank and let XI be a left inverse of X. Then XIAX is a Rayleigh quotient of A.

Theorem 4

Let X be orthonormal, A be Hermitian and R = AX − XL.

If `1, . . . , `kare the eigenvalues of L, then there are eigenvalues λj1, . . . , λjk of A such that

|`i− λji| ≤ kRk2 and v u u t

k

X

i=1

(`i− λji)2≤√

2kRkF.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 8 / 42

(9)

### Jacobi’s orthogonal component correction, 1846

Consider the eigenvalue problem A

 1 z



 α cT b F

  1 z



= λ

 1 z



, (1)

where A is diagonal dominant and α is the largest diagonal element.

(1) is equivalent to

 λ = α + cTz, (F − λI)z = −b.

Jacobi iteration: (with z1 = 0)

 θk= α + cTzk,

(D − θkI)zk+1 = (D − F )zk− b (2) where D = diag(F ).

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 9 / 42

(10)

### Davidson’s method (1975)

Algorithm 1 (Davidson’s method) Given unit vector v, set V = [v]

Iterate until convergence

Compute desired eigenpair (θ, s) ofVTAV. Computeu = V sandr = Au − θu.

If (k r k2< ε), stop.

Solve(DA− θI)t = r.

Orthog. t ⊥ V → v, V = [V, v]

end

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 10 / 42

(11)

Let uk= (1, zkT)T. Then

rk = (A − θkI)uk=

 α − θk+ cTzk (F − θkI)zk+ b



Substituting the residual vector rkinto linear systems (DA− θkI)tk= −rk, where DA=

 α 0

0 D

 , we get

(D − θkI)yk = −(F − θkI)zk− b

= (D − F )zk− (D − θkI)zk− b From (2) and above equality, we see that

(D − θkI)(zk+ yk) = (D − F )zk− b = (D − θkI)zk+1.

This implies that zk+1= zk+ ykas one step of JOCC starting with zk.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 11 / 42

(12)

### Jacobi-Davidson method (1996)

k, uk): approx. eigenpair of A, θk≈ λ, with

uk= Vksk, VkTAVksk= θksk and k skk2= 1.

Definition 5

k, uk)is called a Ritz pair of A. θkis called a Ritz value and ukis a Ritz vector.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 12 / 42

(13)

### Jacobi-Davidson method (1996)

k, uk): approx. eigenpair of A, θk≈ λ, with

uk= Vksk, VkTAVksk= θksk and k skk2= 1.

Then

uTkrk ≡ uTk(A − θkI)uk = sTkVkTAVksk− θksTkVkTVksk = 0 ⇒ rk⊥ uk Find the correctiont ⊥ uksuch that

A(uk+ t) = λ(uk+ t).

That is

(A − λI)t = λuk− Auk= −rk+ (λ − θk)uk.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 12 / 42

(14)

Since (I − ukuTk)rk= rk, we have

I − ukuTk (A − λI)t = −rk.

Correction equation

(I − ukuTk)(A − θkI) I − ukuTk t = −rk and t ⊥ uk,

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 13 / 42

(15)

### Solving correction vector t

Correction equation:

I − ukuTk (A − θkI)(I − ukuTk)t = −rk, t ⊥ uk. (3) Scheme SOneLS:

Use preconditioning iterative approximations, e.g., GMRES, to solve (3).

Use a preconditioner

Mp ≡ I − ukuTk M I − ukuTk , where M is an approximation of A − θkI. In each of the iterative steps, it needs to solve

Mpy = b, y ⊥ uk (4)

for a given b.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 14 / 42

(16)

Since y ⊥ uk, Eq. (4) can be rewritten as

I − ukuTk My = b ⇒ My = uTkMy uk+ b ≡ ηkuk+ b.

Hence

y = M−1b + ηkM−1uk, where

ηk= − uTkM−1b uTkM−1uk

.

SSOR preconditioner: Let A − θkI = L + D + U. Then M = (D + ωL)D−1(D + ωU ).

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 15 / 42

(17)

Scheme ST woLS: Since t ⊥ uk, Eq. (3) can be rewritten as

(A − θkI)t = uTk(A − θkI)t uk− rk≡ εuk− rk. (5) Let t1and t2 be approximated solutions of the following linear systems:

(A − θkI)t = −rk and (A − θkI)t = uk, respectively. Then the approximated solution ˜tfor (5) is

˜t = t1+ εt2 for ε = −uTkt1 uTkt2

. Scheme SOneStep: The approximated solution ˜tfor (5) is

˜t = −M−1rk+ εM−1uk for ε = u>kM−1rk u>kM−1uk, where M is an approximation of A − θkI.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 16 / 42

(18)

Algorithm 2 (Jacobi-Davidson Method) Choose an n-by-m orthonormal matrix V0

Do k = 0, 1, 2, · · ·

Compute all the eigenpairs ofVkTAVks = λs.

Select the desired (target) eigenpair (θk, sk)with kskk2 = 1.

Computeuk = Vkskandrk = (A − θkI)uk. If (krkk2< ε), λ = θk, x = uk, Stop

Solve (approximately) atk⊥ uk from

(I − ukuTk)(A − θkI)(I − ukuTk)t = −rk. Orthogonalize tk⊥ Vk→ vk+1, Vk+1= [Vk, vk+1]

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 17 / 42

(19)

Locking:

Vkwith VkV = Ikare convergentSchur vectors, i.e., AVk = VkTk

for someupper triangularTk. SetV = [Vk, Vq]with VV = Ik+q in k + 1-th iteration of Jacobi-Davidson Algorithm. Then

VAV =

 VkAVk VkAVq VqAVk VqAVq



=

 Tk VkAVq VqVkTk VqAVq



=

 Tk VkAVq 0 VqAVq

 . Restarting:

Keep the locked Schur vectors as well as the Schur vectors of interest in the subspace and throw away those we are not interested.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 18 / 42

(20)

Remark 1

If ε = 0, we obtain Davidson method with t1 = −(DA− θkI)−1r.

( ˜tis NOT orthogonal to uk )

If the linear systems in (6) are exactly solved, then the vector t becomes

t = ε(A − θkI)−1uk− uk. (6) Since t is made orthogonal to uk, (6) is equivalent to

t = (A − θkI)−1uk

which is equivalent to shift-invert power iteration. Hence it is quadratic convergence.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 19 / 42

(21)

Consider Ax = λx and assume that λ is simple.

Lemma 6

Consider w with wTx 6= 0. Then the map Fp



I −xwT wTx



(A − λI)



I −xwT wTx



is a bijection from w to w.

Proof: Suppose y⊥w and Fpy = 0. That is



I − xwT wTx



(A − λI)



I −xwT wTx

 y = 0.

Then it holds that

(A − λI)y = εx.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 20 / 42

(22)

Therefore, both y and x belong to the kernel of (A − λI)2. The

simplicity of λ implies that y is a scale multiple of x. The fact that y⊥w and xTw 6= 0implies y = 0, which proves the injectivity of Fp. An obvious dimension argument implies bijectivity.

Extension Fpt ≡



I − uuT uTu



(A − θI)



I − uuT uTu



t = −r, t ⊥ u, r ⊥ u.

Then

t ∈ u

Fp

r ∈ u.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 21 / 42

(23)

Theorem 7

Assume that the correction equation

I − uuT (A − θI) I − uuT t = −r, t ⊥ u (7) is solved exactly in each step of Jacobi-Davidson Algorithm. Assume uk= u → xand uTkxhas non-trivial limit. Then if ukis sufficiently chosen to x, then uk→ x with locally quadratical convergence and

θk= uTkAuk → λ.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 22 / 42

(24)

Proof: Suppose Ax = λx with x such that x = u + z for z ⊥ u. Then (A − θI) z = − (A − θI) u + (λ − θ) x = −r + (λ − θ) x. (8) Consider the exact solution z1 ⊥ u of (7):

(I − P ) (A − θI) z1 = − (I − P ) r, (9) where P = uuT. Note that (I − P ) r = r since u⊥r. Since

x − (u + z1) = z − z1 and z = x − u, for quadratic convergence, it suffices to show that

kx − (u + z1) k = kz − z1k = O kzk2 . (10) Multiplying (8) by (I − P ) and subtracting the result from (9) yields

(I − P ) (A − θI) (z − z1) = (λ − θ) (I − P ) z + (λ − θ) (I − P ) u. (11)

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 23 / 42

(25)

Multiplying (8) by uT and using r ⊥ u leads to λ − θ = uT(A − θI) z

uTx . (12)

Since uTkxhas non-trivial limit, we obtain k (λ − θ) (I − P ) zk =

uT (A − θI) z

uTx (I − P ) z

. (13)

From (11), Lemma 6 and (I − P ) u = 0, we have kz − z1k =

h

(I − P ) (A − θkI) |u

k

i−1

(λ − θ) (I − P ) z

= h

(I − P ) (A − θkI) |u

k

i−1 uTk (A − θkI) z

uTkx (I − P ) z

= O kzk2 .

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 24 / 42

(26)

### Jacobi Davidson method as on accelerated Newton Scheme

Consider Ax = λx and assume that λ is simple. Choose wTx = 1.

Consider nonlinear equation

F (u) = Au − θ(u)u = 0 with θ (u) = wTAu wTu ,

where kuk = 1 or wTu = 1. Then F : {u|wTu = 1} → w. In particular, r ≡ F (u) = Au − θ (u) u ⊥ w.

Suppose uk≈ x and the next Newton approximation uk+1: uk+1 = uk− ∂F

∂u u=uk

!−1

F (uk) is given by uk+1= uk+ t, i.e., t satisfies that

∂F

∂u u=uk

!

t = F (uk) = −r.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 25 / 42

(27)

Since 1 = uTk+1w = (uk+ t)Tw = 1 + tTw, it implies that wTt = 0. By the definition of F , we have

∂F

∂u = A − θ (u) I −− wTAu uwT + wTu uwTA (wTu)2

= A − θI + wTAu

(wTu)2uwT −uwTA wTu =



I −ukwT wTuk



(A − θkI) . Consequently, the Jacobian of F acts on wand is given by

∂F

∂u u=uk

! t =



I − ukwT wTuk



(A − θkI) t, t ⊥ w.

Hence the correction equation of Newton method read as t ⊥ w,



I − ukwT wTuk



(A − θkI) t = −r,

which is the correction equation of Jacobi-Davidson method in (7) with w = u.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 26 / 42

(28)

### Polynomial Jacobi-Davidson method

k, uk): approx. eigenpair of A(λ) ≡Pτ

k=0λkAk, θk≈ λ, with uk= Vksk, Vk>A(λ)Vksk= 0 and k skk2= 1.

Let

rk= A(θk)uk. Then

u>krk= u>kA(θk)uk= s>kVk>A(θk)Vksk= 0 ⇒ rk ⊥ uk Find the correctiontsuch that

A(λ)(uk+ t) = 0.

That is

A(λ)t = −A(λ)uk = −rk+ (A(θk) − A(λ))uk.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 27 / 42

(29)

Let

pk = A0k)uk

τ

X

i=1

i−1k Ai

! uk.

A(λ) = A − λI: pk= −uk,

(A(θk) − A(λ))uk = (λ − θk)uk= (θk− λk)pk

A(λ) = A − λB:

pk= −Buk,

(A(θk) − A(λ))uk= (λ − θk)Buk= (θk− λ)pk A(λ) =Pτ

i=0λiAiwith τ ≥ 2:

(A(θk) − A(λ))uk =



k− λ)A0k) −1

2(θk− λ)2A00k)

 uk

= (θk− λ)pk−1

2(θk− λ)2A00k)uk

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 28 / 42

(30)

Hence

A(λ)t = −rk+ (θk− λ)pk−1

2(θk− λ)2A00k)uk Since rk⊥ uk, we have



I − pku>k u>kpk



A(λ)t = −rk−1

2(θk− λ)2



I − pku>k u>kpk



A00k)uk. Correction equation:



I − pku>k u>kpk



A(θk)(I − uku>k)t = −rk and t ⊥ uk,

or



I −pku>k u>kpk



(A − θkB)



I −ukp>k p>kuk



t = −rk and t ⊥Buk, with symmetric positive definite matrix B.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 29 / 42

(31)

Algorithm 3 (Jacobi-Davidson Algorithm for solving A(λ)x = 0) Choose an n-by-m orthonormal matrix V0

Do k = 0, 1, 2, · · ·

Compute all the eigenpairs ofVk>A(λ)Vk = 0.

Select the desired (target) eigenpair (θk, sk)with kskk2 = 1.

Computeuk = Vksk,rk= A(θk)ukandpk= A0k)uk. If (krkk2< ε), λ = θk, x = uk, Stop

Solve (approximately) atk⊥ uk from (I −pukTuTk

kpk)A(θk)(I − ukuTk)t = −rk. Orthogonalize tk⊥ Vk→ vk+1, Vk+1= [Vk, vk+1]

pk = A0k)uk

τ

X

i=1

i−1k Ai

! uk.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 30 / 42

(32)

### Non-equivalence deflation of quadratic eigenproblems

Let λ1 be a real eigenvalue of Q(λ) ≡ λ2M + λC + K and x1, z1be the associated right and left eigenvectors, respectively, with zT1Kx1= 1.

Let

θ1 = (z1TM x1)−1. We introduce a deflated quadratic eigenproblem

Q(λ)x ≡e h

λ2M + λ ef C + eK i

x = 0, where

Mf = M − θ1M x1z1TM, Ce = C + θ1

λ1(M x1z1TK + Kx1z1TM ), Ke = K − θ1

λ21Kx1zT1K.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 31 / 42

(33)

### Complex deflation

Let λ1 = α1+ iβ1be a complex eigenvalue of Q(λ) and

x1= x1R+ ix1I, z1 = z1R+ iz1I be the associated right and left eigenvectors, respectively, such that

Z1TKX1= I2, where X1 = [x1R, x1I]and Z1= [z1R, z1I]. Let

Θ1= (Z1TM X1)−1.

Then we introduce a deflated quadratic eigenproblem with Mf = M − M X1Θ1Z1TM,

Ce = C + M X1Θ1Λ−T1 Z1TK + KX1Λ−11 ΘT1Z1TM, Ke = K − KX1Λ−11 Θ1Λ−T1 Z1TK

in which Λ1 =

 α1 β1

−β1 α1

 .

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 32 / 42

(34)

Theorem 8

(i) Let λ1 be a simple real eigenvalue of Q(λ). Then the spectrum of Q(λ)e is given by

(σ(Q(λ)){λ1}) ∪ {∞}

provided that λ216= θ1.

(ii) Let λ1 be a simple complex eigenvalue of Q(λ). Then the spectrum of eQ(λ)is given by

σ(Q(λ)){λ1, ¯λ1} ∪ {∞, ∞}

provided that Λ1ΛT1 6= Θ1.

Furthermore, in both cases (i) and (ii), if λ26= λ1 and (λ2, x2)is an eigenpair of Q(λ) then the pair (λ2, x2)is also an eigenpair of eQ(λ).

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 33 / 42

(35)

Suppose that M, C, K are symmetric. Given an eigenmatrix pair (Λ1, X1) ∈ Rr×r× Rn×r of Q(λ), where Λ1 is nonsingular and X1 satisfies

X1TKX1= Ir, Θ1 := (X1TM X1)−1. We define ˜Q(λ) := λ2M + λ ˜˜ C + ˜K, where

M˜ := M − M X1Θ1X1TM,

C˜ := C + M X1Θ1Λ−T1 X1TK + KX1Λ−11 Θ1X1TM, K˜ := K − KX1Λ−11 Θ1Λ−T1 X1TK.

Theorem 9

Suppose that Θ1− Λ1ΛT1 is nonsingular. Then the eigenvalues of the real symmetric quadratic pencil ˜Q(λ)are the same as those of Q(λ) except that the eigenvalues of Λ1, which are closed under complex conjugation, are replaced by r infinities.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 34 / 42

(36)

Proof: Since (Λ1, X1)is an eigenmatrix pair of Q(λ), i.e., M X1Λ21+ CX1Λ1+ KX1 = 0,

we have

Q(λ)˜ = Q(λ) + [M X1(λIr+ Λ1) + CX1] Θ1Λ−T1 (X1TK − λΛT1X1TM )

= Q(λ) + Q(λ)X1(λIr− Λ1)−1Θ1Λ−T1 (X1TK − λΛT1X1TM ).

By using the identity

det(In+ RS) = det(Im+ SR), where R, ST ∈ Rn×m, we have

det[ ˜Q(λ)]

= det[Q(λ)]det[I + X1(λIr− Λ1)−1Θ1Λ−T1 (X1TK − λΛT1X1TM )]

= det[Q(λ)]det[Ir+ (λIr− Λ1)−1Θ1Λ−T1 (Ir− λΛT1Θ−11 )]

= det[Q(λ)]

det(λIr− Λ1)det(Θ1Λ−T1 − Λ1).

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 35 / 42

(37)

Since (Θ1− Λ1ΛT1) ∈ Rr×ris nonsingular, we have det(Θ1Λ−T1 − Λ1) 6= 0.

Therefore, ˜Q(λ)has the same eigenvalues as Q(λ) except that r eigenvalues of Λ1 are replaced by r infinities.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 36 / 42

(38)

### Non-equiv. deflation for cubic poly. eigenproblems

Let (Λ, Vu) ∈ Rr×r× Rn×r be an eigenmatrix pair of

A(λ) ≡ λ3A3+ λ2A2+ λA1+ A0 (14) with VuTVu= Irand 0 /∈ σ(Λ), i.e.,

A3VuΛ3+ A2VuΛ2+ A1VuΛ + A0Vu = 0. (15) Define a new deflated cubic eigenvalue problem by

A(λ)u = (λe 33+ λ22+ λ ˜A1+ ˜A0)u = 0, (16) where





0 = A0,

1 = A1− (A1VuVuT + A2VuΛVuT + A3VuΛ2VuT), A˜2 = A2− (A2VuVuT + A3VuΛVuT),

3 = A3− A3VuVuT.

(17)

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 37 / 42

(39)

Lemma 10

Let A(λ) and eA(λ)be cubic pencils given by (14) and (16), respectively. Then it holds

A(λ) = A(λ) Ie n− λVu(λIr− Λ)−1VuT . (18) Theorem 11

Let (Λ, Vu)be an eigenmatrix pair of A(λ) with VuTVu = Ir. Then (i) (σ(A(λ))σ(Λ)) ∪ {∞} = σ( eA(λ)).

(ii) Let (µ, z) be an eigenpair of A(λ) with k z k2= 1and µ /∈ σ(Λ).

Define

˜

z = (In− µVuΛ−1VuT)z ≡ T (µ)z. (19) Then (µ, ˜z)is an eigenpair of eA(λ).

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 38 / 42

(40)

Proof of Lemma: Using (17) and (15), and the fundamental matrix calculation, we have

A(λ)e = A(λ) − λ λ2A3VFVFT + λA2VFVFT + λA3VFΛVFT +A1VFVFT + A2VFΛVFT + A3VFΛ2VFT

= A(λ) − λ A3VF(λIr− Λ)3(λIr− Λ)−1VFT +3A3VFΛ(λIr− Λ)2(λIr− Λ)−1VFT +3A3VFΛ2(λIr− Λ)(λIr− Λ)−1VFT +A2VF(λIr− Λ)2(λIr− Λ)−1VFT +2A2VFΛ(λIr− Λ)(λIr− Λ)−1VFT

+A1VF(λIr− Λ)(λIr− Λ)−1VFT

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 39 / 42

(41)

A(λ)e = A(λ) − λA3VF3Ir− Λ3) + A2VF2Ir− Λ2) +A1VF(λIr− Λ) + A0VF − A0VF] (λIr− Λ)−1VFT

= A(λ) − λA(λ)VF(λIr− Λ)−1VFT

= A(λ)In− λVF(λIr− Λ)−1VFT .

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 40 / 42

(42)

Proof of Theorem: (i) Using the identity

det(In+ RS) = det(Im+ SR) and Lemma 10, we have

det( eA(λ)) = det(A(λ)) det In− λVF(λIr− Λ)−1VFT

= det(A(λ)) det In− λ(λIr− Λ)−1

= det(A(λ)) det(λIr− Λ)−1det(−Λ).

Since 0 /∈ σ(Λ), det(−Λ) 6= 0. Thus, eA(λ)and A(λ) have the same finite spectrum except the eigenvalues in σ(Λ). Furthermore, dividing Eq. (16) by λ3 and using the fact that

Ae3VF = (A3− A3VFVFT)VF = 0,

we see that (diagr{∞, · · · , ∞}, VF)is an eigenmatrix pair of eA(λ) corresponding to infinite eigenvalues.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 41 / 42

(43)

(ii) Since µ /∈ σ(Λ), the matrix T (µ) = (I − µVFΛ−1VFT)in (19) is invertible with the inverse

T (µ)−1= In− µVF(µIr− Λ)−1VFT. (20) From Lemma 10, we have

A(µ)˜e z = A(µ)In− µVF(µIr− Λ)−1VFT In− µVFΛ−1VFT z = 0.

This completes the proof.

T.M. Huang (Taiwan Normal Univ.) Poly. JD method for PEP April 28, 2011 42 / 42

The main tool in our reconstruction method is the complex geometri- cal optics (CGO) solutions with polynomial-type phase functions for the Helmholtz equation.. This type of

Tsung-Min Hwang, Wei-Cheng Wang and Weichung Wang, Numerical schemes for three dimensional irregular shape quantum dots over curvilinear coordinate systems, accepted for publication

Arbenz et al.[1] proposed a hybrid preconditioner combining a hierarchical basis preconditioner and an algebraic multigrid preconditioner for the correc- tion equation in the

Keywords: Finite volume method; Heterostructure; Large scale polynomial eigenvalue problem; Semiconductor pyramid quantum dot;.. Schr€

In this chapter we develop the Lanczos method, a technique that is applicable to large sparse, symmetric eigenproblems.. The method involves tridiagonalizing the given

In section 4, based on the cases of circular cone eigenvalue optimization problems, we study the corresponding properties of the solutions for p-order cone eigenvalue

Theorem 5.6.1 The qd-algorithm converges for irreducible, symmetric positive definite tridiagonal matrices.. It is necessary to show that q i are in

We have also discussed the quadratic Jacobi–Davidson method combined with a nonequivalence deﬂation technique for slightly damped gyroscopic systems based on a computation of