Structured Algorithms for Palindromic Quadratic Eigenvalue Problems - - PowerPoint PPT Presentation

Structured Algorithms for Palindromic Quadratic Eigenvalue Problems : Vibration of Fast Trains

Wen-Wei Lin

Department of Mathematics Tsing Hua University, Taiwan A joint work with J. Qian and T.-M. Huang

Jan., 2008

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 1 / 54

Outline

1 Palindromic quadratic eigenvalue problem 2 Structured Algorithm I 3 Structured Algorithm II vs. URV-based structured Algorithm 4 T-skew-Hamiltonian Arnoldi algorithm 5 Generalized ⊤-skew-Hamiltonian Arnoldi Algorithm 6 Numerical results 7 Conclusions

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 2 / 54

Resonances in rail tracks excited by high speed trains

With new ICE trains crossing Europe at speeds of up to 300 km/h, sound and vibration levels in the trains are an important issue. Hilliges/Mehrmann/Mehl(2004) first proposed this problem on a project with company SFE GmbH in Berlin.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 3 / 54

Finite Element Model

(a). The rails are straight and infinitely long: A 3D finite element discretization of the rail with linear isoparametric tetrahedron elements(Chu/Huang/Lin/Wu, JCAM, 2007) produces an infinite-dimensional system of O.D.E.: M ¨ x + D ˙ x + Kx = F, where M, D and K are block-tridiagonal matrices.

Figure: A 3D rail model.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 4 / 54

(b). The rail sections of track between two ties are identical: The system is periodic and leads to a Palindromic QEP P(λ)x ≡ (λ2A⊤

1 + λA0 + A1)x = 0,

(1) where A0, A1 ∈ Cn×n, and A⊤

0 = A0.

◮ The coefficient matrices A0 and A1 in (1) depend on some ω

associated with the excitation frequency of the external force.

◮ Eq. (1) is called a palindromic QEP problem.

Symplectic property

If λ is an eigenvalue of P(λ), then so is λ−1.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 5 / 54

Numerical difficulties

Need to compute all finite, nonzero eigenvalues and eigenvectors for all ω; Lots of eigenvalues are 0 and ∞; Problem size 1,000–100,000; Range of eigenvalues |λ| ∈ [10−20, 1020]; Problem is badly scaled;

20 40 60 80 100 120 10-20 10

• 15

10

• 10

10

• 5

10 10

5

10

10

10

15

10

20

| λ |

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 6 / 54

Find all eigenpairs of a 2n × 2n linearized form (M − λL)z ≡

• A1

−A0 −I

• − λ
• I

A⊤

1

x λx

• = 0

by using QZ algorithm (standard approach). ⇒The symplectic structure is not preserved which produces large numerical errors.

• D. Mackey/N. Mackey/Mehl/Mehrmann (SIMAX, 2006 (2)) linearize
• Eq. (1) to a palindromic linear pencil λZ ± Z⊤.

Hilliges/Mehl/Mehrmann (2005) propose a hybrid method (Laub’s trick+ Jacobi-like method) to solve λZ ± Z⊤. ⇒ It work well if there are no EWs near ±1. and needs O(n3 log(n)) flops.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 7 / 54

Schr¨

• der (2001) proposed an URV based structured method for

solving λZ ± Z⊤. ⇒ It is equivalent to the structured Algorithm I (SA I) applying to µ2Z⊤ + µ0 + Z with µ2 = λ. Chu/Huang/L/Wu (JCAM, 2007) proposed a structure-preserving doubling alg. for solving NME: X + A⊤

1 X−1A1 = −A0 associated

with (1). ⇒ The numerical results show much promise, but the convergence theory holds under assumption that the sequence by SDA is well-defined.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 8 / 54

Consider the linearization: (M − λL)z ≡ A1 −A0 −I

• − λ

I A⊤

1

x λx

• = 0.

(2) M and L satisfy MJ M⊤ = LJ L⊤ with J = In −In

• .

(3) M − λL is called a ⊤-symplectic pencil. If ν is an eigenvalue of M − λL, so is 1/ν.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 9 / 54

Structured Algorithm I (SA I)

