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 Gramian Based Model Reduction: Balanced 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 Domain Decomposition: Couple Linear with Nonlinear Domains Neural Modeling: Local Reduction ⇒ Many Interactions

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 P = ∞

• x(τ)x(τ)Tdτ

Eigenvectors of P P = VS2VT Orthogonal Basis x(t) = VSw(t)

D.C. Sorensen 4

PCA or POD Reduced Basis

Low Rank Approximation x ≈ Vkˆ xk(t) Galerkin condition – Global Basis ˙ ˆ xk = VT

k f(Vkˆ

xk(t), u(t)) Global Approximation Error (H2 bound for LTI) x − Vkˆ xk2 ≈ σk+1 Snapshot Approximation to P P ≈ 1 m

m

• j=1

x(tj)x(tj)T = XXT Truncate SVD : X = VSUT ≈ VkSkUT

k

Image Compression - Feature Detection

• riginal

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 WX2 = X1 + E where W = blockdiag(I − 2wwT) min

Wˆ X2=ˆ X1

• X1

X2

• −

ˆ X1 ˆ X2

• 2

F

and ˆ X1 ˆ X2

• = USVT

Solved by: U = 1 √ 2 U1 U2

• ,

S = √ 2S1, V = V1. and U2 = WU1, with

U1S1VT

1 = 1 2 (X1 + WX2)

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 ∈ SV = Range(V) k-diml. subspace i.e. Put x = Vˆ x, and then force WT[V˙ ˆ x − (AVˆ x + Bu)] = 0 ˆ y = CVˆ x If WTV = Ik, then the k dimensional reduced model is ˙ ˆ x = ˆ Aˆ x + ˆ Bu ˆ y = ˆ Cˆ x where ˆ

A = WTAV, ˆ B = WTB, ˆ C = CV.

POD for LTI systems

Impulse Response: H(t) = C(tI − A)−1B, t ≥ 0 Input to State Map: x(t) = eAtB Controllability Gramian: P = ∞

• x(τ)x(τ)Tdτ =

∞

• eAτBBTeAT τdτ

State to Output Map: y(t) = CeAtx(0) Observability Gramian: Q = ∞

• eAT τCTCeAτdτ

D.C. Sorensen 12

Balanced Reduction (Moore 81)

Lyapunov Equations for system Gramians AP + PAT + BBT = 0 ATQ + QA + CTC = 0 With P = Q = S : Want Gramians Diagonal and Equal States Difficult to Reach are also Difficult to Observe Reduced Model Ak = WT

k AVk , Bk = WT k B , Ck = CkVk

◮ PVk = WkSk

QWk = VkSk

◮ Reduced Model Gramians Pk = Sk and Qk = Sk. D.C. Sorensen 13

Hankel Norm Error estimate (Glover 84)

Why Balanced Truncation?

◮ Hankel singular values =

• λ(PQ)

◮ Model reduction H∞ error (Glover)

y − ˆ y2 ≤ 2 × (sum neglected singular values)u2

◮ 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

10 10

1

10

2

10

3

10

4

10

5

10

6

10

7

10

−10

10

−8

10

−6

10

−4

10

−2

10 10

2

Freq−Response CD−Player : τ = 0.001, n = 120 , k = 12

|G(jω)| Frequency ω

Original Reduced

CD Player Frequency Response

10 10

1

10

2

10

3

10

4

10

5

10

6

10

7

10

−10

10

−8

10

−6

10

−4

10

−2

10 10

2

Freq−Response CD−Player : τ = 1e−005, n = 120 , k = 37

|G(jω)| Frequency ω

Original Reduced

CD Player - Hankel Singular Values

• λ(PQ)

20 40 60 80 100 120 10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 10

2

Hankel Singular Values

Approximate Balancing

AP + PAT + BBT = 0 ATQ + QA + CTC = 0

• Sparse Case: Iteratively Solve in Low Rank Factored Form,

P ≈ UkUT

k ,

