[PPT] - Krylov Model-Order Reduction Techniques for Time- and PowerPoint Presentation

SLIDE 1

Krylov Model-Order Reduction Techniques for Time- and Frequency-Domain Wavefjeld Problems

Rob Remis Delft University of Technology

November 6 – 10, 2017 – ICERM Brown University 1

SLIDE 2

Acknowledgment

This is joint work with Mikhail Zaslavsky, Schlumberger-Doll Research Jörn Zimmerling, Delft University of Technology Vladimir Druskin, Schlumberger-Doll Research

2

SLIDE 3

Outline

Basic equations and symmetry properties Polynomial Krylov reduction 1 Perfectly matched layers Polynomial Krylov reduction 2 Laurent polynomial (extended) Krylov reduction Introducing rational and preconditioned rational Krylov reduction

More on this in J. Zimmerling’s talk

3

SLIDE 4

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Basic equations

First-order wavefjeld system

acoustics, seismics, electrodynamics

(D + M∂t) F = −w(t)Q Plus initial conditions Dirichlet boundary conditions Signature matrix δ− δ− = diag(1, −1, −1, −1) (acoustics)

4

SLIDE 5

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Basic equations

Spatial discretization (D + M∂t) f = −w(t)q Order of this system can be very large especially in 3D Discretized counterpart of δ− is denoted by d−

5

SLIDE 6

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Basic equations

Solution f (t) = −w(t) ∗ η(t) exp(−At)M−1q η(t) Heaviside unit step function System matrix A = M−1D

6

SLIDE 7

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Symmetry

System matrix A is skew-symmetric w.r.t. WM Evolution operator exp(−At) is orthogonal w.r.t. WM Inner product and norm x, y = yHWMx x = x, x1/2 Stored fjeld energy in the computational domain

(sum of fjeld energies)

1 2f 2 Initial-value problem: norm of f is preserved

7

SLIDE 8

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Symmetry

System matrix A is symmetric w.r.t. WMd− Bilinear form x, y = yHWMd−x Free fjeld Lagrangian

(difgerence of fjeld energies)

1 2f , f Symmetry property related to reciprocity

8

SLIDE 9

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Symmetry

Introduce dp = 1 2(I + d−) and dm = 1 2(I − d−) Write f = f (q) to indicate that the fjeld is generated by a source q Reciprocity Source vector: q = dpq, receiver vector r = dpr f (q), r = q, f (r) Source vector: q = dpq, receiver vector r = dmr f (q), r = −q, f (r)

9

SLIDE 10

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Polynomial Krylov reduction

Exploit symmetry of system matrix in a Lanczos reduction algorithm For lossless media both symmetry properties lead to the same reduction algorithm First symmetry property is lost for lossy media with a system matrix of the form A = M−1(D + S) Second symmetry property still holds

10

SLIDE 11

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Polynomial Krylov reduction

For lossless media, FDTD can be written in a similar form as Lanczos algorithm recurrence relation for FDTD = recurrence relation for Fibonacci polynomials Stability of FDTD and numerical dispersion can be studied using this connection

11

SLIDE 12

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Polynomial Krylov reduction

Lanczos recurrence coeffjcients: βi Comparison with FDTD: 1/βi = time step of Lanczos Automatic time step adaptation – no Courant condition Lanczos reduction hardly provides any speedup compared with FDTD Both are polynomial fjeld approximations

Lanczos: fjeld is approximated by a Lanczos polynomial in A FDTD: fjeld is approximated by a Fibonacci polynomial in A

12

SLIDE 13

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

PML

No outward wave propagation has been included yet Implementation via Perfectly Matched Layers (PML) Coordinate stretching (Laplace domain) ∂k ← → χ−1

k ∂k

k = x, y, z Stretching function χk(k, s) = αk(k) + βk(k) s

13

SLIDE 14

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

PML

Stretched fjrst-order system

D(s) + S + sM

ˆ F = −ˆ w(s)Q Direct spatial discretization

D(s) + S + sM

ˆ f = −ˆ w(s)q Leads to nonlinear eigenproblems for spatial dimensions > 1

14

SLIDE 15

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

PML

Linearization of the PML Spatial fjnite-difgerence discretization using complex PML step sizes (Dcs + S + sM) fcs = −w(s)q System matrix Acs = M−1(Dcs + S)

V. Druskin and R. F. Remis, “A Krylov stability-corrected coordinate stretching method to simulate wave

propagation in unbounded domains,” SIAM J. Sci. Comput., Vol. 35, 2013, pp. B376 – B400.

V. Druskin, S. Güttel, and L. Knizhnerman, “Near-optimal perfectly matched layers for indefjnite