(S + S−1)-transform for ⊤-symplectic M − λL (L, LAA, 1987): Ms ≡ MJ L⊤ + LJ M⊤, Ls ≡ LJ L⊤. (4) µ is a double eigenvalue of Ms − λLs ⇔ ν, 1

ν are eigenvalues of

M − λL, where ν, 1

ν are two roots of λ + 1 λ = µ.

Theorem

Let M − λL be the T-symplectic pencil of (2). If zs = [z⊤

1 , z⊤ 2 ]⊤ is an EV

• f Ms − λLs corresp. to µ = ν + 1

ν (µ = ±2), then z1 + 1 ν z2 and z1 + νz2

are EVs of P(λ) in (1) corresp. to ν and 1

ν , resp..

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 10 / 54

From (2), we have Ms − λLs = A1 − A⊤

1

A0 −A0 A1 − A⊤

1

• − λ

−A1 A⊤

1

• =
• A0

A⊤

1 − A1

A1 − A⊤

1

A0

• − λ

−A1 −A⊤

1

• J

≡ (K − λN)J , (5) where K, N are both ⊤-skew-Hamiltonian, i.e., (KJ )⊤ = −KJ and (NJ )⊤ = −NJ . Note that a matrix S is ⊤-skew-Hamiltonian, then S = S11 S12 S21 S⊤

11

• with S⊤

12 = −S12, and S⊤ 21 = −S21.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 11 / 54

Remark

The (S + S−1)-transform Ms − λLs, in general, is a nonlinear transform, e.g., the (S + S−1)-transform of M − λL ≡ A −H I

• − λ

I G A⊤

• arising in discrete-time optimal

controls produces a quadratic (S + S−1)-transform. The special form of the T-symplectic pencil (2) leads to a linear (S + S−1)-transform. This is an advantage by applying patel’s trick (Patel, LAA, 1993) to K − λN to find all eigenpair of P(λ).

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 12 / 54

If K and N are both real skew-Hamiltonian, Patel (LAA, 1993) developed an efficient structured algorithm for reducing K − µN to a block triangular condensed form QH(K − µN)Z = K11 K12 K⊤

11

• − µ

N11 N12 N⊤

11

• ,

where K11, N11 are upper Hessenberg and triangular, resp.. respectively, K12 and N12 are skew-symm., Y and Z are orthogonal matrices satisfying Z = J ⊤QJ .

Goal (Straightforward)

Extend Patel’s approach to complex ⊤-skew-Hamiltonian pencil (K, N) in (5).

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 13 / 54

Definition

A matrix U ∈ C2n×2n is unitary ⊤-symplectic, if UUH = I2n and UJ U⊤ = J . The set of unitary ⊤-symplectic matrices in C2n×2n is denoted by UTS2n.

Lemma

(i) Any matrix U ∈ UTS2n is of the form U = U1 ¯ U2 −U2 ¯ U1

• , where

U1, U2 ∈ Cn×n. (ii) For any matrix U ∈ UTS2n, if S is ⊤-skew-Hamiltonian, then UHSU is still ⊤-skew-Hamiltonian.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 14 / 54

A standard Givens rotation in C2n×2n is given by G(i, j, c, s) =       Ii−1 c ¯ s Ij−i−1 −s ¯ c I2n−j       , (6) where c¯ c + s¯ s = 1.

Three types of unitary ⊤-symplectic matrices

1 Gs1(i, c, s) := G(i, n + i, c, s) 2 Gs2(i, j, c, s) :=

G(i, j, c, s) ¯ G(i, j, c, s)

• 3 Hs(k, v) :=

H(k, v) ¯ H(k, v)

• , where

H(k, v) = Ik ⊕ (In−k − 2vvH

vHv) with v ∈ Cn−k.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 15 / 54

Two types of transformations in the structured algorithm

1 Using Givens ⊤-symplectic matrices Gs1(i, c, s):

K := GH

s1KGs1,

N := GH

s1NGs1.

(7)

2 Let U, V ∈ Cn×n be unitary representing the application of Givens

• rotations. Set

K := QH

0 KZ0,

N := QH

0 NZ0,

(8) where QH

0 =

U⊤ V ⊤

