Invariant Subspace Computation in Scientific Computing: Gramian - PowerPoint PPT Presentation

Invariant Subspace Computation in Scientific Computing: Gramian Based Model Order Reduction RANMEP 2008 D.C. Sorensen ◮ M. Heinkenschloss, A.C. Antoulas , M. Shah and K. Sun ◮ Support: NSF and AFOSR NCTS Taiwan, Jan 2008

SVD Type Projection Methods for MOR Brief Intro to SVD in Model Reduction Balanced Reduction Gramian Based Model Reduction: Symmetry Preserving SVD: Symmetric Data Solving Large Lyapunov Equations: Approximate Balancing 1M vars now possible Balanced Reduction of Oseen Eqns: Extension to Descriptor System Couple Linear with Domain Decomposition: Nonlinear Domains Local Reduction ⇒ Many Interactions Neural Modeling:

Example: VLSI circuits 1960’s: IC 1971: Intel 4004 2001: Intel Pentium IV 10 µ details 0.18 µ details 2300 components 42M components 64KHz speed 2GHz speed 2km interconnect 7 layers D.C. Sorensen 3

Gramian Based Model Reduction Proper Orthogonal Decomposition (POD) Principal Component Analysis (PCA) x ( t ) = f ( x ( t ) , u ( t )) , ˙ y = g ( x ( t ) , u ( t )) The Gramian � ∞ x ( τ ) x ( τ ) T d τ P = o Eigenvectors of P P = VS 2 V T Orthogonal Basis x ( t ) = VSw ( t ) D.C. Sorensen 4

PCA or POD Reduced Basis Low Rank Approximation x ≈ V k ˆ x k ( t ) Galerkin condition – Global Basis ˙ x k = V T ˆ k f ( V k ˆ x k ( t ) , u ( t )) Global Approximation Error ( H 2 bound for LTI) � x − V k ˆ x k � 2 ≈ σ k +1 Snapshot Approximation to P m � P ≈ 1 x ( t j ) x ( t j ) T = XX T m j =1 Truncate SVD : X = VSU T ≈ V k S k U T k

Image Compression - Feature Detection original rank = 10 rank = 30 rank = 50

POD in CFD Extensive Literature Karhunen-Lo´ eve, L. Sirovich Burns, King Kunisch and Volkwein Gunzburger Many, many others Incorporating Observations – Balancing Lall, Marsden and Glavaski K. Willcox and J. Peraire D.C. Sorensen 7

Symmetry Preserving SVD M. Shah Backbone representation of HIV-1 protease (from M. Moll) Rotational symmetry should be present ◮ Get matrix representation of symmetry group generators Now possible for all interesting symmetries (Weyl 7 inf. series) ◮ Determine best axis (axes) of symmetry ◮ Compute best (low rank) symmetric approximation (SPSVD) Iteratively using ARPACK

The Symmetric SVD Approximation If WX 2 = X 1 + E where W = blockdiag ( I − 2 ww T ) � � X 1 � � ˆ �� � ˆ � 2 � � X 1 X 1 � � = USV T min − and � � ˆ ˆ X 2 X 2 X 2 W ˆ X 2 =ˆ X 1 F Solved by: � U 1 � √ 1 √ U = , S = 2 S 1 , V = V 1 . and U 2 = WU 1 , U 2 2 with 1 = 1 U 1 S 1 V T 2 ( X 1 + WX 2 ) D.C. Sorensen 9

SPSVD - 2Fold Dihedral Symmetry Figure: Different approximations of 1IDS (6K atoms, rank 6 approx.)

LTI Model Reduction by Projection ˙ x = Ax + Bu = y Cx Approximate x ∈ S V = Range ( V ) k -diml. subspace i.e. Put x = V ˆ x , and then force W T [ V ˙ ˆ − ( AV ˆ x + Bu )] = 0 x ˆ y = CV ˆ x If W T V = I k , then the k dimensional reduced model is ˙ ˆ x + ˆ ˆ x = A ˆ Bu ˆ ˆ = C ˆ y x where ˆ A = W T AV , ˆ B = W T B , ˆ C = CV .