Q ≈ LkLT

k

[X, S, Y] = svd(UT

k Lk)

Wk = LYkS−1/2

k

and Vk = UXkS−1/2

k

. Now: PWk ≈ VkSk and QVk ≈ WkSk

D.C. Sorensen 18

Low Rank Smith = ADI

Convert to Stein Equation: AP + PAT + BBT = 0 ⇐ ⇒ P = AµPAT

µ + BµBT µ ,

where Aµ = (A − µI)(A + µI)−1, Bµ =

• 2|µ|(A + µI)−1B.

Solution: P =

∞

• j=0

Aj

µBµBT µ (Aj µ)T = LLT,

where L = [Bµ, AµBµ, A2

µBµ, . . . ]

Factored Form

D.C. Sorensen 19

Multi-Shift (Modified) Low Rank Smith

LR - Smith: Update Factored Form Pm = LmLT

m:

(Penzl, Li, White) Lm+1 = [AµLm, Bµ] = [Am+1

µ

Bµ, Lm] 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, Lm]); Lm+1 ← VkSk; (σk+1 < tol · σ1)

D.C. Sorensen 20

Approximate Power Method (Hodel)

APU + PA

T U + BB T U = 0

APU + PUU

T A T U + BB T U + P(I − UU T )A T U = 0

Thus APU + PUH

T + BB T U ≈ 0 where H = U T AU

Solving AZ + ZH

T + BB T U = 0

gives approximation to Z ≈ PU Iterate ⇒ Approximate Power Method Zj → US with PU = US (also see Vasilyev and White 05)

A Parameter Free Synthesis (P ≈ US2U

T)

Step 1: Solve the reduced order Lyapunov equation Solve H ˆ P + ˆ PHT + ˆ Bˆ BT = 0. with

H = Uk TAUk, ˆ B = Uk TB.

Step 2: (APM step) Solve a projected Sylvester equation AZ + ZHT + Bˆ BT = 0, Step 3: Modify B Update B ← (I − Z ˆ P−1UT)B. Step 4: (ADI step) Update factorization and basis Uk Re-scale Z ← Z ˆ P−1/2. Update (and truncate) [U, S] ← svd[US, Z]. Uk ← 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

−0.07 −0.06 −0.05 −0.04 −0.03 −0.02 −0.01 0.01 −0.03 −0.02 −0.01 0.01 0.02 0.03 Comparison Eigenvalues A vs H

eigenvalues A eigenvalues small H eigenvalues full H

Convergence History , Supg, n = 32, N = 1024

Laptop Iter

P+−P P+

Bj ˆ Bj 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 Pf is rank k = 59 Comptime(Pf ) = 16.2 secs

P−Pf P

= 8.8e − 9 Comptime(P) = 810 secs

D.C. Sorensen 24

Convergence History , Supg, n = 800, N = 640,000

CaamPC Iter

P+−P P+

Bj ˆ Bj 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 Pf is rank k = 120 Comptime(Pf ) = 157 mins = 2.6 hrs

D.C. Sorensen 25

Oseen Equations: A Descriptor System

E11 d dt v(t) = A11v(t) + A12p(t) + B1g(t), = AT

12v(t),

v(0) = v0, y(t) = C1v(t) + C2p(t) + Dg(t). E d dt v(t) = Av(t) + Bu(t), y(t) = Cv(t) + Du(t) Note E is singular index 2 See Stykel (LAA 06) – general theory and approach

Notation

Put

• E = ΠE11ΠT,
• A = ΠA11ΠT,
• B = ΠB1,
• C = CΠT.

With this notation,

• A P

E + E P AT + B BT = 0,

• AT Q

E + E Q A + CT C = 0. where Π = I − A12(AT

12E−1 11 A12)−1AT 12E−1 11

= ΘlΘT

r

ΠT Projector onto Null(AT

12)

ADI Derivation Step 1

Begin with

• AP
• E + ¯

µ A T = −

• E − ¯

µ A

• P