• ,

Z0 = V U

• .

(9)

Remark

Let (K, N) be ⊤-skew-Hamiltonian pencil. Then the resulting K and N in (7) and (8) are still ⊤-skew-Hamiltonian.

Sketch W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 16 / 54

Condensed form of the structured algorithm I

Let K and N be ⊤-skew-Hamiltonian. Then K − λN can be reduced to a block triangular structure, that is, K := QHKZ = K11 K12 K⊤

11

• ,

N := QHNZ = N11 N12 N⊤

11

• ,

(10) where K11 ∈ Cn×n is upper Hessenberg, N11 ∈ Cn×n is upper triangular, and Q, Z are unitary satisfying Z = J ⊤QJ .

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 17 / 54

Patel’s Algorithm: Reduction of K and N (n = 4)

            × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × ×             ⇓             × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × ×            

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 18 / 54

Structured Algorithm I (SA I)

Input : A palindromic QEP: P(λ) ≡ λ2A⊤

1 + λA0 + A1 with

A⊤

0 = A0;

Output : All eigenvalues and eigenvectors of P(λ). Step 1. Form the pencil K − λN in (5); Step 2. Reduce K − λN to block upper triangular forms in (10). Step 3. Use QZ alg. to compute eigenpairs {(µi, yi)}n

i=1 of

K11 − λN11; Step 4. Compute zi = J Z yi

• ≡

zi1 zi2

• , i = 1, 2, . . . , n;

Step 5. Compute eigenvalues νi and 1

νi of P(λ) by solving

ν2 − µiν + 1 = 0; Compute EVs xi1 ≡ zi1 + 1

νi zi2,

xi2 ≡ zi1 + νizi2 to νi, 1

νi resp., for i = 1, 2, . . . , n.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 19 / 54

Structured Algorithm II vs. URV-based structured Algorithm

A “good” linearization of the palindromic QEP [D. Mackey/N. Mackey/ Mehl/ Mehrmann SIMAX, 2006] is proposed by λZ⊤ + Z ≡ λ A⊤

1

A0 − A1 A⊤

1

A⊤

1

• +
• A1

A1 A0 − A⊤

1

A1

• (11)

We rewrite (11) into a new palindromic QEP by Q(µ) ≡ µ2Z⊤ + µ02n + Z, (12) where λ = µ2 and apply SA I algorithm to solve (12).

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 20 / 54

Structured Algorithm II (SA II)

As in (5) we form

• K =
• Z⊤ − Z

Z − Z⊤

• ,
• N =

−Z −Z⊤

• .

By (10) there are unitary Ua, Va ∈ C4n×4n with Ua = J ⊤

4nVaJ4n s.t.

UH

a

KVa = Ka

11

Ka

12

Ka

11 ⊤

• ,

UH

a

NVa = Na

11

Na

12

Na

11 ⊤

• ,

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 21 / 54

where Ka

11 =

      × × × × × × × × × ×       , N a

11 =

      × × × × × × × × ×       Ka

12; Na 12 =

      × × × × × × × × × × × ×       .

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 22 / 54

Denote P2n = [e1, e3, . . . , e2n−1|e2, e4, . . . , e2n] and let UH = (P⊤

2n ⊕ P⊤ 2n)UH a , V = Va(P2n ⊕ P2n).

We have UH ˆ KV =     R1 T1 R1 −T2 R⊤

2

R⊤

1

    , UH ˆ NV =     R3 −T3 R4 T ⊤

3

R⊤

3

R⊤

4

    Applying the periodic QZ to solve R1R−1

4 R2 − λR3 which gives

{(γi, yi)}n

i=1.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 23 / 54

Let µi = √γi, ηi = µiR−1

1 R3yi and ˆ

yi = [y⊤

i , η⊤ i ]⊤. Compute

ˆ zi = V ˆ yi

• ≡

ˆ zi1 ˆ zi2

• and solve νi, 1

νi for ν2 + (2 − γi)ν + 1 = 0. Compute

