Computing the quasipotential for nongradient SDEs in 3D Samuel F. - - PowerPoint PPT Presentation

Computing the quasipotential for nongradient SDEs in 3D

Samuel F. Potter

Joint work with Shuo Yang and Maria K. Cameron

July 22, 2019

Setting

System evolving according to SDE (Itˆ

• ):

dx = b(x)dt + √ǫdw, x ∈ R3 standard Brownian motion small noise parameter C1 vector field For attractor A ⊆ R3, compute quasipotential: U(x) = min

ψ,L

L

0 (b(ψ)ψ′ − b(ψ) · ψ′) ds

geometric action

(obtained from Freidlin-Wentzell action)

Only makes sense to compute U if b is nongradient

1/29

attractor A

Computing the Quasipotential on a Grid

goal: march U forward causally in semi-Lagrangian fashion (approximate locally minimum action paths) U|A ≡ 0 2/29

Dijkstra-like Solvers for HJ Equations

Computing quasipotential involves solving Hamilton-Jacobi equation of the form: F

• x,

∇u(x) ∇u(x)

• ∇u(x) = 1

Previous approach: ordered upwind method (OUM)—requires bounded anisotropy (Sethian & Vladimirsky ‘01): Υ = Fmax/Fmin May be unbounded—OUM for quasipotential handles this (Cameron ‘12). Additional error is O(h2) OLIMs (Dahiya and Cameron ‘17) reduce error constant by 102–103 by higher-order quadrature rules

3/29

Description of the Algorithm

For grid points x, four states: attractor A unknown considered accepted front accepted

4/29

Description of the Algorithm

Update considered neighbors of new accepted front point attractor A unknown considered accepted front accepted

Kh 4/29

Description of the Algorithm

Update considered neighbors of new accepted front point attractor A unknown considered accepted front accepted

4/29

Motivation

Straightforward extension of the OUM (Sethian & Vladimirsky ‘01) to the quasipotential in 3D is possible but far too slow Hence, the OUM-based 2D quasipotential (Cameron ‘12) suffers from the same problem The 2D quasipotential OLIM (Dahiya & Cameron ‘17) can be extended to 3D thanks to its hierarchical update strategy Despite this, the extension of this method 3D on a mesh of size 6123 was taking days. . .

5/29

Contribution

Thanks to some technical innovations, the runtime of the 3D quasipotential OLIM was reduced from days to hours: The hierarchical update strategy was enhanced The hierarchical updates were combined with update skipping using the KKT optimality conditions Improving the process of simplex selection reduced the

• verall number of updates significantly

6/29

Simplex updates

In 3D, three types of updates: line, triangle, tetrahedron E.g., for the tetrahedron update:

Q3(x0, x1, x2, x) = min

λ {Uλ + Q(xλ, x)}, where

Uλ = U(x0) + λ1(U(x1) − U(x0)) + λ(U(x2) − U(x0)) xλ = x0 + λ1(x1 − x0) + λ2(x2 − x0) Q(xλ, x) = bmλx − xλ − bmλ · (x − xλ) bmλ = bm0 + λ1(bm1 − bm0) + λ2(bm2 − bm0) bmi = b xi + x 2

• ,

i = 0, 1, 2, subject to λ1 ≥ 0, λ2 ≥ 0, λ1 + λ2 ≤ 1. 7/29

Simplex update minimization problem

Cost function for simplex update: f(λ) = Uλ + x − xλbmλ − (x − xλ) · bmλ ∇f = δU + x − xλ bmλ B⊤bmλ + bmλ x − xλX⊤(x − xλ) − B⊤(x − xλ) − X⊤bmλ H = X⊤(x − xλ)

• B⊤bmλ

⊤ +

• X⊤(x − xλ)
• B⊤bmλ

⊤⊤ x − xλbmλ + x − xλ bmλ B⊤B + bmλ x − xλX⊤X − x − xλ bmλ3 B⊤bmλ