AT + B BT . and derive

• E + µ

A

• P
• E + ¯

µ A T =

• E − ¯

µ A

• P
• E − µ

A T −2 Re(µ) B BT. Problem:

• E + µ

A

• is Singular

Key Pseudo-Inverse Lemma

Suppose ΘT

r E11Θr + µ ΘT r A11Θr is invertible.

Then the matrix

• E + µ

A I ≡ Θr

• ΘT

r E11Θr + µ ΘT r A11Θr

−1 ΘT

r

satisfies

• E + µ

A I

• E + µ

A

• = ΠT

and

• E + µ

A

• E + µ

A I = Π.

Projected Stein Equation

P = AµP A∗

µ − 2 Re( µ)

Bµ B∗

µ.

where

• Aµ ≡
• E + µ

A I

• E − ¯

µ A

• ,

and Bµ ≡

• E + µ

A I B Solution: P = −2 Re(µ)

∞

• j=0
• Aj

µ

Bµ B∗

µ

• A∗

µ

j . Convergent for stable pencil with Real(µ) < 0

Key Implementation Lemma

If M = ΠTM, then the computation Z =

• E + µ

A I

• E − ¯

µ A

• M

may be accomplished with the following steps.

• 1. Put

F = (E11 − ¯ µ A11) M.

• 2. Solve

E11 + µ A11 A12 AT

12

Z Λ

• =

F

• .

Note that Z satisfies Z = ΠTZ Similar result holds for computing Bµ

Algorithm:Single Shift ADI

• 1. Solve

E11 + µ A11 A12 AT

12

Z Λ

• =

B

• ;
• 2. U = Z;
• 3. while ( ’not converged’)

3.1 Z ← (E11 − ¯ µ A11) Z; 3.2 Solve (in place) E11 + µ A11 A12 AT

12

Z Λ

• =

Z

• ;

3.3 U ← [U, Z] ;

end

• 4. U ←
• 2| Re(µ)| U.

Derivation Multi-Shift ADI

Easy to see (P − Pk) = Ak

µP

• A∗

µ

k . Hence

• A(P − Pk)

E + E(P − Pk) A∗ =

• A

Ak

µP

• Ak

µ

∗ E + E Ak

µP

• Ak

µ

∗ A∗ =

• Ak

µ

• AP

E + EP A∗

• Ak

µ

∗ . To get

• A (P − Pk)

E + E (P − Pk) A∗ = − Ak

µ

• B

B∗

• Ak

µ

∗ . Where

• Aµ ≡
• E − ¯

µ A

• E + µ

A I .

Algorithm:Multi-Shift ADI

• 1. U = [ ];
• 2. while ( ’not converged’)

for i = 1:m, 2.1 Solve E11 + µi A11 A12 AT

12

Z Λ

• =

B

• ;

2.2 U0 = Z; 2.3 for j = 1:k-1,

2.3.1 Z ← (E11 − ¯ µi A11) Z; 2.3.2 Solve (in place) ✒ E11 + µi A11 A12 AT

12

✓ ✒ Z Λ ✓ = ✒ Z ✓ ; 2.3.3 U0 ← [U0, Z] ;

end 2.4 U ←

• U,
• 2| Re(µi)| U0
• ;

% Update and truncate SVD(U); 2.5 B ← (E11 − ¯ µi A11) Z. end

end

Model Problem: Oseen Equations

∂ ∂t v(x, t) + (a(x) · ∇)v(x, t) − ν∆v(x, t) + ∇p(x, t) = χΩg (x)gΩ(x, t) in Ω × (0, T), ∇ · v(x, t) = in Ω × (0, T), (−p(x, t)I + ν∇v(x, t))n(x) =

• n Γn × (0, T),

v(x, t) =

• n Γd × (0, T),

v(x, t) = gΓ(x, t)

• n Γg × (0, T),

v(x, 0) = v0(x) in Ω,

Channel Geometry and Grid

1 2 3 4 5 6 7 8 0.5 1

