師大

### 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**

師大

### 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**

師大

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**

師大

### 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 = 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.)** **Poly. JD method for PEP** **April 28, 2011** **4 / 42**

師大

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

yi = 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.

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

師大

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

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

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

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

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

師大

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}AXand

min kRk = kGk = kX_{⊥}^{H}AXk.

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.

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

師大

Definition 3

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 4

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

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

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

k

X

i=1

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

2kRkF.

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

師大

### Jacobi’s orthogonal component correction, 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 z_{1} = 0)

θk= α + c^{T}zk,

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

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

師大

### Davidson’s method (1975)

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

Iterate until convergence

Compute desired eigenpair (θ, s) ofV^{T}AV.
Computeu = V sandr = 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.)** **Poly. JD method for PEP** **April 28, 2011** **10 / 42**

師大

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)t_{k}= −r_{k}, where DA=

α 0

0 D

, we get

(D − θkI)yk = −(F − θ_{k}I)zk− b

= (D − F )zk− (D − θ_{k}I)zk− 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}.

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

師大

### Jacobi-Davidson method (1996)

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

u_{k}= V_{k}s_{k}, V_{k}^{T}AV_{k}s_{k}= θ_{k}s_{k} and k s_{k}k_{2}= 1.

Definition 5

(θk, uk)is called a Ritz pair of A. θ_{k}is called a Ritz value and u_{k}is a
Ritz vector.

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

師大

### Jacobi-Davidson method (1996)

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

u_{k}= V_{k}s_{k}, V_{k}^{T}AV_{k}s_{k}= θ_{k}s_{k} and k s_{k}k_{2}= 1.

Then

u^{T}_{k}r_{k} ≡ u^{T}_{k}(A − θ_{k}I)u_{k} = s^{T}_{k}V_{k}^{T}AV_{k}s_{k}− θ_{k}s^{T}_{k}V_{k}^{T}V_{k}s_{k} = 0 ⇒ r_{k}⊥ u_{k}
Find the correctiont ⊥ uksuch that

A(u_{k}+ t) = λ(u_{k}+ t).

That is

(A − λI)t = λu_{k}− Au_{k}= −r_{k}+ (λ − θ_{k})u_{k}.

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

師大

Since (I − u_{k}u^{T}_{k})rk= rk, we have

I − u_{k}u^{T}_{k} (A − λI)t = −r_{k}.

Correction equation

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

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

師大

### Solving correction vector t

Correction equation:

I − uku^{T}_{k} (A − θ_{k}I)(I − uku^{T}_{k})t = −rk, t ⊥ uk. (3)
Scheme S_{OneLS}:

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

Use a preconditioner

M_{p} ≡ I − u_{k}u^{T}_{k} M I − u_{k}u^{T}_{k} ,
where M is an approximation of A − θ_{k}I.
In each of the iterative steps, it needs to solve

M_{p}y = b, y ⊥ u_{k} (4)

for a given b.

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

師大

Since y ⊥ u_{k}, Eq. (4) can be rewritten as

I − uku^{T}_{k} My = b ⇒ My = u^{T}_{k}My u_{k}+ b ≡ ηkuk+ b.

Hence

y = M^{−1}b + η_{k}M^{−1}u_{k},
where

ηk= − u^{T}_{k}M^{−1}b
u^{T}_{k}M^{−1}uk

.

SSOR preconditioner: Let A − θ_{k}I = 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**

師大

Scheme S_{T woLS}: Since t ⊥ u_{k}, Eq. (3) can be rewritten as

(A − θ_{k}I)t = u^{T}_{k}(A − θ_{k}I)t u_{k}− r_{k}≡ εu_{k}− r_{k}. (5)
Let t_{1}and t_{2} be approximated solutions of the following linear systems:

(A − θ_{k}I)t = −r_{k} and (A − θ_{k}I)t = u_{k},
respectively. Then the approximated solution ˜tfor (5) is

˜t = t1+ εt2 for ε = −u^{T}_{k}t_{1}
u^{T}_{k}t2