POD for LTI systems H ( t ) = C ( t I − A ) − 1 B , Impulse Response: t ≥ 0 x ( t ) = e A t B Input to State Map: Controllability Gramian: � ∞ � ∞ e A τ BB T e A T τ d τ x ( τ ) x ( τ ) T d τ = P = o o y ( t ) = C e A t x (0) State to Output Map: Observability Gramian: � ∞ e A T τ C T C e A τ d τ Q = o D.C. Sorensen 12

Balanced Reduction (Moore 81) Lyapunov Equations for system Gramians A P + P A T + BB T = 0 A T Q + Q A + C T C = 0 With P = Q = S : Want Gramians Diagonal and Equal States Difficult to Reach are also Difficult to Observe Reduced Model A k = W T k AV k , B k = W T k B , C k = C k V k ◮ P V k = W k S k Q W k = V k S k ◮ Reduced Model Gramians P k = S k and Q k = S k . D.C. Sorensen 13

Hankel Norm Error estimate (Glover 84) Why Balanced Truncation? � ◮ Hankel singular values = λ ( PQ ) ◮ Model reduction H ∞ error (Glover) � y − ˆ y � 2 ≤ 2 × ( sum neglected singular values ) � u � 2 ◮ Extends to MIMO ◮ Preserves Stability Key Challenge ◮ Approximately solve large scale Lyapunov Equations in Low Rank Factored Form D.C. Sorensen 14

CD Player Frequency Response Freq − Response CD − Player : τ = 0.001, n = 120 , k = 12 2 10 Original Reduced 0 10 − 2 10 |G(j ω )| − 4 10 − 6 10 − 8 10 − 10 10 0 1 2 3 4 5 6 7 10 10 10 10 10 10 10 10 Frequency ω

CD Player Frequency Response Freq − Response CD − Player : τ = 1e − 005, n = 120 , k = 37 2 10 Original Reduced 0 10 − 2 10 |G(j ω )| − 4 10 − 6 10 − 8 10 − 10 10 0 1 2 3 4 5 6 7 10 10 10 10 10 10 10 10 Frequency ω

� CD Player - Hankel Singular Values λ ( PQ ) Hankel Singular Values 2 10 0 10 −2 10 −4 10 −6 10 −8 10 −10 10 −12 10 −14 10 0 20 40 60 80 100 120

Approximate Balancing A P + P A T + BB T = 0 A T Q + Q A + C T C = 0 • Sparse Case: Iteratively Solve in Low Rank Factored Form, P ≈ U k U T Q ≈ L k L T k , k [ X , S , Y ] = svd ( U T k L k ) W k = LY k S − 1 / 2 and V k = UX k S − 1 / 2 . k k Now: P W k ≈ V k S k and Q V k ≈ W k S k D.C. Sorensen 18

