Hamiltonian and Symplectic Lanczos Processes David S. Watkins - - PDF document

Hamiltonian and Symplectic Lanczos Processes David S. Watkins

Department of Mathematics Washington State University Pullman, WA 99164-3113 U.S.A. PIMS Workshop, Vancouver BC August 2003 Collaborators: V. Mehrmann, T. Apel,

• P. Benner, H. Faßbender, . . .

1

Problem: Linear Elasticity

• (λ2M + λG + K)v = 0

MT = M > 0 GT = −G KT = K ≤ 0

• quadratic eigenvalue problem
• large, sparse (finite elements)
• Find few eigenvalues nearest imaginary axis

(and corresponding eigenvectors).

2

Problem: Optimal Control

• A

BBT CTC −AT

• − λ
• E

ET

• (large and sparse)
• Hamiltonian/skew-Hamiltonian
• multiply by J =
• I

−I

• CTC

−AT −A −BBT

• − λ
• ET

−E

• symmetric/skew-symmetric

3

Hamiltonian Structure

• Our matrices are real.
• λ, λ, −λ, −λ occur together.
• seen also in Hamiltonian matrices

4

Hamiltonian Matrices

• H ∈ R2n×2n
• J =
• I

−I

• ∈ R2n×2n
• H is Hamiltonian iff JH is symmetric.
• H =
• A

K N −AT

• ,

where K = KT and N = NT

5

Linearization

• λ2Mv + λGv + Kv = 0
• w = λv,

Mw = λMv

• −K

−M v w

• −λ
• G

M −M v w

• = 0
• Ax − λNx = 0
• symmetric/skew-symmetric

6

Reduction to Hamiltonian Matrix

• A − λN

(symmetric/skew-symmetric)

• N = RTJR
• J =
• I

−I

• sometimes easy, always possible
• Transform:

A − λRTJR R−TAR−1 − λJ JTR−TAR−1 − λI

• H = JTR−TAR−1 is Hamiltonian.

7

Example

• N =
• G

M −M

• N = RTJR =
• I

−1

2G

M I −I I

1 2G

M

• H

= JTR−TAR−1 =

• I

−1

2G

I M−1 −K I −1

2G

I

• 8

Sparse Representation of H

• Krylov subspace methods
• We just need to apply the operator.

(M = LLT) H =

• I

−1

2G

I M−1 −K I −1

2G

I

• H−1

=

• I

1 2G

I (−K)−1 M I

1 2G

I

• 9

Exploitable Structures

• Hamiltonian

H−1 H−1(H − τI)−1(H + τI)−1

• skew-Hamiltonian

H−2 (H − τI)−1(H + τI)−1

• symplectic

(H − τI)−1(H + τI) τ = target shift Note: (H −τI)−1 has none of these structures.

10

Unsymmetric Lanczos Process

• Standard unsymmetric Lanczos effects a

(partial) similarity transformation A

• u1

· · · un

• =
• u1

· · · un

• 

 

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅

  

U−1AU =

  

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅

  

• partial similarity transformation:

A

• u1

· · · uk

• =
• u1

· · · uk

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅

• +uk+1βkeT

k

11

• short Lanczos runs

(breakdowns!!, no look-ahead) A

• u1

· · · uk

• =
• u1

· · · uk

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅

• +uk+1βkeT

k

• Get eigenvalues of
• ❅

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅

• Restart (implicitly)

IRA (Sorensen 1991), ARPACK Restart Lanczos with HR (Grimme/Sorensen/Van Dooren 1996)

12

Structured Lanczos Methods

• Similarity transformation: S−1AS = ˆ

A

• S symplectic ⇒ structure preserved
• symplectic (Lie group)
• Hamiltonian (Lie algebra)
• skew-Hamiltonian (Jordan algebra)
• Conclusion: A “Lanczos” process that builds

a symplectic similarity transformation will preserve structure. Vectors produced should be columns of a symplectic matrix.

13

Symplectic Matrices

• S ∈ R2n×2n
• STJS = J
• J =
• I