.
Scheme S_{OneStep}: The approximated solution ˜tfor (5) is

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

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

師大

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

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

Compute all the eigenpairs ofV_{k}^{T}AVks = λs.

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

Computeu_{k} = V_{k}s_{k}andr_{k} = (A − θ_{k}I)u_{k}.
If (kr_{k}k_{2}< ε), λ = θ_{k}, x = u_{k}, Stop

Solve (approximately) at_{k}⊥ u_{k} from

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

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

師大

Locking:

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

for someupper triangularT_{k}. SetV = [V_{k}, V_{q}]with V^{∗}V = I_{k+q} in
k + 1-th iteration of Jacobi-Davidson Algorithm. Then

V^{∗}AV =

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

=

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

=

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

. 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**

師大

Remark 1

If ε = 0, we obtain Davidson method with
t_{1} = −(D_{A}− θ_{k}I)^{−1}r.

( ˜tis NOT orthogonal to u_{k} )

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

t = ε(A − θ_{k}I)^{−1}u_{k}− u_{k}. (6)
Since t is made orthogonal to u_{k}, (6) is equivalent to

t = (A − θkI)^{−1}uk

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**

師大

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

Lemma 6

Consider w with w^{T}x 6= 0. Then the map
Fp ≡

I −xw^{T}
w^{T}x

(A − λI)

I −xw^{T}
w^{T}x

is a bijection from w^{⊥} to w^{⊥}.

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

I − xw^{T}
w^{T}x

(A − λI)

I −xw^{T}
w^{T}x

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**

師大

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 x^{T}w 6= 0implies y = 0, which proves the injectivity of Fp. An
obvious dimension argument implies bijectivity.

Extension
F_{p}t ≡

I − uu^{T}
u^{T}u

(A − θI)

I − uu^{T}
u^{T}u

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**

師大

Theorem 7

Assume that the correction equation

I − uu^{T} (A − θI) I − uu^{T} t = −r, t ⊥ u (7)
is solved exactly in each step of Jacobi-Davidson Algorithm. Assume
uk= u → xand u^{T}_{k}xhas non-trivial limit. Then if u_{k}is sufficiently
chosen to x, then u_{k}→ x with locally quadratical convergence and

θk= u^{T}_{k}Auk → λ.

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

師大

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 z_{1} ⊥ u of (7):

(I − P ) (A − θI) z_{1} = − (I − P ) r, (9)
where P = uu^{T}. Note that (I − P ) r = r since u⊥r. Since

x − (u + z_{1}) = z − z_{1} and z = x − u, for quadratic convergence, it
suffices to show that

kx − (u + z_{1}) k = kz − z1k = O kzk^{2} . (10)
Multiplying (8) by (I − P ) and subtracting the result from (9) yields

(I − P ) (A − θI) (z − z_{1}) = (λ − θ) (I − P ) z + (λ − θ) (I − P ) u. (11)

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

師大

Multiplying (8) by u^{T} and using r ⊥ u leads to
λ − θ = u^{T}(A − θI) z

u^{T}x . (12)

Since u^{T}_{k}xhas non-trivial limit, we obtain
k (λ − θ) (I − P ) zk =

u^{T} (A − θI) z

u^{T}x (I − P ) z

. (13)

From (11), Lemma 6 and (I − P ) u = 0, we have
kz − z_{1}k =

h

(I − P ) (A − θkI) |_{u}^{⊥}

k

i−1

(λ − θ) (I − P ) z

= h

(I − P ) (A − θkI) |_{u}^{⊥}

k

i−1 u^{T}_{k} (A − θ_{k}I) z

u^{T}_{k}x (I − P ) z

= O kzk^{2} .

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

師大

### Jacobi Davidson method as on accelerated Newton Scheme

Consider Ax = λx and assume that λ is simple. Choose w^{T}x = 1.

Consider nonlinear equation

F (u) = Au − θ(u)u = 0 with θ (u) = w^{T}Au
w^{T}u ,

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

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