• B⊤bmλ
• −

bmλ x − xλX⊤(x − xλ)

• X⊤(x − xλ)

⊤ − B⊤X −

• B⊤X

⊤ Where: B = b(xm1) − b(xm0) b(xm2) − b(xm0) X = x0 − x1 x0 − x2

• δU = (U(x1) − U(x0), U(x2) − U(x0))

8/29

Simplex update minimization problem

Cost function for simplex update: f(λ) = Uλ + x − xλbmλ − (x − xλ) · bmλ ∇f = δU + x − xλ bmλ B⊤bmλ + bmλ x − xλX⊤(x − xλ) − B⊤(x − xλ) − X⊤bmλ H = X⊤(x − xλ)

• B⊤bmλ

⊤ +

• X⊤(x − xλ)
• B⊤bmλ

⊤⊤ x − xλbmλ + x − xλ bmλ B⊤B + bmλ x − xλX⊤X − x − xλ bmλ3 B⊤bmλ

• B⊤bmλ
• −

bmλ x − xλX⊤(x − xλ)

• X⊤(x − xλ)

⊤ − B⊤X −

• B⊤X

⊤ Where: B = b(xm1) − b(xm0) b(xm2) − b(xm0) X = x0 − x1 x0 − x2

• δU = (U(x1) − U(x0), U(x2) − U(x0))

away from equilibria of b(x), H is positive definite for h small enough if K ∼ h−α, where 0 < α < 1

H

8/29

A fast search for admissible updates

Neighborhoods used to organize different types of updates: N1(x0) =

• x ∈ Z3 : x − x01 = 1
• N2(x0) =
• x ∈ Z3 : x − x01 = 2 and x − x0∞ = 1
• N3(x0) =
• x ∈ Z3 : x − x01 = 3 and x − x0∞ = 1
• N(x0) = N1(x0) ∪ N2(x0) ∪ N3(x0)

N K far(x0) ≈

• x ∈ Z3 : x − x02 ≤ K
• 9/29

A fast search for admissible updates

Neighborhoods used to organize different types of updates: N1(x0) =

• x ∈ Z3 : x − x01 = 1
• N2(x0) =
• x ∈ Z3 : x − x01 = 2 and x − x0∞ = 1
• N3(x0) =
• x ∈ Z3 : x − x01 = 3 and x − x0∞ = 1
• N(x0) = N1(x0) ∪ N2(x0) ∪ N3(x0)

N K far(x0) ≈

• x ∈ Z3 : x − x02 ≤ K
• considered point

being updated 9/29

A fast search for admissible updates

Motivation for using these simplexes: interpolation error

E.g., assume U(x) = U0 + g⊤x + 1

2x⊤Hx

In 1D, on interval of size h, max interpolation error = 1

8h2H

In 2D or higher, more complicated, but same trend: small h = ⇒ smaller interpolation error Main idea: use minimizing line update to guide search for area which is to be refined into triangle and tetrahedron updates 10/29

Skipping updates using the KKT conds.

Constrained minimization problem in canonical form: minimize f(λ) = Uλ + bmλx − xλ − bmλ · (x − xλ) subject to λ1 ≥ 0 λ2 ≥ 0 1 − λ1 − λ2 ≥ 0 Lagrangian: L(λ, µ) = f(λ) − µ1λ1 − µ2λ2 − µ3(1 − λ1 − λ2) KKT conditions: ∇λL = ∇f − µ1 1

• − µ2

1

• − µ3

−1 −1

• = 0,

µi ≥ 0 for i = 1, 2, 3, λ1 ≥ 0, λ2 ≥ 0 1 − λ1 − λ2 ≥ 0, λ1µ1 = 0, λ2µ2 = 0 , (1 − λ1 − λ2)µ3 = 0. 11/29

Skipping updates using the KKT conds.