ˆ xi1 := (ˆ zi2 − 1 √νi ˆ zi1), ˆ xi2 := (ˆ zi2 − √νiˆ zi1). Then xi1 := ˆ xi1(1 : n) + ˆ xi1(n + 1 : 2n), xi2 := ˆ xi2(1 : n) + ˆ xi2(n + 1 : 2n) are EVs of P(λ) corresp. to νi and 1

νi , resp..

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 24 / 54

URV-based structued method [Schr¨

• der, 2007]:

U⊤ZV =     , V⊤(Z − Z⊤)V =     , U⊤(Z⊤ − Z)U =       . Let UH

0 =

   

• n
• n

I I     U⊤ V⊤

• ,

V0 = J ⊤

4nU0J4n,

where

n =

   1 ... 1   

n×n

.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 25 / 54

Define UH

b = (P2n ⊕ P2n)UH 0 ,

Vb = J ⊤

4nUbJ4n,

where UH

b ˆ

KVb =

• Kb

11

Kb

12

Kb

11 ⊤

• ,

UH

b

ˆ NVb =

• Nb

11

Nb

12

Nb

11 ⊤

• .

Theorem

If Ka

11 and Kb 11 are invertible, Na 11 and Nb 11 are nonsingular, then the SA II

algorithm is equivalent to the URV-based structured method.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 26 / 54

T-skew-Hamiltonian Arnoldi algorithm [Mehrmann/Watkins, SISC, 2001]

Let λ0 / ∈ σ(M, L), then µ0 ≡ λ0 + 1

λ0 /

∈ σ(K, N). Define the shift-invert transform K − µ N for K − µN with ˜ µ =

1 µ−µ0 by

• K

≡ −λ0N = −λ0LJ LT J T = λ0 A⊤

1

A1

• ,
• N

≡ −λ0(K − µ0N) = (M − λ0L) J

• MT − λ0LT

J T ≡ N1N2, where N1 = M − λ0L, N2 = J (M⊤ − λ0L⊤)J ⊤, Note that N1 and N2 are nonsingular and satisfy N ⊤

2 J = J N1.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 27 / 54

The generalized eigenvalue problem

• Kz =

µ Nz ≡ µN1N2z is then equivalent to the eigenvalue problem By ≡ N −1

1

• KN −1

2 y =

µy (13) with y = N2z. Using KJ = J K⊤ and N ⊤

2 J = J N1, we find that B

satisfies J B⊤ = BJ , and hence B is again ⊤-skew-Hamiltonian.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 28 / 54

• Two Theorems [Mehrmann/Watkins, SISC, 2001]

Theorem

Let B ∈ C2n×2n be ⊤-skew-Hamiltonian and let Kj ≡ Kj[B, u1] (1 ≤ j ≤ n) be a Krylov matrix with rank(Kj) = j. If Kj = UjRj is a QR-decomp., then BUj = UjHj + uj+1e⊤

j ,

where Hj ∈ Cj×j is unreduced upper Hessenberg Uj is ⊤-isotropic orthog. matrix and uj+1 ∈ C2n is a suitable vector s.t. UH

j

uj+1 = 0 and U⊤

j J

uj+1 = 0.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 29 / 54

Theorem

Let B ∈ C2n×2n be ⊤-skew-Hamiltonian. If rank(Kn[B, u1]) = n, then there is a unitary ⊤-symplectic matrix U with Ue1 = u1 such that UHBU = Hn Nn H⊤

n

• ,

where Hn is unreduced upper Hessenberg and Nn is ⊤-skew-symmetric.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 30 / 54

jth step of Arnoldi process

hj+1,juj+1 = Buj −

j

• i=1

hijui, (14) where hij = uH

i Buj, i = 1, . . . , j and hj+1,j > 0 is chosen so that

uj+12 = 1. In exact arithmetic, uj+1 computed by (14) should be orthog. against the vector J ¯ u1, . . . , J ¯ uj automatically. However, roundoff error causes the inner products of uH

j+1J ¯

ui, i = 1, . . . , j, to be tiny values. The jth step of Arnoldi process for zj should be modified.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 31 / 54

jth step of ⊤-isotropic Arnoldi process

hj+1,juj+1 = Buj −

j

• i=1

hijui −

j

• i=1

tijJ ¯ ui, (15) where hij = uH