∂u u=uk

!−1

F (u_{k})
is given by u_{k+1}= u_{k}+ t, i.e., t satisfies that

∂F

∂u
u=u_{k}

!

t = F (u_{k}) = −r.

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

師大

Since 1 = u^{T}_{k+1}w = (u_{k}+ t)^{T}w = 1 + t^{T}w, it implies that w^{T}t = 0. By
the definition of F , we have

∂F

∂u = A − θ (u) I −− w^{T}Au uw^{T} + w^{T}u uw^{T}A
(w^{T}u)^{2}

= A − θI + w^{T}Au

(w^{T}u)^{2}uw^{T} −uw^{T}A
w^{T}u =

I −u_{k}w^{T}
w^{T}u_{k}

(A − θ_{k}I) .
Consequently, the Jacobian of F acts on w^{⊥}and is given by

∂F

∂u u=uk

! t =

I − u_{k}w^{T}
w^{T}u_{k}

(A − θ_{k}I) t, t ⊥ w.

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

I − ukw^{T}
w^{T}u_{k}

(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**

師大

### Polynomial Jacobi-Davidson method

(θ_{k}, u_{k}): approx. eigenpair of A(λ) ≡Pτ

k=0λ^{k}A_{k}, θ_{k}≈ λ, with
u_{k}= V_{k}s_{k}, V_{k}^{>}A(λ)V_{k}s_{k}= 0 and k s_{k}k_{2}= 1.

Let

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

u^{>}_{k}rk= u^{>}_{k}A(θk)uk= s^{>}_{k}V_{k}^{>}A(θk)Vksk= 0 ⇒ rk ⊥ u_{k}
Find the correctiontsuch that

A(λ)(u_{k}+ t) = 0.

That is

A(λ)t = −A(λ)u_{k} = −r_{k}+ (A(θ_{k}) − A(λ))u_{k}.

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

師大

Let

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

τ

X

i=1

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

!
u_{k}.

A(λ) = A − λI:
p_{k}= −u_{k},

(A(θk) − A(λ))uk = (λ − θk)uk= (θk− λ_{k})pk

A(λ) = A − λB:

p_{k}= −Bu_{k},

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

i=0λ^{i}Aiwith τ ≥ 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.)** **Poly. JD method for PEP** **April 28, 2011** **28 / 42**

師大

Hence

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

2(θ_{k}− λ)^{2}A^{00}(ξ_{k})u_{k}
Since r_{k}⊥ u_{k}, we have

I − p_{k}u^{>}_{k}
u^{>}_{k}pk

A(λ)t = −rk−1

2(θk− λ)^{2}

I − p_{k}u^{>}_{k}
u^{>}_{k}pk

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

I − pku^{>}_{k}
u^{>}_{k}p_{k}

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

or

I −pku^{>}_{k}
u^{>}_{k}p_{k}

(A − θ_{k}B)

I −ukp^{>}_{k}
p^{>}_{k}u_{k}

t = −r_{k} and t ⊥_{B}u_{k},
with symmetric positive definite matrix B.

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

師大

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

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

Compute all the eigenpairs ofV_{k}^{>}A(λ)V_{k} = 0.

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

Computeuk = Vksk,rk= A(θk)ukandpk= A^{0}(θk)uk.
If (kr_{k}k_{2}< ε), λ = θ_{k}, x = u_{k}, Stop

Solve (approximately) at_{k}⊥ u_{k} from
(I −^{p}_{u}^{k}T^{u}^{T}^{k}

kpk)A(θ_{k})(I − u_{k}u^{T}_{k})t = −r_{k}.
Orthogonalize t_{k}⊥ V_{k}→ v_{k+1}, V_{k+1}= [Vk, vk+1]

pk = A^{0}(θk)uk≡

τ

X

i=1

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

! uk.

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

師大

### Non-equivalence deflation of quadratic eigenproblems

Let λ1 be a real eigenvalue of Q(λ) ≡ λ^{2}M + λC + K and x1, z1be the
associated right and left eigenvectors, respectively, with z^{T}_{1}Kx_{1}= 1.