Constrained minimization problem in canonical form: minimize f(λ) = Uλ + bmλx − xλ − bmλ · (x − xλ) subject to λ1 ≥ 0 λ2 ≥ 0 1 − λ1 − λ2 ≥ 0 Lagrangian: L(λ, µ) = f(λ) − µ1λ1 − µ2λ2 − µ3(1 − λ1 − λ2) KKT conditions: ∇λL = ∇f − µ1 1

• − µ2

1

• − µ3

−1 −1

• = 0,

µi ≥ 0 for i = 1, 2, 3, λ1 ≥ 0, λ2 ≥ 0 1 − λ1 − λ2 ≥ 0, λ1µ1 = 0, λ2µ2 = 0 , (1 − λ1 − λ2)µ3 = 0. to verify, ends up requiring only simple comparisons involving ∇f 11/29

Skipping updates using the KKT conds.

Order updates from low dimension to high dimension

Start with line updates, follow with triangle updates, then tetra updates Lower dim updates are incident on the boundaries of higher dim updates When doing a triangle update, if we get an interior point solution, since f is strictly convex for h small enough (and with properly chosen K), we can determine if λ∗ is the constrained optimum of higher dim adjacent tetrahedron updates Important part of getting a speed-up in 3D and verified numerically x0 x x1 xλ∗ x2 12/29

The hierarchical update strategy

A brute force update algorithm:

• 1. x0 ← considered point with smallest U value
• 2. x0 becomes accepted front
• 3. accepted front points in N(x0) with no considered

neighbors become accepted

• 4. for all considered x ∈ N K

far(x0), set U(x) to minimum Q3(x, x0, x1, x2) value where x1, x2 are accepted front points in N1(x0) ∪ N2(x0)

• 5. for all unknown x ∈ N(x0)

(a) x becomes considered (b) for all accepted front y0 ∈ N K far(x), set U(x) to minimum Q3(y0, y1, y2, x) value where y1, y2 are accepted front points in N1(y0) ∪ N2(y0)

13/29

A brute force update algorithm:

• 1. x0 ← considered point with smallest U value
• 2. x0 becomes accepted front
• 3. accepted front points in N(x0) with no considered

neighbors become accepted

• 4. for all considered x ∈ N K

far(x0), set U(x) to minimum Q3(x, x0, x1, x2) value where x1, x2 are accepted front points in N1(x0) ∪ N2(x0)

• 5. for all unknown x ∈ N(x0)

(a) x becomes considered (b) for all accepted front y0 ∈ N K far(x), set U(x) to minimum Q3(y0, y1, y2, x) value where y1, y2 are accepted front points in N1(y0) ∪ N2(y0)

The hierarchical update strategy

this is the most expensive step—in this form, it leads to a huge number of updates (line, tri, and tetra) 13/29

A brute force update algorithm:

• 1. x0 ← considered point with smallest U value
• 2. x0 becomes accepted front
• 3. accepted front points in N(x0) with no considered

neighbors become accepted

• 4. for all considered x ∈ N K

far(x0), set U(x) to minimum Q3(x, x0, x1, x2) value where x1, x2 are accepted front points in N1(x0) ∪ N2(x0)

• 5. for all unknown x ∈ N(x0)

(a) x becomes considered (b) for all accepted front y0 ∈ N K far(x), set U(x) to minimum Q3(y0, y1, y2, x) value where y1, y2 are accepted front points in N1(y0) ∪ N2(y0)

The hierarchical update strategy

this is the most expensive step—in this form, it leads to a huge number of updates (line, tri, and tetra)

goal: use continuity of f in terms of origin of local characteristic to guide search & combine with update skipping

13/29

The hierarchical update strategy

The cost function: f(xλ) = U(xλ)+

• b

x + xλ 2

• x−xλ−b

x + xλ 2

• ·(x−xλ)

