師大

### Jacobi Davidson method and its applications

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University, Taiwan

October 31, 2008

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 1 / 71

師大

### Outline

1 Some basic theorems

2 Jacobi-Davidson method for linear eigenvalue problems

3 Polynomial eigenvalue problems

4 Locking method for computing λ2, . . . , λ`

5 Non-equivalence deflation of quadratic eigenproblems

6 Non-equivalence deflation for cubic polynomial eigenproblems

7 Applications

8 Numerical experiments for linear eigenproblems

9 Numerical experiments for polynomial eigenproblems

師大

### Some basic theorems

Eigenproblems

Standard eigenproblems: Ax = λx Generalized eigenproblems: Ax = λBx

Higher order poly. eigenproblems: (A0+ λA1+ .... + λ^{n}An)x = 0
Eigenproblems of λ-matrices: F (λ)x = 0

What do we care ?

(i) In theory: eigenstructure, spectrum decomposition, canonical form, . . . , etc.

(ii) In computation: eigenvalues, eigenvectors, invariant subspaces, . . . , etc.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 3 / 71

師大

Theorem (Fischer)

Let the Hermitian matrix A have eigenvalues
λ_{1} ≥ λ_{2} ≥ · · · ≥ λ_{n}.
Then

λi = min dim(W)=n−i +1

w ∈W,kw kmax2=1w^{H}Aw

and

λ_{i} = max
dim(W)=i

min

w ∈W,kw k2=1w^{H}Aw

.

師大

Theorem

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 = X^{I}AX ,
where X^{I} is a matrix satisfying X^{I}X = I .

If (λ, x ) is an eigenpair of A with x ∈ X , then (λ, X^{I}x ) 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.) JD method and applications October 31, 2008 5 / 71

師大

### Proof:

Let

X = [x_{1}· · · x_{k}] and Y = AX = [y_{1}· · · y_{k}] .

Since y_{i} ∈ X and X is a basis for X , there is a unique vector `_{i} such that
y_{i} = X `_{i}.

If we set L = [`_{1}· · · `_{k}], then AX = XL and
L = X^{I}XL = X^{I}AX .

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

λx = Ax = AXu = XLu ⇒ λu = λX^{I}x = Lu.

Conversely, if Lu = λu, then

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

師大

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

R = AX − XL and S^{H} = X^{H}A − LX^{H}.
Then kRk and kS k are minimized when

L = X^{H}AX ,
in which case

(a) kRk = kX_{⊥}^{H}AX k,
(b) kSk = kX^{H}AX⊥k,
(c) X^{H}R = 0.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 7 / 71

師大

### Proof:

Set

X^{H}
X_{⊥}^{H}

A

X X⊥ =

ˆL H

G M

. Then

X^{H}
X_{⊥}^{H}

R =

ˆL H

G M

X^{H}
X_{⊥}^{H}

X −

X^{H}
X_{⊥}^{H}

XL =

L − Lˆ G

. It implies that

kRk =

X^{H}
X_{⊥}^{H}

R

=

ˆL − L G

,
which is minimized when L = ˆL = X^{H}AX and

min kRk = kG k = kX_{⊥}^{H}AX k.

The proof for S is similar. If L = X^{H}AX , then

X^{H}R = X^{H}AX − X^{H}XL = X^{H}AX − L = 0.

師大

Definition

Let X be of full column rank and let X^{I} be a left inverse of X . Then
X^{I}AX is a Rayleigh quotient of A.

Theorem

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

If `1, . . . , `_{k} are the eigenvalues of L, then there are eigenvalues
λ_{j}_{1}, . . . , λ_{j}_{k} of A such that

|`_{i} − λ_{j}_{i}| ≤ kRk_{2} and
v
u
u
t

k

X

i =1