Low Rank Smith = ADI Convert to Stein Equation: A P + P A T + BB T = 0 P = A µ P A T µ + B µ B T ⇐ ⇒ µ , where � A µ = ( A − µ I )( A + µ I ) − 1 , 2 | µ | ( A + µ I ) − 1 B . B µ = Solution: ∞ � µ ) T = LL T , A j µ B µ B T µ ( A j P = j =0 where L = [ B µ , A µ B µ , A 2 µ B µ , . . . ] Factored Form D.C. Sorensen 19

Multi-Shift (Modified) Low Rank Smith LR - Smith: Update Factored Form P m = L m L T m : ( Penzl , Li , White ) L m +1 = [ A µ L m , B µ ] [ A m +1 = B µ , L m ] µ Multi-Shift LR - Smith: (Gugercin, Antoulas, and S.) Update and Truncate SVD Re-Order and Aggregate Shift Applications Much Faster and Far Less Storage B ← A µ B ; [ V , S , Q ] = svd ([ A µ B , L m ]); L m +1 ← V k S k ; ( σ k +1 < tol · σ 1 ) D.C. Sorensen 20

Approximate Power Method (Hodel) T U + BB T U = 0 A P U + P A T ) A T A T U + BB T U + P ( I − UU T U = 0 A P U + P UU Thus T + BB T U ≈ 0 where H = U T AU A P U + P UH Solving T + BB T U = 0 AZ + ZH gives approximation to Z ≈ P U Iterate ⇒ Approximate Power Method Z j → US with P U = US (also see Vasilyev and White 05)

T ) A Parameter Free Synthesis ( P ≈ US 2 U Step 1: Solve the reduced order Lyapunov equation P H T + ˆ B T = 0 . H ˆ P + ˆ B ˆ Solve ˆ with H = U k T AU k , B = U k T B . Step 2: (APM step) Solve a projected Sylvester equation AZ + ZH T + B ˆ B T = 0 , Step 3: Modify B B ← ( I − Z ˆ P − 1 U T ) B . Update Step 4: (ADI step) Update factorization and basis U k Z ← Z ˆ P − 1 / 2 . Re-scale Update (and truncate) [ U , S ] ← svd [ US , Z ] . U k ← U (: , 1 : k ), basis for dominant subspace. D.C. Sorensen 22

Automatic Shift Selection - Placement? SUPG discretization advection-diffusion operator on square grid k = 32, m = 59 , n = 32*32, Thanks Embree Comparison Eigenvalues A vs H 0.03 0.02 0.01 0 −0.01 −0.02 −0.03 −0.07 −0.06 −0.05 −0.04 −0.03 −0.02 −0.01 0 0.01 eigenvalues A eigenvalues small H eigenvalues full H

Convergence History , Supg, n = 32, N = 1024 Laptop �P + −P� � ˆ Iter � B j � B j � �P + � 1 2.7e-1 1.6e+0 4.7e+0 2 7.2e-2 1.6e-1 1.5e+0 3 6.6e-3 1.1e-2 1.3e-1 4 3.7e-4 3.5e-7 1.3e-2 5 9.3e-9 1.4e-11 3.4e-7 P f is rank k = 59 Comptime( P f ) = 16.2 secs �P−P f � = 8 . 8 e − 9 Comptime( P ) = 810 secs �P� D.C. Sorensen 24

Convergence History , Supg, n = 800, N = 640,000 CaamPC �P + −P� � ˆ Iter � B j � B j � �P + � 6 1.3e-01 2.5e+00 2.4e+00 7 7.5e-02 1.1e+00 1.2e+00 8 3.5e-02 6.74e-01 5.0e-01 9 2.0e-02 1.2e-02 6.7e-01 10 2.0e-04 7.1e-07 1.2e-02 11 1.0e-08 2.3e-11 6.4e-07 P f is rank k = 120 Comptime( P f ) = 157 mins = 2.6 hrs D.C. Sorensen 25

Oseen Equations: A Descriptor System d dt v ( t ) = A 11 v ( t ) + A 12 p ( t ) + B 1 g ( t ) , E 11 A T 0 = 12 v ( t ) , v (0) = v 0 , y ( t ) = C 1 v ( t ) + C 2 p ( t ) + Dg ( t ) . E d dt v ( t ) = A v ( t ) + B u ( t ) , y ( t ) = C v ( t ) + D u ( t ) Note E is singular index 2 See Stykel (LAA 06) – general theory and approach

Notation Put � � � � E = ΠE 11 Π T , A = ΠA 11 Π T , C = CΠ T . B = ΠB 1 , With this notation, A T + � A P � � E + � E P � B � B T = 0 , A T Q � C T � � E + � E Q � A + � C = 0 . where 12 E − 1 11 A 12 ) − 1 A T 12 E − 1 I − A 12 ( A T Π = 11 Θ l Θ T = r Π T Projector onto Null( A T 12 )

ADI Derivation Step 1 Begin with � � T �� � B T � A T + � � � µ � � µ � P � B � AP E + ¯ A = − E − ¯ A . and derive � � � � T � � � � T � E + µ � � µ � � µ � E − µ � � − 2 Re( µ ) � B � B T . A P E + ¯ A = E − ¯ A P A � � E + µ � � Problem: A is Singular

Key Pseudo-Inverse Lemma Suppose Θ T r E 11 Θ r + µ Θ T r A 11 Θ r is invertible. Then the matrix � � I � � − 1 E + µ � � Θ T r E 11 Θ r + µ Θ T Θ T ≡ Θ r A r A 11 Θ r r satisfies � � I � � � E + µ � E + µ � � = Π T A A and � � � � I � E + µ � E + µ � � A A = Π .