i Buj, tij = −u⊤ i J Buj, i = 1, . . . , j

and hj+1,j > 0 is chosen so that uj+12 = 1.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 32 / 54

Algorithm (⊤SHIRA Algorithm for B = N −1

1

• KN −1

2 )

Generate a ℓth step of ⊤-isotropic Arnoldi factor.: BUℓ = UℓHℓ + hℓ+1,ℓuℓ+1e⊤

ℓ .

For k = 1, 2, . . . , until ℓ eigenpairs of B are convergent, Extend the ℓth step of ⊤-isotropic Arnoldi factor. to the (ℓ + p)th step of ⊤-isotropic Arnoldi factor.: BUℓ+p = Uℓ+pHℓ+p + hℓ+p+1,ℓ+puℓ+p+1e⊤

ℓ+p.

Reform a new ℓth step of ⊤-isotropic Arnoldi factor. by implicitly restarting processes. End Compute the eigenpairs ( µk, zk), k = 1, . . . , ℓ of Hℓ. Set µk := µ−1

k

+ λ0 + λ−1 and solve N2zk = Uℓ zk, k = 1, . . . , ℓ.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 33 / 54

Generalized ⊤-skew-Hamiltonian Arnoldi Algorithm

Goal

Develop a structured Arnoldi algorithm to compute the desired eigenpairs

• f the generalized eigenvalue problem
• Kz =

µ Nz, where K and N are ⊤-skew-Hamiltonian.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 34 / 54

Patel’s Condensed form:

Let K and N be ⊤-skew-Hamiltonian. Compute unitary U and V with U = J ⊤VJ such that VH KU = Hn Sn HT

n

• ,

VH NU = Rn Tn RT

n

• ,

where Hn is unreduced upper Hessenberg, Rn is nonsing. upper triang., Sn and Tn are ⊤-skew-symmetric.

⊤-bi-isotropic

X and Y are called ⊤-bi-isotropic if Y ⊤J X = 0.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 35 / 54

Theorem

Let B = N −1

1

• KN −1

2

be ⊤-skew-Hamiltonian defined in (13) with

• N = N1N2. Let Kj ≡ Kj[B, u1] be a Krylov matrix with rank(Kj) = j. If

N −1

2 Kj = ZjR2,j

and N1Kj = YjR1,j are QR-factor., then we have

• KZj = YjHj +

yj+1e⊤

j

and

• NZj = YjRj,

where Hj is unreduced upper Hessenberg, Rj is nonsing. upper triang., Yj and Zj are ⊤-bi-isotropic s.t. Y H

j

yj+1 = 0 and Z⊤

j J

yj+1 = 0, for a suitable yj+1 ∈ C2n.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 36 / 54

Theorem

Let B = N −1

1

• KN −1

2

be ⊤-skew-Hamiltonian in (12) with N = N1N2. If rank(Kn[B, u1]) = n, then ∃ unitary U and V satisfying U = J ⊤VJ and Ve1 = N1u1/N1u12 s.t. VH KU = Hn Sn HT

n

• ,

VHNU = Rn Tn RT

n

• ,

where Hn is unreduced upper Hessenberg, Rn is nonsing. upper triang. and Sn and Tn are ⊤-skew-symmetric.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 37 / 54

Give ⊤-bi-isotropic matrices Zj−1 and Yj satisfying

• KZj−1 = Yj−1Hj−1 +

yje⊤

j−1

and

• NZj = YjRj.

Comparing the jth column of both sides gives

• Nzj =

j−1

• i=1

rijyi + rjjyj. (16) For the (j − 1)th step of generalized ⊤-isotropic Arnoldi process generates

• NZj−1 = Yj−1Rj−1,

the vector zj in (16) can be reformulated as r−1

jj zj =

N −1yj −

j−1

• i=1
• rijzi,

(17) where [ r1j, . . . , rj−1,j]⊤ := −r−1

jj R−1 j−1[r1j, . . . , rj−1,j]⊤.

(18)

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 38 / 54

Using the orthogonality of Zj, the coefficient ˜ rij in (18) can be evaluated by ˜ rij = zH

j