continuously depends on xλ, the base of the characteristic arriving at x. Find minimum line update, and then do surrounding triangle and tetrahedron updates to significantly reduce the number

• f higher dim updates

Number of line updates still depends on K, but number of tri and tetra updates is now O(1) with respect to K, big savings

14/29

The cost function: f(xλ) = U(xλ)+

• b

x + xλ 2

• x−xλ−b

x + xλ 2

• ·(x−xλ)

continuously depends on xλ, the base of the characteristic arriving at x. Find minimum line update, and then do surrounding triangle and tetrahedron updates to significantly reduce the number

• f higher dim updates

Number of line updates still depends on K, but number of tri and tetra updates is now O(1) with respect to K, big savings

The hierarchical update strategy

combine this idea with using the KKT conditions to skip higher dim updates—reduce the number of tri and tetra updates even more (reject 50%–70% of leftover tri and tetra updates)

14/29

Numerical tests

Five examples to test accuracy: Three examples using linear SDEs Two nonlinear examples with an analytically available quasipotential Two examples by Tao (‘18)—two point attractors with a hyperbolic periodic orbit serving as a transition state A genetic switch model w/ two asymptotically stable equilibria

Solver affected by relative magnitudes of rotational (l(x) = 1

2∇U(x) + b(x)) and potential components of b(x) (− 1 2∇U(x))

Define: ξ(x) = l(x)/ 1

2∇U(x)

15/29

Numerical tests (experimental setup)

Problem size:

• mesh size: N 3
• grid point spacing: h ∼ CN −1
• N = 2p + 1, where p = 5, 6, 7, 8, 9
• i.e., N = 33, 65, 129, 257, 513

The parameter K ranges from 1 to 30 Workstation:

• 4.2 GHz quad-core Intel Core i7 (run on single thread)
• 64 GB 2400MHz DDR4 SODIMM SDRAM

16/29

Numerical tests (examples)

Example 1 (a basic test)

b(x, y, z) =   −3 −4 −1 3 −4 −1 3 4 −1     x y z   U(x) = x⊤Qx Q = diag(3, 4, 1) ξ(x) ∈ [0, √ 3] Note: eigenvectors of Q aligned with coordinate axes (Ω = [−1, 1]3) Domain: Vector field: Exact quasipotential: Quasipotential matrix: Ratio of rotational component to potential:

17/29

Numerical tests (examples)

Examples 2 and 3 (parametrized ξ)

J =   −1

−1 2

−ρ ρ − 1

2

  , R = Rx1( 2π

3 )Rx2( π 8 )Rx3( π 5 )

ξ(x) ∈ [0, 2ρ] Domain: Ω = [−1, 1]3 SDE: dx = R⊤JRx + √ǫdw where: Exact quasipotential: x⊤R⊤QRx, where Q = diag(1, 1

2, 1 2)

Ratio of rotational component to potential: Note: eigenvectors of Q unaligned with coordinate axes

18/29

Numerical tests (examples)

Examples 2 and 3 (parametrized ξ)

J =   −1

−1 2

−ρ ρ − 1

2

  , R = Rx1( 2π

3 )Rx2( π 8 )Rx3( π 5 )

ξ(x) ∈ [0, 2ρ] Domain: Ω = [−1, 1]3 SDE: dx = R⊤JRx + √ǫdw where: Exact quasipotential: x⊤R⊤QRx, where Q = diag(1, 1

2, 1 2)

Ratio of rotational component to potential: Note: eigenvectors of Q unaligned with coordinate axes

Example 2: ρ = √ 3/2 = ⇒ ξ ∈ [0, √ 3] Example 3: ρ = 5 = ⇒ ξ ∈ [0, 10] 18/29

Numerical tests

Examples 4 and 5 (double-well potential)

Domain: Ω = [−2, 0] × [−1, 1]2 Vector field: b(x) =   −2(x3

1 − x1) − ρ(x2 + x3)

−x2 + 2ρ(x3

1 − x1)