(`_{i} − λ_{j}_{i})^{2} ≤√

2kRk_{F}.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 9 / 71

師大

### Jacobi’s orthogonal component correction (JOCC), 1846

Consider the eigenvalue problem A

1 z

≡

α c^{T}
b F

1 z

= λ

1 z

, (1)

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

(1) is equivalent to

λ = α + c^{T}z,
(F − λI )z = −b.

Jacobi iteration : (with z1 = 0)

θ_{k} = α + c^{T}z_{k},

(D − θ_{k}I )z_{k+1} = (D − F )z_{k} − b (2)
where D = diag(F ).

師大

### Davidson’s method (1975)

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

Compute desired eigenpair (θ, s) ofV^{T}AV.
Computeu = Vs andr = Au − θu.

If (k r k_{2}< ε), stop.

Solve (D_{A}− θI )t = r.

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

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 11 / 71

師大

Let u_{k} = (1, z_{k}^{T})^{T}. Then

rk = (A − θkI )uk =

α − θ_{k}+ c^{T}z_{k}
(F − θ_{k}I )z_{k} + b

Substituting the residual vector r_{k} into linear systems
(DA− θ_{k}I )tk = −rk, where DA =

α 0

0 D

, we get

(D − θ_{k}I )y_{k} = −(F − θ_{k}I )z_{k} − b

= (D − F )z_{k} − (D − θ_{k}I )z_{k}− b
From (2) and above equality, we see that

(D − θ_{k}I )(z_{k}+ y_{k}) = (D − F )z_{k} − b = (D − θ_{k}I )z_{k+1}

This implies that z_{k+1} = z_{k}+ y_{k} as one step of JOCC starting with z_{k}.

師大

### Polynomial eigenvalue problems

Polynomial eigenvalue problems:

A(λ)x ≡ λ^{τ}A_{τ} + · · · + +λ^{2}A_{2}+ λA_{1}+ A_{0} x = 0. (3)
Enlarged linear eigenvalue problem: (e.g. cubic polynomial)

0 I 0

0 0 I

A_{0} A_{1} A_{2}

x
λx
λ^{2}x

= λ

I 0 0

0 I 0

0 0 −A_{3}

x
λx
λ^{2}x

. Disadvantages:

I The order of the larger matrices are tripled

I The condition number of eigenvalues and eigenvectors may increase

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 13 / 71

師大

Consider the quadratic eigenvalue problem^{1}
Q(λ)x ≡ λ^{2}M + λC + K x = 0
with

M = 1

2In, K = diag_{1≤j ≤n} j^{2}π^{2}(j^{2}π^{2}+ τ − κv^{2})/2 ,
C = −C^{T} = (c_{ij}) with c_{ij} =

( _{4ij}

j^{2}−i^{2}v if i + j is odd
0 otherwise.

Take v = 10, κ = 0.8 and τ = 77.9.

Enlarged linear eigenvalue problem:

K 0 0 I

x λx

= λ

−C −M

I 0

x λx

1F. Tisseur and K. Meerbergen, SIAM Rev. Vol. 43, pp. 235-286, 2001

師大

−1 0 1

x 10^{−9}

−4

−3

−2

−1 0 1 2 3

4x 10^{4} n = 60

polyeig exact sol

Figure: The spectrum of Q(λ)

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 15 / 71

師大

### Polynomial Jacobi-Davidson method

(θ_{k}, u_{k}): approx. eigenpair of A(λ), θ_{k} ≈ λ, with

u_{k} = V_{k}s_{k}, V_{k}^{T}A(λ)V_{k}s_{k} = 0 and k s_{k} k_{2}= 1.

Let

r_{k} = A(θ_{k})u_{k}.
Then

u_{k}^{T}r_{k} = u^{T}_{k}A(θ_{k})u_{k} = s_{k}^{T}V_{k}^{T}A(θ_{k})V_{k}s_{k} = 0 ⇒ r_{k} ⊥ u_{k}
Find the correction t such that

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

That is

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

師大

Let

p_{k} = A^{0}(θ_{k})u_{k} ≡

τ

X

i =1

i θ^{i −1}_{k} Ai

!
u_{k}.

A(λ) = A − λI :

pk = −uk,

(A(θ_{k}) − A(λ))u_{k} = (λ − θ_{k})u_{k} = (θ_{k} − λ_{k})p_{k}
A(λ) = A − λB:

pk = −Buk,

(A(θ_{k}) − A(λ))u_{k} = (λ − θ_{k})Bu_{k} = (θ_{k}− λ)p_{k}
A(λ) =Pτ

i =0λ^{i}A_{i} with τ ≥ 2:

(A(θ_{k}) − A(λ))u_{k} =

(θ_{k}− λ)A^{0}(θ_{k}) −1

2(θ_{k} − λ)^{2}A^{00}(ξ_{k})

u_{k}

= (θk − λ)p_{k} −1

2(θk − λ)^{2}A^{00}(ξk)uk

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 17 / 71

師大

Hence

A(λ)t = −rk+ (θk − λ)p_{k}− 1

2(θk− λ)^{2}A^{00}(ξk)uk

Since r_{k} ⊥ u_{k}, we have

I −p_{k}u_{k}^{T}
u_{k}^{T}pk

A(λ)t = −r_{k} −1

2(θ_{k}− λ)^{2}

I − p_{k}u_{k}^{T}
u_{k}^{T}pk

A^{00}(ξ_{k})u_{k}.
Correction equation:

I − pku^{T}_{k}
u_{k}^{T}p_{k}

A(θ_{k})(I − u_{k}u_{k}^{T})t = −r_{k} and t ⊥ u_{k},

or

I − pku_{k}^{T}
u_{k}^{T}p_{k}

(A − θ_{k}B)

I − ukp^{T}_{k}
p_{k}^{T}u_{k}

t = −r_{k} and t ⊥_{B} u_{k},

with symmetric positive definite matrix B.

師大

### Solving correction vector t

Correction equation:

I − pku^{T}_{k}
u_{k}^{T}p_{k}

A(θ_{k})(I − u_{k}u_{k}^{T})t = −r_{k}. (4)
Method I:

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

Use a preconditioner
M_{p}≡

I − p_{k}u^{T}_{k}
u_{k}^{T}pk

M

I − uku^{T}_{k}

≈

I − p_{k}u_{k}^{T}
u_{k}^{T}pk

A(θk)

I − uku_{k}^{T}

,
where M is an approximation of A(θ_{k}).

In each of the iterative steps, it needs to solve the linear system

M_{p}t = y , t ⊥ u_{k} (5)

for a given y .

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 19 / 71

師大

Since t ⊥ u_{k}, Eq. (5) can be rewritten as

I − p_{k}u^{T}_{k}
u_{k}^{T}p_{k}

Mt = y ⇒ Mt = u_{k}^{T}Mt

u_{k}^{T}p_{k} p_{k}+ y ≡ η_{k}p_{k} + y .
Hence

t = M^{−1}y + η_{k}M^{−1}p_{k},
where

ηk = − u_{k}^{T}M^{−1}y
u_{k}^{T}M^{−1}p_{k}.

SSOR preconditioner: Let A(θ_{k}) = L + D + U. Then
M = (D + ωL)D^{−1}(D + ωU).

師大

Method II: Since t ⊥ u_{k}, Eq. (4) can be rewritten as
A(θk)t = u^{T}_{k}A(θ_{k})t

u_{k}^{T}p_{k} pk − r_{k} ≡ εp_{k} − r_{k}. (6)
Let t1 and t2 be approximated solutions of the following linear

systems:

A(θ_{k})t = −r_{k} and A(θ_{k})t = p_{k},
respectively. Then the approximated solution ˜t for (6) is

˜t = t_{1}+ εt_{2} for ε = −u_{k}^{T}t_{1}
u_{k}^{T}t2

. The approximated solution ˜t for (6) is

˜t = −M^{−1}r_{k}+ εM^{−1}p_{k} for ε = u^{T}_{k}M^{−1}r_{k}
u_{k}^{T}M^{−1}p_{k},
where M is an approximation of A(θk).

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 21 / 71

師大

Method III:

Eq. (6) implies that

t = εA(θ_{k})^{−1}p_{k} − A(θ_{k})^{−1}r_{k} = εA(θ_{k})^{−1}p_{k} − u_{k}.
Let t1 be approximated solution of the following linear system:

A(θ_{k})t = p_{k}.
Then the approximated solution ˜t for (6) is

˜t = εt_{1}− u_{k} for ε =

u_{k}^{T}t_{1}−1

.

師大

Algorithm (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 of V_{k}^{T}A(λ)Vk = 0.

Select the desired (target) eigenpair (θ_{k}, s_{k}) with ks_{k}k_{2}= 1.

Computeu_{k} = V_{k}s_{k},r_{k} = A(θ_{k})u_{k} andp_{k} = A^{0}(θ_{k})u_{k}.
If (krkk_{2} < ε), λ = θk, x = uk, Stop

Solve (approximately) a t_{k} ⊥ u_{k} from
(I − ^{p}^{k}^{u}

T k

u^{T}_{k}pk)A(θk)(I − uku^{T}_{k})t = −rk.
Orthogonalize t_{k} ⊥ V_{k} → v_{k+1}, V_{k+1} = [V_{k}, v_{k+1}]

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 23 / 71

師大

### Restarting

A double precision real variable needs to 8 bytes to save in Memory.

125 double precision real variables ≈ 1 KB 125, 000 double precision real variables ≈ 1 MB

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

師大

(v) Solve correct equation (approximately) to obtain a t ⊥ u_{k}
by the method determined below.

If ( krkk_{2}> 0.1 and k ≤ 9 ) then

Use {BiCGSTAB, No precond., 7, 10^{−3}}
else

Use {GMRES, SSOR, 30, 10^{−3}}
End if

Figure: The heuristic strategy for computing the first target eigenvalue.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 25 / 71

師大

### Locking for Ax = λx

V_{k} with V_{k}^{∗}V = I_{k} are convergentSchur vectors, i.e.,
AV_{k} = V_{k}T_{k}

for someupper triangular Tk. SetV = [Vk, Vq]with V^{∗}V = Ik+q in
k + 1-th iteration of Jacobi-Davidson Algorithm. Then

V^{∗}AV =

V_{k}^{∗}AV_{k} V_{k}^{∗}AV_{q}
V_{q}^{∗}AVk V_{q}^{∗}AVq

=

T_{k} V_{k}^{∗}AV_{q}
V_{q}^{∗}VkTk V_{q}^{∗}AVq

=

T_{k} V_{k}^{∗}AV_{q}
0 V_{q}^{∗}AVq

.

師大

Locking Polynomial Jacobi–Davidson Method (0) Given A(λ) =Pτ

i =0λ^{i}A_{i} and the number of desired eigenvalues σ.

(1) Initialize V = [V_{ini}] as an orthonormal matrix and V_{x} = [].

(2) For j = 1, 2, ..., σ

(2.1) Iterate until convergence

(i) Compute thej th desired eigenpairs (θ, s) of V^{T}A(θ)V .
(ii) Compute u = Vs, p = A^{0}(θ)u, and r = A(θ)u.

(iii) If (||r ||2 < ε), Set λj = θ, xj = u, Goto Step (2.2).

(iv) Solve (approximately) a t ⊥ u from
(I − ^{pu}_{u}T^{T}p)A(θ)(I − uu^{T})t = −r .
(v) Orthogonalize t ⊥ V → v , V = [V , v ]
(2.2) Orthogonalizexj ⊥ V_{x} → x_{j}; Vx = [Vx, xj]

(2.3) Choose an orthonormal matrixVini ⊥ V_{x}; Set V = [Vx, Vini]
Ref: G. L. G. Sleijpen, G. L. Booten, D. R. Fokkema and H. A. van der
Vorst, BIT, 36:595-633, 1996

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 27 / 71

師大

(v) Solve correct equation (approximately) to obtain a t ⊥ u_{k}
by the method determined below.

If ( k r_{k} k_{2}> 0.1 and k < 10 ) then

Use {BiCGSTAB, No precond., 7, 10^{−3}}
else if ( k r_{k} k_{2}≥ 0.1 and k > 14 ) then

Use {GMRES, SSOR, 30, 10^{−3}}

else if ( k r_{k} k_{2}< 0.1 and k r_{k−1} k_{2}/ k r_{k} k_{2}< 4) then
Set j = min(30, j + 2) and use {GMRES, SSOR, j, 10^{−3}}
else

Use {GMRES, SSOR, j, 10^{−3}}
End if

Figure: The heuristic strategy for computing eigenvalues other than the first convergent eigenvalue.

師大

### Non-equivalence deflation of quadratic eigenproblems

Let λ1 be a real eigenvalue of Q(λ) and x1, z1 be the associated right and
left eigenvectors, respectively, with z_{1}^{T}Kx_{1} = 1. Let

θ1 = (z_{1}^{T}Mx1)^{−1}.
We introduce a deflated quadratic eigenproblem

Q(λ)x ≡e h

λ^{2}M + λ ee C + eKi
x = 0,
where

Me = M − θ_{1}Mx_{1}z_{1}^{T}M,
Ce = C + θ1

λ1

(Mx1z_{1}^{T}K + Kx1z_{1}^{T}M),
Ke = K − θ1

λ^{2}_{1}Kx1z_{1}^{T}K .

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 29 / 71

師大

### Complex deflation

Let λ_{1} = α_{1}+ i β_{1} be a complex eigenvalue of Q(λ) and x_{1} = x_{1R}+ ix_{1I},
z_{1} = z_{1R}+ iz_{1I} be the associated right and left eigenvectors, respectively,
such that

Z_{1}^{T}KX1 = I2,
where X1 = [x1R, x1I] and Z1 = [z1R, z1I]. Let

Θ_{1} = (Z_{1}^{T}MX_{1})^{−1}.

Then we introduce a deflated quadratic eigenproblem with
Me = M − MX_{1}Θ_{1}Z_{1}^{T}M,

Ce = C + MX_{1}Θ_{1}Λ^{−T}_{1} Z_{1}^{T}K + KX_{1}Λ^{−1}_{1} Θ^{T}_{1}Z_{1}^{T}M,
Ke = K − KX1Λ^{−1}_{1} Θ1Λ^{−T}_{1} Z_{1}^{T}K

in which Λ_{1} =

α_{1} β_{1}

−β_{1} α_{1}

.

師大

Theorem

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

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

provided that λ^{2}_{1} 6= θ_{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}Λ^{T}_{1} 6= Θ_{1}.

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

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

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 31 / 71

師大

Suppose that M, C , K are symmetric. Given an eigenmatrix pair

(Λ1, X1) ∈ R^{r ×r} × R^{n×r} of Q(λ), where Λ1 is nonsingular and X1 satisfies
X_{1}^{T}KX1 = Ir, Θ1:= (X_{1}^{T}MX1)^{−1}.

We define ˜Q(λ) := λ^{2}M + λ ˜˜ C + ˜K , where
M˜ := M − MX_{1}Θ_{1}X_{1}^{T}M,

C˜ := C + MX_{1}Θ_{1}Λ^{−T}_{1} X_{1}^{T}K + KX_{1}Λ^{−1}_{1} Θ_{1}X_{1}^{T}M,
K˜ := K − KX1Λ^{−1}_{1} Θ1Λ^{−T}_{1} X_{1}^{T}K .

Theorem

Suppose that Θ_{1}− Λ_{1}Λ^{T}_{1} 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.

師大

### Proof:

Since (Λ_{1}, X_{1}) is an eigenmatrix pair of Q(λ), i.e.,
MX_{1}Λ^{2}_{1}+ CX_{1}Λ_{1}+ KX_{1} = 0,
we have

Q(λ)˜ = Q(λ) + [MX1(λIr + Λ1) + CX1] Θ1Λ^{−T}_{1} (X_{1}^{T}K − λΛ^{T}_{1}X_{1}^{T}M)

= Q(λ) + Q(λ)X1(λIr − Λ_{1})^{−1}Θ1Λ^{−T}_{1} (X_{1}^{T}K − λΛ^{T}_{1}X_{1}^{T}M).

By using the identity

det(I_{n}+ RS ) = det(I_{m}+ SR),
where R, S^{T} ∈ R^{n×m}, we have

det[ ˜Q(λ)]

= det[Q(λ)]det[I + X1(λIr − Λ_{1})^{−1}Θ1Λ^{−T}_{1} (X_{1}^{T}K − λΛ^{T}_{1}X_{1}^{T}M)]

= det[Q(λ)]det[Ir + (λIr − Λ_{1})^{−1}Θ1Λ^{−T}_{1} (Ir − λΛ^{T}_{1}Θ^{−1}_{1} )]

= det[Q(λ)]

det(λIr − Λ_{1})det(Θ1Λ^{−T}_{1} − Λ_{1}).

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 33 / 71

師大

Since (Θ1− Λ_{1}Λ^{T}_{1}) ∈ R^{r ×r} is nonsingular, we have
det(Θ_{1}Λ^{−T}_{1} − Λ_{1}) 6= 0.

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

師大

### Non-equivalence deflation for cubic polynomial eigenproblems

Let (Λ, V_{u}) ∈ R^{r ×r} × R^{n×r} be an eigenmatrix pair of A(λ) with
V_{u}^{T}Vu= Ir and 0 /∈ σ(Λ), i.e.,

A_{3}V_{u}Λ^{3}+ A_{2}V_{u}Λ^{2}+ A_{1}V_{u}Λ + A_{0}V_{u}= 0. (7)
Define a new deflated cubic eigenvalue problem by

A(λ)u = (λe ^{3}A˜3+ λ^{2}A˜2+ λ˜A1+ ˜A0)u = 0, (8)
where

A˜0 = A0,

A˜_{1} = A_{1}− (A_{1}V_{u}V_{u}^{T}+ A_{2}V_{u}ΛV_{u}^{T} + A_{3}V_{u}Λ^{2}V_{u}^{T}),
A˜_{2} = A_{2}− (A_{2}V_{u}V_{u}^{T}+ A_{3}V_{u}ΛV_{u}^{T}),

A˜3 = A3− A_{3}VuV_{u}^{T}.

(9)

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 35 / 71

師大

Lemma

Let A(λ) and eA(λ) be cubic pencils given by (3) and (8), respectively.

Then it holds

A(λ) = A(λ)e

I_{n}− λV_{u}(λI_{r} − Λ)^{−1}V_{u}^{T}

. (10)

Theorem

Let (Λ, V_{u}) be an eigenmatrix pair of A(λ) with V_{u}^{T}V_{u}= I_{r}. Then
(i) (σ(A(λ))σ(Λ)) ∪ {∞} = σ(eA(λ)).

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

Define

˜

z = (I_{n}− µV_{u}Λ^{−1}V_{u}^{T})z ≡ T (µ)z. (11)
Then (µ, ˜z) is an eigenpair of eA(λ).

師大

### Proof of Lemma:

Using (9) and (7), and the fundamental matrix calculation, we have A(λ)e = A(λ) − λ

λ^{2}A3VFV_{F}^{T} + λA2VFV_{F}^{T} + λA3VFΛV_{F}^{T}+ A1VFV_{F}^{T}
+A2V_{F}ΛV_{F}^{T} + A3V_{F}Λ^{2}V_{F}^{T}

= A(λ) − λ

A_{3}V_{F}(λI_{r} − Λ)^{3}(λI_{r} − Λ)^{−1}V_{F}^{T}
+3A_{3}V_{F}Λ(λI_{r} − Λ)^{2}(λI_{r} − Λ)^{−1}V_{F}^{T}
+3A3V_{F}Λ^{2}(λIr − Λ)(λI_{r} − Λ)^{−1}V_{F}^{T}
+A2VF(λIr − Λ)^{2}(λIr − Λ)^{−1}V_{F}^{T}
+2A_{2}V_{F}Λ(λI_{r} − Λ)(λI_{r} − Λ)^{−1}V_{F}^{T}

+A_{1}V_{F}(λI_{r} − Λ)(λI_{r} − Λ)^{−1}V_{F}^{T}

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 37 / 71

師大

A(λ)e = A(λ) − λA_{3}V_{F}(λ^{3}I_{r} − Λ^{3}) + A_{2}V_{F}(λ^{2}I_{r} − Λ^{2})
+A1VF(λIr − Λ) + A_{0}VF − A_{0}VF] (λIr − Λ)^{−1}V_{F}^{T}

o

= A(λ) − λh

A(λ)V_{F}(λI_{r} − Λ)^{−1}V_{F}^{T}i

= A(λ)h

I_{n}− λV_{F}(λI_{r} − Λ)^{−1}V_{F}^{T}i
.

師大

### Proof of Theorem

: (i) Using the identity

det(I_{n}+ RS ) = det(I_{m}+ SR)
and Lemma 8, we have

det( eA(λ)) = det(A(λ)) det

In− λV_{F}(λIr − Λ)^{−1}V_{F}^{T}

= det(A(λ)) det In− λ(λI_{r} − Λ)^{−1}

= det(A(λ)) det(λI_{r} − Λ)^{−1}det(−Λ).

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

Ae3V_{F} = (A3− A_{3}V_{F}V_{F}^{T})V_{F} = 0,

we see that (diag_{r}{∞, · · · , ∞}, V_{F}) is an eigenmatrix pair of eA(λ)
corresponding to infinite eigenvalues.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 39 / 71

師大

(ii) Since µ /∈ σ(Λ), the matrix T (µ) = (I − µV_{F}Λ^{−1}V_{F}^{T}) in (11) is
invertible with the inverse

T (µ)^{−1}= I_{n}− µV_{F}(µI_{r} − Λ)^{−1}V_{F}^{T}. (12)
From Lemma 8, we have

A(µ)˜e z = A(µ) h

In− µV_{F}(µIr − Λ)^{−1}V_{F}^{T}
i h

In− µV_{F}Λ^{−1}V_{F}^{T}
i

z = 0.

This completes the proof.

師大

### Applications

Quantum well (2 dim.)

Quantum wire (1 dim.)

Quantum dot (0 dim.)

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 41 / 71

師大

### Nanometer

1nm = 10^{−9}m

Nano-scale ≈ 1–100 nm A semiconductor QD ≈ 10 nm QD : hair ≈ 1 : 10000

Why consider quantum effects?

Small devices imply significant quantum effect

師大

### Quantum dot

Cross-sections of hetero-structure InAs/GaAs QDs by Transmission Electron Microscope [Schoenfeld, 00]

**150 Å**

**300 Å**
**150 Å**

**300 Å**
**150 Å**

**300 Å**
**150 Å**

**300 Å**

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 43 / 71

師大

### Molecular beam epitaxy

師大

### Quantum dot

Cross-sections of hetero-structure InAs/GaAs QDs by Transmission Electron Microscope [Schoenfeld, 00]

**150 Å**

**300 Å**
**150 Å**

**300 Å**
**150 Å**

**300 Å**
**150 Å**

**300 Å**

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 45 / 71

師大

### Energy levels (eigenvalues)

Confined Discrete Energy Levels Confined Discrete

Energy Levels

師大

### Numerical experiments for linear eigenproblems

TheSchr¨odinger equationfor Semiconductor:

−∇ · (α∇u) + Vu = λu,

where α =

( α^{−}≡ _{2m}^{~}^{2}

1 inside,
α^{+}≡ _{2m}^{~}^{2}

2 outside, V =

V^{−}= V1 inside,
V^{+}= V_{2} outside

~: Plank constant m_{`}: parabolic effective mass
V_{`}: confinement potential λ: total energy

Interface condition:

α^{−}∂u

∂n ∂D−

= α^{+}∂u

∂n ∂D+

Dirichlet boundary conditions

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 47 / 71

師大

Symmetric eigenvalue problem

Ax = λx where A is a symmetric matrix.

Reference:

Tsung-Min Hwang, Wen-Wei Lin, Wei-Cheng Wang and Weichung Wang, Numerical simulation of three dimensional pyramid quantum dot, Journal of Computational Physics, Vol 196, pp. 208-232, 2004.

Figure: PRB, 54, 8743, (1996)

師大

### Three dimensional pyramid quantum dot

H

W W

Pyramid Dot Cuboid Matrix i

2 1. 5 1 0. 5 0 0.5 1 1.5 2

1 0. 5 0 0.5 1 1.5 2

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 2 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 2 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 2 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 1 2 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 1 1 2 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 1 1 1 5 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 1 1 3 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 1 3 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 1 3 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 1 3 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 1 3 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 4 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 49 / 71

師大

Finite volume discretized scheme Symmetric eigenvalue problems:

Ax = λx

Second order convergent rate Correction vector with SSOR preconditioner:

t = −M_{A}^{−1}r + εM_{A}^{−1}p

with

ε = u^{T}M_{A}^{−1}r
u^{T}M_{A}^{−1}p,

whereM_{A}= (D + ωL)D^{−1}(D + ωU).

H

W W

Pyramid Dot Cuboid Matrix i

Figure: Structure schema of a pyramid quantum dot.

師大

The width of the QD (matrix) base is 12.4 nm (24.8 nm); the height of the QD (matrix) is 6.2 nm (18.6 nm).

InAs QD: α1 = 0.024me and V1 = 0.0 GaAs matrix: α2= 0.067me and V2= 0.70.

Stopping criteria: residual < 10^{−10}
Convergent rate

(L,M,N) Mtx. dim. λ_{1} Rate λ_{2} Rate

(16, 16, 12) 2,475 0.4226 - 0.6527 -
(32, 32, 24) 22,103 0.4001 - 0.6423 -
(64, 64, 48) 186,543 0.3934 1.744 0.6391 1.708
(128,128, 96) 1,532,255 0.3916 1.905 0.6383 1.866
(256,256,192) 12,419,775 0.3911 1.954 0.6380 1.912
Table: Convergent rate =log_{2} (λ^{(4h)}− λ^{(2h)})/(λ^{(2h)}− λ^{(h)})

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 51 / 71

師大

dim(A) = 1, 532, 255on PC with P4 1.8 GHz CPUand1 GB of main memory.

Value Ite. no. CPU time (sec.)

λ1 0.3916 72 278.0

λ_{2} 0.6383 72 284.3

λ3 0.6383 125 521.4

dim(A) = 32, 401, 863on HP workstation with a 1.0 GHz Intel Itanium II CPUand 12 GBsof main memory.

Value Ite. no. CPU time (sec.)

λ_{1} 0.3910 138 5,852

λ2 0.6380 133 5,354

λ3 0.6380 220 8,511

師大

Unsymmetric eigenvalue problem

Ax = λx where A is a unsymmetric matrix.

Reference:

Tsung-Min Hwang, Wei-Hua Wang and Weichung Wang, Efficient numerical schemes for electronic states in coupled quantum dots, accepted for publication in Journal of Nanoscience and

Nanotechnology.

JAP, 90-12, (2001)-

Figure: JAP, 90-12, (2001)

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 53 / 71

師大

### Vertically aligned quantum dot array

(a) Nonuniform mes h (b) Uniform mesh

師大

Reduce three-dimensional systems to two-dimensional systems by Fourier transformation.

Finite volume discretized scheme Unsymmetric eigenvalue problems:

Ax = λx

Second order convergent rate

Method II with SSOR preconditioner:

M_{A}= (D + ωL)D^{−1}(D + ωU).

Figure: Structure schema of a cylindrical vertically aligned quantum dot array.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 55 / 71

師大

Reduce three-dimensional systems to two-dimensional systems by Fourier transformation.

Finite volume discretized scheme Unsymmetric eigenvalue problems:

Ax = λx

Second order convergent rate

Method II with SSOR preconditioner:

M_{A}= (D + ωL)D^{−1}(D + ωU). _{Figure:} Structure schema of a
cylindrical vertically aligned
quantum dot array.

師大

Reduce three-dimensional systems to two-dimensional systems by Fourier transformation.

Finite volume discretized scheme Unsymmetric eigenvalue problems:

Ax = λx

Second order convergent rate

Method II with SSOR preconditioner:

M_{A}= (D + ωL)D^{−1}(D + ωU).

Figure: Structure schema of a cylindrical vertically aligned quantum dot array.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 55 / 71

師大

gap = 6nm and dim(A) = 12, 288, 000on HP workstation with a 1.0 GHz Intel Itanium II CPU and12 GBsof main memory.

r2 = 3.198 r2 = 3.223

Value Ite. no. Time (sec.) Value Ite. no. Time (sec.)

λ1 0.1587 83 20,287 0.1587 84 23,840

λ_{2} 0.35553 109 26,030 0.3526 91 31,612

λ3 0.3558 53 13,831 0.3555 61 25,706

gap = 3nm and dim(A) = 11, 059, 200on HP workstation with a 1.0 GHz Intel Itanium II CPU and12 GBsof main memory.

r2 = 3.198 r2 = 3.223

Value Ite. no. Time (sec.) Value Ite. no. Time (sec.)

λ_{1} 0.1586 85 20,072 0.1586 83 18,455

λ_{2} 0.3540 257 56,811 0.3515 130 29,506

λ3 0.3564 53 12,257 0.3558 54 13,082

師大

Generalized eigenvalue problem

Ax = λBx

where A is symmetric positive definiteand B is a positive diagonalmatrix.

Reference:

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 in Journal of Computational Physics.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 57 / 71

師大

### Three dimensional arbitrary shape quantum dots

Quantum dot Matrix

Z_{A}
ZB

ZC

1 M

Rmtx

0 H

N1

N2

N

H HC

HRmtx

R adial coordinate

Axial coordinate

師大

Curvilinear coordinate system

Reduce three-dimensional systems to two-dimensional systems by Fourier transformation.

Jump condition capturing scheme (nine points)

Generalized eigenvalue problems:

Ax = λBx , A issymmetric positive definite and B is apositive diagonal matrix.

Second order convergent rate

Method I with SSOR preconditioner:

M_{A}= (D + ωL)D^{−1}(D + ωU).

**Quantum dot**
**Matrix**

Figure: Structure schema of the quantum dot model.

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 59 / 71