N −1yj, i = 1, . . . , j − 1, and rjj in (17) is chosen so that zj2 = 1. Substituting [ r1j, . . . , rj−1,j]⊤ into (18), we obtain the coefficient vector [r1j, . . . , rj−1,j]⊤. In exact arithmetic, zj is orthogonal to J ¯ Yj automatically. Roundoff error causes z⊤

j J yi, i = 1, . . . , j, to be tiny values.

The jth step of generalized ⊤-isotropic Arnoldi process for zj should be modified.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 39 / 54

jth step of generalized ⊤-isotropic Arnoldi process for zj

r−1

jj zj =

N −1yj −

j−1

• i=1
• rijzi −

j

• i=1

sijJ ¯ yi, where

• rij

= zH

j

N −1yj, i = 1, . . . , j − 1, sij = y⊤

i J ⊤(

N −1yj −

j−1

• i=1
• rijzi),

i = 1, . . . , j.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 40 / 54

On the other hand, similar to (15), the jth step of generalized ⊤-isotropic Arnoldi process for yj+1 is given by hj+1,jyj+1 = Kzj −

j

• i=1

hijyi −

j

• i=1

tijJ ¯ zi, where hij = yH

i

Kzj, tij = z⊤

i J ⊤

Kzj, i = 1, . . . , j, and hj+1,j > 0 is chosen so that yj+12 = 1.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 41 / 54

Algorithm (G⊤SHIRA)

Generate a ℓth step of generalized ⊤-isotropic Arnoldi factor.:

• KZℓ = YℓHℓ + hℓ+1,ℓyℓ+1e⊤

ℓ ,

• NZℓ = YℓRℓ

For k = 1, 2, . . . , until ℓ eigenpairs of are convergent, Extend the ℓth step of generalized ⊤-isotropic Arnoldi

• factor. to the (ℓ + p)th step:
• KZℓ+p = Yℓ+pHℓ+p + hℓ+p+1,ℓ+pyℓ+p+1e⊤

ℓ+p,

• NZℓ+p = Yℓ+pRℓ+p.

Reform a new ℓth step of generalized ⊤-isotropic Arnoldi

• factor. by implicitly restarted processes.

End

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 42 / 54

Theorem

Let V ∈ Cn×r be unitary and A, B ∈ Cn×n. Then AV − BV C2 ≥ AV − BV P2, ∀C ∈ Cr×r, where P = (UHBV )−1(UHAV ) with BV = US a QR-factor. of BV.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 43 / 54

Numerical results (Small size)

Example for vibration of high-speed trains and rails: Using the finite element discretization with an uniform mesh, we

• btain a palindromic QEP:

(λ2A⊤

1 + λA0 + A1)x = 0,

where A0, A1 ∈ C5757×5757. After deflated zero and infinity eigenvalues, a small size palindromic quadratic eigenvalue problem: (λ2 A⊤

1 + λ

A0 + A1)x = 0, with A0, A1 ∈ C303×303 is given. To measure accuracy of an approximate eigenpair (λ, x), we use the relative residual RRes ≡ λ2 A⊤

1 x + λ

A0x + A1x2 (|λ2| A12 + |λ| A02 + A12)x2 .

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 44 / 54

10-20 10

• 10

10 10

10

10

20

10

• 18

10

• 16

10

• 14

10

• 12

10

• 10

10

• 8

| λ | R elative residuals of eigenpairs S +S

• 1

QZ UR V

The relative residuals of computed eigenpairs for the eigenvalues with absolute values in [10−20, 1020] and ω = 1000 For eigenvalues with small modulus, SA I performs better than SA II/URV and QZ algorithm. For eigenvalues near the unit circle or with large modulus, all three algorithm have similar accuracy.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 45 / 54

20 40 60 80 100 120 140 10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 index reciprocity of eigenvalues

r(λi) = min{|λiλj − 1|; j = i} SA I and SA II preserve the essential reciprocity property as expected For QZ algorithm,

◮ has only less than 18 pairs of computed eigenvalues near the unit circle

with reciprocity near zero.

◮ the smallest, the average and the maximal value of all reciprocities are

2.96 × 10−14, 0.223 and 0.952, respectively.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 46 / 54