−x3 + 2ρ(x3

1 − x1)

  Ratio ξ(x) = ρ everywhere except the equilibria Two asymptotically stable equilibria: (±1, 0, 0) One Morse index 1 saddle: (0, 0, 0) Exact quasipotential within level set passing through the

• rigin: U(x) = x4

1 − 2x2 1 + x2 2 + x2 3 + 1

19/29

Exact quasipotential within level set passing through the

• rigin: U(x) = x4

1 − 2x2 1 + x2 2 + x2 3 + 1

Two asymptotically stable equilibria: (±1, 0, 0)

Numerical tests

Examples 4 and 5 (double-well potential)

Domain: Ω = [−2, 0] × [−1, 1]2 Vector field: b(x) =   −2(x3

1 − x1) − ρ(x2 + x3)

−x2 + 2ρ(x3

1 − x1)

−x3 + 2ρ(x3

1 − x1)

  Ratio ξ(x) = ρ everywhere except the equilibria One Morse index 1 saddle: (0, 0, 0) Example 4: ρ = 1 Example 5: ρ = 10

19/29

Numerical tests

Example 1

20/29

Numerical tests

Example 4

21/29

Numerical tests

Example 5

22/29

Numerical tests

A guideline for choosing K for mesh size N 3

N 33 65 129 257 513 K 4 6 8 10 14

Example #1 #2 #3 T vs K (N = 513) [s] 31.4K1.79 8.57K1.80 93.8K1.78 E∞ vs h (K opt.) 0.181h0.750 0.0808h0.658 12.3h1.49 T vs N (K opt.) [ns] 27.4N 4.05 61.0N 4.11 131.0N 3.98 Example #4 #5 T vs K (N = 513) [s] 66.3K1.80 1.94K1.80 E∞ vs h (K opt.) 0.635h0.989 19.4h1.54 T vs N (K opt.) [ns] 49.2N 4.10 92.6N 4.00

Least squares fits

23/29

Numerical tests

24/29

Numerical tests

Example 6 (systems with hyperbolic periodic orbits)

25/29

Numerical tests

Example 7 (a genetic switch model)

26/29

Conclusions

Extension of previous OLIM in 2D to 3D, preserving higher accuracy compared to OUM-based approach Presented a hierarchical update strategy which reduces runtime in 3D on large mesh from days to hours Conducted an extensive numerical study validating the proposed solver on many examples Present guidelines for choosing parameter K

27/29

Future Work

Extend solver to efficiently handle systems where stochastic dynamics are effectively restricted to a neighborhood of a low-dimensional manifold (e.g. genetic switch model) Explore methods of computing the prefactor densely (see the Maier & Stein papers for related ideas) Computing the prefactor requires ∆U—can a quasipotential solver be made higher-order accurate? There are several approaches for solving static Hamilton-Jacobi equations in the literature, but they would likely need to be modified to compute the quasipotential

28/29

References

Computing the quasipotential for nongradient SDEs in 3D, Shuo Yang, Samuel F. Potter, and Maria K. Cameron, Journal of Computational Physics, 379 (2019) 325-350, https://doi.org/10.1016/j.jcp. 2018.12.005, arXiv: 1808.00562 Ordered Line Integral Methods for Solving the Eikonal Equation, Samuel

• F. Potter and Maria K. Cameron, in revision (2019), arXiv:1902.06825

Computing the quasipotential for highly dissipative and chaotic SDEs. An application to stochastic Lorenz’63, Maria Cameron and Shuo Yang, CAMCoS, accepted (2019), arXiv:1809.09987v2 An Ordered Line Integral Method for Computing the Quasi-potential in the case of Variable Anisotropic Diffusion, Daisy Dahiya and Maria Cameron, 2018, Physica D 382-383 (2018), 33–45, https://doi.org/10.1016/j.physd. 2018.07.002, arXiv:1806.05321 29/29