Figure: The channel geometry and coarse grid

EXAMPLE 1. y(t) =

• Ωobs

−∂x2v1(x, t) + ∂x1v2(x, t)dx

• ver the subdomain Ωobs = (1, 3) × (0, 1/2).

Model Reduction Results

nv np k 1352 205 13 5520 761 14 12504 1669 15 22304 2929 15

Table: Number nv of semidiscrete velocities v(t), number np of semidis- crete pressures p(t), and and size k of the reduced order velocities v(t) for various uniform refinements of the coarse grid.

Hankel S-vals and Convergence ADI

5 10 15 20 25 30 10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

The largest Hankel singular values

10 20 30 40 50 60 10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 Iteration Residual (Normalized) P Q

Convergence of the multishift ADI Algorithm

Figure: The left plot shows the largest Hankel singular values and the threshold τσ1. The right plot shows the normalized residuals Bk2 gen- erated by the multishift ADI Algorithm . for the approximate solution of the controllability Lyapunov equation (◦) and of the observability Lyapunov equation (∗).

Time (left) and Frequency (right) Responses

1 2 3 4 5 6 7 −150 −100 −50 50 100 150 Time full reduced 10

−4

10

−2

10 10

2

10

4

75.5 76 76.5 77 77.5 78 78.5 79 79.5 80 Frequency full reduced

Figure: Time response (left) and frequency response (right) for the full

• rder model (circles) and for the reduced order model (solid line).

Velocity Profile

full 22K dof reduced 15 dof

Figure: Velocities generated with the full order model (left col- umn) and with the reduced order model (right column) at t = 1.0996, 2.9845, 3.7699, 4.86965, 6.2832 (top to bottom).

Domain Decomposition Model Reduction

Systems with Local Nonlinearities

Linear Model in Large Domain Ω1

coupled to

Non-Linear Model in Small Domain Ω2

Linear Model provides B.C.’s to Non-Linear Model with R. Glowinski

D.C. Sorensen 41

Simple 1-D Model Problem

Ω1 ∪ Ω2 ∪ Ω3 Ω1 = (−10, −1), Ω2 = (−1, 1), Ω3 = (1, 10) Ω1 : Convection-Diffusion Ω2 : Burgers Equation Ω3 : Convection-Diffusion

D.C. Sorensen 42

Equations

ρk ∂yk ∂t (x, t) − µk ∂2yk ∂x2 (x, t) = Sk(x, t), (x, t) ∈ Ωk × (0, T), yk(x, 0) = yk0(x), x ∈ Ωk, k = 1, 3, ∂y1 ∂x (−10, t) = 0, ∂y3 ∂x (10, t) = 0 t ∈ (0, T), ρ2 ∂y2 ∂t (x, t) − µ2 ∂2y2 ∂x2 + y2 ∂y2 ∂x (x, t) = 0, (x, t) ∈ Ω2 × (0, T), y2(x, 0) = y20(x), x ∈ Ω2, with appropriate interface conditions

Dimension Reduction

Table: Dimension of the full and of the reduced order models for various discretization parameters N1, N2, N3 and τ = 10−4.

N1 = N3 N2 size of full model size of ROM 10 10 201 41 20 20 401 63 40 40 801 107 20 10 381 43 40 20 761 67

D.C. Sorensen 44

Time Response

5 10 15 2 4 t

Figure: Outputs 1, 2, 3 of the full order system corresponding to the discretization N1 = N2 = N3 = 10 are given by dotted, dashed and solid lines, respectively. Outputs 1, 2, 3 of the reduced order system are given by ∗, ◦ and , respecitively.

D.C. Sorensen 45

State Approximation

!10 !5 5 10 5 10 15 !5 5 10 x t !10 !5 5 10 5 10 15 !0.02 0.02 0.04 x t

Figure: Solution of the reduced order discretized PDE (left) and error between the solution of the discretized PDE and the reduced order system (right) for discretization N1 = N2 = N3 = 10.