QZ Algorithm SA_I Algorithm SA_II/URV Algorithm

Figure: The residuals of eigenvalues vs. ω.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 47 / 54

Numerical results (Large size)

Example for vibration of high-speed trains and rails: Using the finite element discretization with an uniform mesh, we have a QEP: (λ2A⊤

1 + λA0 + A1)x = 0,

where A0, A1 ∈ C5757×5757. To measure accuracy of an approximate eigenpair (λ, x), we use the residual RRes ≡ λ2A⊤

1 x + λA0x + A1x2

(|λ2|A12 + |λ|A02 + A12)x2 .

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 48 / 54

• Alg. 1
• Alg. 2

Solving linear system 2 2 Matrix-vector product for A1 4 4 Matrix-vector product for A0 2 2 Inner products 6j 8j − 2 Saxpy operations 6j + 6 8j + 4

Table: Computational cost of one step of ⊤-isotropic Arnoldi process in Alg. 1 and Alg. 2.

• Alg. 1 ≡ ⊤SHIRA
• Alg. 2 ≡ G⊤SHIRA

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 49 / 54

10-5 10 10

5

10

• 16

10

• 15

10

• 14

10

• 13

10

• 12

10

• 11

10

• 10

10

• 9

10

• 8

(a) | λ | R R es 10

• 5

10 10

5

10

• 17

10

• 16

10

• 15

10

• 14

10

• 13

10

• 12

10

• 11

(b) | λ | R R es

Figure: The RRes of the eigenpairs computed by Alg. 1 and 2. The notations “∆” and “×” are denoted the results computed by Alg. 1 and 2 , resp.. In (a) and (b), the frequency ω are equal to 50 and 2000, resp..

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 50 / 54

100 300 500 2 4 6 8 10 12 14 16 18 20 ω S tacked value (a.1) 500 3000 5000 2 4 6 8 10 12 14 16 18 20 ω S tacked value (a.2) 100 300 500 2 4 6 8 10 12 14 16 18 20 ω S tacked value (b.1 ) 500 3000 5000 2 4 6 8 10 12 14 16 18 20 ω S tacked value (b.2 )

1.0e-14 1.0e-15 1.0e-13 1.0e-12 1.0e-11 1.0e-10 1.0e-9

Figure: The stacked bars of ℓω,k for first 20 EWs close to one with different ω (k = 1, . . . , 7). The results in (a) and (b) are computed by Alg. 1 and 2 , resp..

ℓω,k = the number of RRes belonging to I1 = (0, 10−15], I2 = (10−15, 10−14], . . . , I7 = (10−10, 10−9). ω = 5, 10, . . . , 500, 550, 600, . . . , 5000.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 51 / 54

100 300 500 10-16 10

• 15

10

• 14

10

• 13

10

• 12

10

• 11

10

• 10

10 -9 10 -8 ω Average of R R es (a) 500 3000 5000 10

• 16

10

• 15

10

• 14

10

• 13

10

• 12

10

• 11

10

• 10

10

• 9

10

• 8

ω Average of R R es (b)

Figure: The average of RRes for the twelve eigenpairs computed by Alg. 1 and 2. The notations “∆” and “×” are denoted the results computed by Alg. 1 and 2 resp..

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 52 / 54

Conclusion

We transform a ⊤-palindromic QEP to a ⊤-skew-Hamiltonian pencil by (S + S−1)-transformation. We extend Patel’s algorithm to solve the ⊤-skew-Hamiltonian pencil by using unitary ⊤-symplectic matrices. The extension does not hold for the complex conjugate transpose case. Structured algorithm is equivalent to URV-based structured algorithm while applying to µ2Z⊤ + µ02n + Z. We develop a structured generalized ⊤-skew-Hamiltonian implicitly restarted Arnoldi algorithm to solve a large-scale ⊤-skew-Hamiltonian pencil. Numerical results show that the accuracy of desired eigenpairs satisfy SA I > SA II/URV > QZ. G⊤SHIRA > ⊤SHIRA.

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 53 / 54

Thank you for your attention!

W.W. Lin (Tsing Hua Univ.) Palindromic QEP Jan., 2008 54 / 54