Helmholtz problems,” SIAM Rev. 58-1 (2016), pp. 90 – 116. 15

SLIDE 16

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

PML

Spectrum of the system matrix Acs

Re(λ) Im(λ) Lossless resonator

16

SLIDE 17

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

PML

Eigenvalues move into the complex plane

Re(λ) Im(λ) Complex scaling

17

SLIDE 18

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

PML

Stable part of the spectrum

Re(λ) Im(λ) Stable part

18

SLIDE 19

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

PML

Stability correction

Re(λ) Im(λ) Stable part

19

SLIDE 20

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

PML

Time-domain stability-corrected wave function f (t) = −w(t) ∗ 2η(t)Re

η(Acs) exp(−Acst)q
Complex Heaviside unit step function

η(z) =

1

Re(z) > 0 Re(z) < 0

20

SLIDE 21

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

PML

Frequency-domain stability-corrected wave function ˆ f (s) = −ˆ w(s)

r(Acs, s) + r(¯

Acs, s)

q

with r(z, s) = η(z) z + s Note that ˆ f (¯ s) = ¯ ˆ f (s) and the stability-corrected wave function is a nonentire function of the system matrix Acs

21

SLIDE 22

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

PML

Symmetry relations are preserved With a step size matrix W that has complex entries These entries correspond to PML locations

22

SLIDE 23

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Polynomial Krylov Reduction

Stability-corrected wave function cannot be computed by FDTD SLDM fjeld approximations via modifjed Lanczos algorithm Reduced-order model fm(t) = −w(t) ∗ 2M−1qη(t)Re [Vmη(Hm) exp(−Hmt)e1]

23

SLIDE 24

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Polynomial Krylov Reduction

m = 300

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10

−13

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 x 10

8

Time [s] Electric Field Strength [V/m]

24

SLIDE 25

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Polynomial Krylov Reduction

m = 400

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10

−13

−8 −6 −4 −2 2 4 6 8 x 10

7

Time [s] Electric Field Strength [V/m]

25

SLIDE 26

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Polynomial Krylov Reduction

m = 500

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10

−13

−8 −6 −4 −2 2 4 6 8 x 10

7

Time [s] Electric Field Strength [V/m]

26

SLIDE 27

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Extended Krylov reduction

Stability-corrected wave function is approximated by a Lanczos polynomial in Acs The wave function is a nonentire function of the system matrix Idea: approximate the stability-corrected function by a Laurent polynomial Perhaps an even better idea: approximate the stability-corrected wave function by rational functions

(more on this later)

27

SLIDE 28

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Extended Krylov reduction

Extended Krylov subspace Km1,m2 = span{A−m1+1q, ..., A−1q, q, Aq, ..., Am2−1q} Elements from this space: Laurent polynomials in matrix A acting on the source vector q

28

SLIDE 29

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Extended Krylov reduction

Original extended Krylov method of Druskin and Knizhnerman generates the sequence of subspaces Km1,1 ⊂ Km1,2 ⊂ ... ⊂ Km1,m2 via short-term recurrence relations A more general approach was proposed by Jagels and Reichel Effjciently generate the sequence of subspaces K1,i+1 ⊂ K2,2i+1 ⊂ ... ⊂ Kk,ki+1 again via short term recurrence relations i is an integer # matvec with A = i · # matvec with A−1

29

SLIDE 30

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Extended Krylov reduction

SEG Salt model/velocity profjle – 3D acoustics, frequency-domain, order ≈ 93 million

30

SLIDE 31

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Extended Krylov reduction

SEG Salt model - generalized EKS implementation

1000 2000 3000 4000 5000 6000 7000 8000 10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

Iteration number L2 Error Frequency range 2.5Hz to 7.5Hz

+/−=3 +/−=5 +/−=7 +/−=Inf

31

SLIDE 32

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Extended Krylov reduction

SEG Salt model - generalized EKS implementation

2000 4000 6000 8000 10000 10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

Iteration number L2 Error Frequency range 1Hz to 7.5Hz

+/−=3 +/−=5 +/−=7 +/−=Inf

32

SLIDE 33

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Extended Krylov reduction

SEG Salt model - generalized EKS implementation

2000 4000 6000 8000 10000 10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

Iteration number L2 Error Frequency range 0.5Hz to 7.5Hz

+/−=3 +/−=5 +/−=7 +/−=Inf

33

SLIDE 34

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Extended Krylov reduction

SEG Salt model - generalized EKS implementation

2000 4000 6000 8000 10000 10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

Frequency range 0.1Hz to 7.5Hz Iteration number L2 Error

+/−=3 +/−=5 +/−=7 +/−=Inf

34