−I

• S =
• U

V

• UT

V T

• J
• U

V

• =
• I

−I

• UTJU = 0, V TJV = 0, UTJV = I
• Subspaces are isotropic.

14

Isotropic Subspaces

• yTJx = 0 for all x, y ∈ U
• U =
• u1

· · · uk

• UTJU = 0
• Structured methods build isotropic subspaces.

15

Skew-Hamiltonian Case

Theorem: B skew Hamiltonian, x = 0 ⇒ span{x, Bx, . . . , Bj−1x} is isotropic.

• Conclusion: Krylov subspace methods

preserve skew-Hamiltonian structure automatically.

• Examples: Arnoldi, unsymmetric Lanczos
• exact vs. floating-point arithmetic

16

Skew-Hamiltonian Arnoldi Process

• Isotropic Arnoldi process

˜ qj+1 = Bqj −

j

• i=1

qihij −

j

• i=1

Jqitij

• produces isotropic subspaces:

Jq1, . . . , Jqk are orthogonal to q1, . . . , qk.

• Theory tij = 0
• Practice tij = ǫ

(roundoff)

• Enforcement of isotropy is crucial.
• Consequence: get each eigenvalue only once.

17

Example

• Method: Implicitly Restarted Arnoldi

(effective combination of Arnoldi and subspace iteration)

• Toy problem (n = 64); asking for 8

eigenvalues (right half-plane).

• Target τ = i (not particularly good)
• After 12 Arnoldi steps (no restart) . . .

18

• After one restart (12 more Arnoldi steps)
• Errors: 10−14, 10−7, 10−6, 10−2
• After 7 iterations (restarts) algorithm stops

with 8 eigenvalues correct to ten decimal places.

• Residuals: (λ2M + λG + K)v ≤ 10−12

(v = 1)

19

Further Experience

• Fortran/C code
• n ≈ 2 × 105
• Disadvantage: Eigenvectors cost extra.

(eigenvectors of H2 vs. H)

• We haven’t done

skew-Hamiltonian Lanczos.

20

Hamiltonian Case

• Bunse-Gerstner/Mehrmann 1986:

S−1HS =

• E

T D −E

• =

   

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅

   

• Further condensation: E = 0,

D = diag{±1 · · · ± 1}.

• S =
• U

V

• H
• U

V

• =
• U

V T D

• 21

Condensed Hamiltonian Lanczos Process

• H
• U

V

• =
• U

V

• 

  

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅

   

U =

• u1

u2 · · ·

• V =
• v1

v2 · · ·

• Huk = vkdk

Hvk = uk−1bk−1 + ukak + uk+1bk

• uk+1bk

= Hvk − ukak − uk−1bk−1 vk+1dk+1 = Huk+1

• Coefficients are chosen so that

S =

• U

V

• is symplectic.
• Collect coefficients.

22

Equivalence

• H2 is skew-Hamiltonian.
• Condensed Hamiltonian Lanczos applied to

H is theoretically equivalent to ordinary Lanczos applied to H2.

• Hamiltonian algorithm costs half as many

matrix-vector multiplies.

23

Isotropy

• STJS = J,

J =

• I

−I

• UT

V T

• J
• U

V

• =
• I

−I

• UTJU = 0,

(isotropic subspaces) V TJV = 0,

• In (floating-point) practice, isotropy must

be enforced by J-reorthogonalization.

• All vectors must be retained.
• short Lanczos runs, restarts

24

Implicitly Restarted Hamiltonian Lanczos Process

• use SR, not QR

(Benner/Fassbender 1997)

• In condensed case, SR = HR

(Benner/Fassbender/W 1998)

• Use of HR yields significant simplification.

25

Symplectic Case Structure

• Eigenvalues of S appear in quartets

µ, µ−1, µ, µ−1.

• Symplectic Lanczos process must extract

these simultaneously.

• This is accomplished

by using both S and S−1.

• S−1 = −JSTJ

26

Symplectic Similarity

• symplectic butterfly form:

(Banse/Bunse-Gerstner 1994) W −1SW =

• D1

T1 D2 T2

• =

   

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅

   

• Further condensation: D1 = 0,

D2 = diag{±1 · · · ± 1}, T1 = −D2, . . .

• W =
• U

V

• S
• U

V

• =
• U

V −D D DT

• 27

Condensed Symplectic Lanczos Process

• S
• U

V

• =
• U

V

• 

    

❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅

     

• Suk = vkdk

Svk = −ukdk + vk−1˜ bk−1 + vk˜ ak + vk+1˜ bk

• vk+1˜

bk = Svk − vk˜ ak − vk−1˜ bk−1 + ukdk uk+1dk+1 = S−1vk+1

• Coefficients are chosen so that
• U

V

• is symplectic.
• Collect coefficients.

28

Equivalence

• S + S−1 is skew-Hamiltonian.
• Condensed symplectic Lanczos applied to

S is theoretically equivalent to ordinary Lanc- zos applied to S + S−1.

• Symplectic algorithm costs half as many

matrix-vector multiplies.

29

Isotropy (rerun)

• UT

V T

• J
• U

V

• =
• I

−I

• UTJU = 0,

(isotropic subspaces) V TJV = 0,

• In (floating-point) practice, isotropy must

be enforced by J-reorthogonalization.

• All vectors must be retained.
• short Lanczos runs, restarts

30

Implicitly Restarted Symplectic Lanczos Process

• use symplectic SR, not QR
• In condensed case, SR = HR

(Benner/Fassbender/W 1998)

• Use of HR yields significant simplification.

31

Remarks on Stability

• Both Hamiltonian and symplectic Lanczos

processes are potentially unstable.

• Breakdowns can occur.
• Are the answers worth anything?
• right and left eigenvectors
• residuals
• condition numbers for eigenvalues
• Don’t skip these tests.

32

Example

• λ2Mv + λGv + Kv = 0
• n = 3423

H =

• I

−1

2G

I M−1 −K I −1

2G

I

• H−1

=

• I

1 2G

I (−K)−1 M I

1 2G

I

• 33

Compare various approaches:

• Hamiltonian(1)

H−1

• Hamiltonian(3)

H−1(H−τI)−1(H+τI)−1

• symplectic

(H − τI)−1(H + τI)

• unstructured

(H − τI)−1 + ordinary Lanczos with implicit restarts Get 6 smallest eigenvalues in right half-plane. Tolerance = 10−8 Take 20 steps and restart with 10.

34

No-Clue Case (τ = 0)

Method Solves Eigvals Max. Found Resid. Hamiltonian(1) 78 9 2 × 10−10 Unstructured 158 7 + 7 5 × 10−7 Unstructured code must find everything twice.

−2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5 −5 −4 −3 −2 −1 1 2 3 4 5 x 10

−3

35

Conservative Shift (τ = 0.5)

Method Solves Eigvals Max. Found Resid. Hamiltonian(1) 78 9 2 × 10−10 Unstructured 138 7 + 2 3 × 10−5 Hamiltonian(3) 174 11 3 × 10−13 Symplectic 156 11 2 × 10−8

−2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5 −5 −4 −3 −2 −1 1 2 3 4 5 x 10

−3

36

Aggressive Shift (τ = 1.5)

Method Solves Eigvals Max. Found Resid. Hamiltonian(1) 78 9 2 × 10−10 Unstructured 96 9 1 × 10−7 Hamiltonian(3) 120 9 2 × 10−12 Symplectic 156 11 2 × 10−11

−2.5 −2 −1.5 −1 −0.5 0.5 1 1.5 2 2.5 −5 −4 −3 −2 −1 1 2 3 4 5 x 10

−3

37

The Last Slide

• We have developed structure-preserving

implicitly-restarted Lanczos methods for Hamiltonian and symplectic eigenvalue problems.

• The structure-preserving methods are more

accurate than a comparable non-structured method.

• By exploiting structure we can solve our

problems more economically.

38