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 • ( λ 2 M + λG + K ) v = 0 M T = M > 0 G T = − G K T = K ≤ 0 • quadratic eigenvalue problem • large, sparse (finite elements) • Find few eigenvalues nearest imaginary axis (and corresponding eigenvectors). 2

Problem: Optimal Control � BB T � � � 0 E A • − λ C T C − A T E T 0 (large and sparse) • Hamiltonian/skew-Hamiltonian � � 0 I • multiply by J = − I 0 � C T C − A T � � � E T 0 • − λ − BB T − E 0 − A • symmetric/skew-symmetric 3

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

Hamiltonian Matrices • H ∈ R 2 n × 2 n � � 0 I ∈ R 2 n × 2 n • J = 0 − I • H is Hamiltonian iff JH is symmetric. � � A K • H = , − A T N where K = K T and N = N T 5

Linearization • λ 2 Mv + λGv + Kv = 0 • w = λv , Mw = λMv � � � � � � � � 0 − K v G M v = 0 • − λ 0 0 − M w − M w • Ax − λNx = 0 • symmetric/skew-symmetric 6

Reduction to Hamiltonian Matrix • A − λN (symmetric/skew-symmetric) � � �� 0 I • N = R T JR J = − I 0 sometimes easy, always possible • Transform: A − λR T JR R − T AR − 1 − λJ J T R − T AR − 1 − λI • H = J T R − T AR − 1 is Hamiltonian. 7

Example � � G M • N = 0 − M − 1 � � � � � � 0 I 0 I 2 G I • N = R T JR = 1 − I 0 0 2 G M M • J T R − T AR − 1 = H � � � M − 1 � � � 0 0 I I 0 = − 1 − 1 0 − K 2 G I 2 G I 8

Sparse Representation of H • Krylov subspace methods • We just need to apply the operator. ( M = LL T ) M − 1 � � � � � � 0 0 I I 0 = H − 1 − 1 − K 0 2 G I 2 G I � � � ( − K ) − 1 � � � 0 0 I I 0 H − 1 = 1 1 0 M 2 G I 2 G 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 u 1 u n u 1 u n  ❅ ❅  ❅ ❅ ❅   ❅ ❅ ❅ ❅ ❅   ❅ ❅ ❅ ❅ ❅ U − 1 AU = ❅  ❅ ❅  ❅ ❅ ❅   ❅ ❅ ❅ ❅ ❅ • partial similarity transformation: � � � � � � ❅ ❅ ❅ = A · · · · · · u 1 u k u 1 u k ❅ ❅ ❅ ❅ ❅ ❅ + u k +1 β k e T k 11

• short Lanczos runs (breakdowns!!, no look-ahead) � � � � � � ❅ ❅ ❅ · · · = · · · A u 1 u k u 1 u k ❅ ❅ ❅ ❅ ❅ ❅ + u k +1 β k e T 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 − 1 AS = ˆ 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 ∈ R 2 n × 2 n � � �� 0 I • S T JS = J J = 0 − I � � • S = U V � U T � � � 0 I � � • = J U V V T − I 0 • U T JU = 0, V T JV = 0, U T JV = I • Subspaces are isotropic. 14

Isotropic Subspaces • y T Jx = 0 for all x , y ∈ U � � • U = · · · u 1 u k • U T JU = 0 • Structured methods build isotropic subspaces. 15

Skew-Hamiltonian Case Theorem: B skew Hamiltonian, x � = 0 ⇒ span { x, Bx, . . . , B j − 1 x } 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 j j � � ˜ q j +1 = Bq j − q i h ij − Jq i t ij i =1 i =1 • produces isotropic subspaces: Jq 1 , . . . , Jq k are orthogonal to q 1 , . . . , q k . • Theory t ij = 0 • Practice t ij = ǫ (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: � ( λ 2 M + λG + K ) v � ≤ 10 − 12 ( � v � = 1) 19

Further Experience • Fortran/C code • n ≈ 2 × 10 5 • Disadvantage: Eigenvectors cost extra. (eigenvectors of H 2 vs. H ) • We haven’t done skew-Hamiltonian Lanczos. 20

Hamiltonian Case • Bunse-Gerstner/Mehrmann 1986:   ❅ ❅ ❅ ❅ � � ❅ ❅ ❅ ❅ E T ❅ ❅ ❅ ❅ S − 1 HS = =     D − E ❅ ❅   ❅ ❅ ❅ ❅ • Further condensation: E = 0, D = diag {± 1 · · · ± 1 } . � � • S = U V � � � 0 T � � � = • H U V U V 0 D 21

Condensed Hamiltonian Lanczos Process   ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ � � � �   • H = U V U V   ❅   ❅ ❅ � � � � U = · · · V = · · · u 1 u 2 v 1 v 2 • Hu k = v k d k Hv k = u k − 1 b k − 1 + u k a k + u k +1 b k • u k +1 b k = Hv k − u k a k − u k − 1 b k − 1 v k +1 d k +1 = Hu k +1 • Coefficients are chosen so that � � S = is symplectic. U V • Collect coefficients. 22

Equivalence • H 2 is skew-Hamiltonian. • Condensed Hamiltonian Lanczos applied to H is theoretically equivalent to ordinary Lanczos applied to H 2 . • Hamiltonian algorithm costs half as many matrix-vector multiplies. 23

Isotropy � � 0 I • S T JS = J , J = 0 − I � U T � � � 0 I � � • = J U V V T − I 0 • U T JU = 0, (isotropic subspaces) V T JV = 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 = − JS T J 26

Symplectic Similarity • symplectic butterfly form: (Banse/Bunse-Gerstner 1994)   ❅ ❅ ❅ ❅ � � ❅ ❅ ❅ ❅ D 1 T 1 ❅ ❅ ❅ ❅ W − 1 SW =   =   D 2 T 2 ❅ ❅ ❅ ❅   ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ • Further condensation: D 1 = 0, D 2 = diag {± 1 · · · ± 1 } , T 1 = − D 2 , . . . � � • W = U V � � � 0 − D � � � = • S U V U V D DT 27

Condensed Symplectic Lanczos Process   ❅ ❅ ❅   ❅ � � � �   • S = U V U V   ❅ ❅ ❅   ❅ ❅ ❅ ❅   ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ • Su k = v k d k Sv k = − u k d k + v k − 1 ˜ a k + v k +1 ˜ b k − 1 + v k ˜ b k • v k +1 ˜ a k − v k − 1 ˜ = Sv k − v k ˜ b k − 1 + u k d k b k u k +1 d k +1 = S − 1 v k +1 • Coefficients are chosen so that � � is symplectic. U V • 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) � U T � � � 0 I � � • = J U V V T − I 0 • U T JU = 0, (isotropic subspaces) V T JV = 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