SLIDE 35

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Rational Krylov reduction

In a rational Krylov method, we approximate the fjeld by the span of snapshots ˆ f (s1), ˆ f (s2), ..., ˆ f (sm) for difgerent frequencies si, i = 1, 2, ..., m The snapshots are obtained by solving discretized wavefjeld systems of the form [D(si) + siM] ˆ f (si) = −q

35

SLIDE 36

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Rational Krylov reduction

No PML linearization is necessary! When using a rational Krylov method, we deal with the above system directly and not with the stability-corrected system/wave function

36

SLIDE 37

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Rational Krylov reduction

Return to the system [D(s) + sM] f = −q Symmetry relation DT(s) ˜ W (s) = ˜ W (s)D(s) ˜ W (s) = W (s)d− is a nonsingular diagonal s-dependent step size matrix

37

SLIDE 38

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Rational Krylov approximations

Multiply by ˜ W to obtain A(s)ˆ f (s) = ˜ q System matrix A(s) = ˜ W (s) [D(s) + sM] Properties AT(s) = A(s) and A∗(s) = A(s∗)

38

SLIDE 39

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Rational Krylov approximations

Solve system for m ≥ 1 difgerent frequencies Construct the subspace Km = span{ˆ f (s1), ˆ f (s2), ..., ˆ f (sm)} and take KR

m = span{Re Km, Im Km}

as an expansion and projection space Note that ˆ f (si) ∈ KR

m

and ˆ f (s∗

i ) ∈ KR m

i = 1, 2, ..., m

39

SLIDE 40

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Rational Krylov approximations

Let Vm be a basis matrix of KR

m

Field approximation fm(s) = Vmˆ am(s) Expansion coefgcients are determined from Galerkin condition Structure-preserving reduced-order model ˆ fm(s) = Vm ˆ R−1

m (s)V T m ˜

q Rm(s) = V T

m A(s)Vm

40

SLIDE 41

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Rational Krylov approximations

Interpolation properties ˆ fm(si) = ˆ f (si) and ˆ fm(s∗

i ) = ˆ

f (s∗

i )

i = 1, 2, ..., m For coinciding source/receiver pairs d ds ˜ qTˆ fm(s)

s=si,s∗

i

= d ds ˜ qTˆ f (s)

s=si,s∗

i

i = 1, 2, ..., m

41

SLIDE 42

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Rational Krylov approximations

Large travel times: frequency-domain wavefjeld highly

scillatory in frequency domain

Rational Krylov method requires many sampling/interpolation points Phase-preconditioned rational Krylov method: factor out the strongly oscillating part using high-frequency asymptotics Much more on this in talk of J. Zimmerling

42

SLIDE 43

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Literature

Complex PML step sizes and stability-correction:

V. Druskin and R. F. Remis, “A Krylov stability-corrected coordinate stretching method to simulate

wave propagation in unbounded domains,” SIAM J. Sci. Comput., Vol. 35, 2013, pp. B376 – B400.

V. Druskin, S. Güttel, and L. Knizhnerman, “Near-optimal perfectly matched layers for indefjnite

Helmholtz problems,” SIAM Rev. 58-1 (2016), pp. 90 – 116.

J. Zimmerling, L. Wei, P. Urbach, and R. Remis, “A Lanczos model-order reduction technique to

effjciently simulate electromagnetic wave propagation in dispersive media,” J. Comp. Phys.,

Vol. 315, 2016, pp. 348 – 362.

43

SLIDE 44

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Literature

Extended Krylov subspace method:

V. Druskin and L. Knizhnerman, “Extended Krylov subspaces: approximation of the matrix square

root and related functions,” SIAM J. Matrix Anal. Appl., Vol. 19, 1998, pp. 755 – 771.

C. Jagels and L. Reichel, “Recursion relations for the extended Krylov subspace method,” Linear

Algebra Appl., Vol. 434, pp. 1716 – 1732, 2011.

V. Druskin, R. Remis, and M. Zaslavsky, “An extended Krylov subspace model-order reduction

technique to simulate wave propagation in unbounded domains,” J. Comp. Phys., Vol. 272, 2014,

pp. 608 – 618.

44

SLIDE 45

Basic equations & symmetry Polynomial Krylov reduction 1 PML Polynomial Krylov Reduction 2 Extended Krylov Reduction (Preconditioned) rational Krylov methods

Literature

Phase-preconditioned rational Krylov method

V. Druskin, R. Remis, M. Zaslavsky, and J. Zimmerling, “Compressing large-scale wave

propagation via phase-preconditioned rational Krylov subspaces,” to appear on ArXiv, 2017. See also Jörn Zimmerling’s talk. 45