Let

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

Q(λ)x ≡e h

λ^{2}M + λ ef C + eK
i

x = 0, where

Mf = M − θ1M x1z_{1}^{T}M,
Ce = C + θ_{1}

λ_{1}(M x_{1}z_{1}^{T}K + Kx_{1}z_{1}^{T}M ),
Ke = K − θ_{1}

λ^{2}_{1}Kx_{1}z^{T}_{1}K.

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

師大

### Complex deflation

Let λ1 = α1+ iβ1be 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}M X_{1})^{−1}.

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

Ce = C + M X_{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}

.

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

師大

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 λ^{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}, x_{2})is also an eigenpair of eQ(λ).

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

師大

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 X_{1}
satisfies

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

M˜ := M − M X_{1}Θ_{1}X_{1}^{T}M,

C˜ := C + M X_{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 9

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.

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

師大

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

we have

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

= Q(λ) + Q(λ)X_{1}(λI_{r}− Λ_{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[I_{r}+ (λI_{r}− Λ_{1})^{−1}Θ_{1}Λ^{−T}_{1} (I_{r}− λΛ^{T}_{1}Θ^{−1}_{1} )]

= det[Q(λ)]

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

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

師大

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.

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

師大

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

Let (Λ, Vu) ∈ R^{r×r}× R^{n×r} be an eigenmatrix pair of

A(λ) ≡ λ^{3}A_{3}+ λ^{2}A_{2}+ λA_{1}+ A_{0} (14)
with V_{u}^{T}V_{u}= I_{r}and 0 /∈ σ(Λ), i.e.,

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

A(λ)u = (λe ^{3}A˜3+ λ^{2}A˜2+ λ ˜A1+ ˜A0)u = 0, (16)
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} = A_{3}− A_{3}V_{u}V_{u}^{T}.

(17)

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

師大

Lemma 10

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

A(λ) = A(λ) Ie n− λV_{u}(λIr− Λ)^{−1}V_{u}^{T} . (18)
Theorem 11

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 k2= 1and µ /∈ σ(Λ).

Define

˜

z = (In− µV_{u}Λ^{−1}V_{u}^{T})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**

師大

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

A(λ)e = A(λ) − λ λ^{2}A3VFV_{F}^{T} + λA2VFV_{F}^{T} + λA3VFΛV_{F}^{T}
+A_{1}V_{F}V_{F}^{T} + A_{2}V_{F}ΛV_{F}^{T} + A_{3}V_{F}Λ^{2}V_{F}^{T}

= A(λ) − λ A3VF(λIr− Λ)^{3}(λIr− Λ)^{−1}V_{F}^{T}
+3A3VFΛ(λIr− Λ)^{2}(λIr− Λ)^{−1}V_{F}^{T}
+3A_{3}V_{F}Λ^{2}(λI_{r}− Λ)(λI_{r}− Λ)^{−1}V_{F}^{T}
+A_{2}V_{F}(λI_{r}− Λ)^{2}(λI_{r}− Λ)^{−1}V_{F}^{T}
+2A2V_{F}Λ(λIr− Λ)(λI_{r}− Λ)^{−1}V_{F}^{T}

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

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

師大

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

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

= A(λ)I_{n}− λV_{F}(λIr− Λ)^{−1}V_{F}^{T} .

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

師大

Proof of Theorem: (i) Using the identity

det(I_{n}+ RS) = det(I_{m}+ SR)
and Lemma 10, 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(λIr− Λ)^{−1}det(−Λ).

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− A_{3}VFV_{F}^{T})VF = 0,

we see that (diag_{r}{∞, · · · , ∞}, V_{F})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**

師大

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

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

A(µ)˜e z = A(µ)I_{n}− µV_{F}(µIr− Λ)^{−1}V_{F}^{T} I_{n}− µV_{F}Λ^{−1}V_{F}^{T} z = 0.

This completes the proof.

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