D.C. Sorensen 46

Neuron Modeling Balanced Truncation on a Compartmental Neuron Model Steve Cox Tony Kellems Undergrads Nan Xiao and Derrick Roos

Quasi-active integrate and fire model. Complex model dimension 6000 Faithfully approximated with 10-20 variable ROM

D.C. Sorensen 47

Neuron Cell

90µm

Neuron Model Full Non-Linear Model

Ij,syn is the synaptic input into branch j aj 2Ri ∂xxvj = Cm∂tvj + GNam3

j hj(vj − ENa)

+ GKn4

j (vj − EK) + Gl(vj − El) + Ij,syn(x, t)

Kinetics of the potassium (n) and sodium (h, m) channels ∂tmj = αm(vj)(1 − mj) − βm(vj)mj ∂thj = αh(vj)(1 − hj) − βh(vj)hj ∂tnj = αn(vj)(1 − nj) − βn(vj)nj.

D.C. Sorensen 49

Cell Response - Lin and Non-Lin

5 10 15 20 25 30 −0.2 −0.1 0.1 0.2 0.3 0.4 0.5 0.6 v Morphology: RC−3−04−04−B.asc t (ms) 5 10 15 20 25 30 0.052 0.0525 0.053 0.0535 0.054 0.0545 0.055 0.0555 0.056 0.0565 0.057 m 5 10 15 20 25 30 0.589 0.59 0.591 0.592 0.593 0.594 0.595 0.596 0.597 0.598 h 5 10 15 20 25 30 0.3175 0.318 0.3185 0.319 0.3195 0.32 0.3205 0.321 0.3215 0.322 n Lin HH BT Full HH

Cell Response - Linear

5 10 15 20 25 30 −0.2 −0.1 0.1 0.2 0.3 0.4 0.5 0.6 v Morphology: RC−3−04−04−B.asc t (ms) 5 10 15 20 25 30 0.052 0.0525 0.053 0.0535 0.054 0.0545 0.055 0.0555 0.056 0.0565 0.057 m 5 10 15 20 25 30 0.589 0.59 0.591 0.592 0.593 0.594 0.595 0.596 0.597 0.598 h 5 10 15 20 25 30 0.3175 0.318 0.3185 0.319 0.3195 0.32 0.3205 0.321 0.3215 0.322 n Lin HH BT Full HH

Cell Response - Near Threshold

5 10 15 20 25 30 −1 −0.5 0.5 1 1.5 2 v Morphology: RC−3−04−04−B.asc t (ms) 5 10 15 20 25 30 0.045 0.05 0.055 0.06 0.065 0.07 m 5 10 15 20 25 30 0.565 0.57 0.575 0.58 0.585 0.59 0.595 0.6 0.605 h 5 10 15 20 25 30 0.315 0.32 0.325 0.33 0.335 0.34 n Lin HH BT Full HH

Neural ROM Results

◮ Interesting example of many-input , single-output system ◮ Simulation time single cell 14 sec Full vs .01 sec Reduced ◮ Ultimate goal is to simulate a few-Million neuron system

• ver a minute of brain-time

◮ Currently limited to a 10K neuron system

• ver a few brain-seconds

◮ Parallel computing required D.C. Sorensen 53

Summary

CAAM TR07-17, SPSVD, M. Shah & DCS, Also SIMAX 06 CAAM TR07-02, M. Heinkenschloss, DCS, & K. Sun CAAM TR07-14, K. Sun, R. Glowinski, M. Heinkenschloss, DCS Gramian Based Model Reduction: Balanced Reduction Solving Large Lyapunov Equations: Approximate Balancing 1M vars now possible Balanced Reduction of Oseen Eqns: Extension to Descriptor System Multi-Shift ADI Without Explicit Projectors: Only need Saddle Point Solver Sparse Direct or Iterative Domain Decomposition - Systems with Local Nonlinearities Neural Modeling - Single Cell ROM ⇒ Many Interactions