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 Maria Cameron and Shuo Yang University of Maryland, College Park, MD Mid-Atlantic Numerical Analysis Day, 2019

What are rare events in SDEs? d x t = b ( x t ) d t + √ ǫ d w t C 1 vector field Brownian motion small parameter stable limit cycles stable equilibria more complex attractors

What do we want to estimate? expected escape times from basins of attraction maximum likelihood escape paths quasi-invariant probability densities in the neighborhood of attractors

Gradient SDEs 2 β − 1 d w t � d x t = −∇ V ( x t ) d t + Saddle Min #1 Min #2

Gradient SDEs Arrhenius equation (1884) � � V saddle − V min reaction rate ∝ exp k B T Gibbs measure � � f ( x ) = Z − 1 exp − V ( x ) k B T Kramer’s/Langer’s formula (1940) (1969) � reaction rate ≅ | λ | | det H saddle | e − β ( V saddle − V min ) det H min 2 π

Large deviation theory Friedlin-Wentzell action functional transition state � T 0 � ˙ S T ( φ ) = 1 φ − b ( φ ) � 2 dt 2 Quasipotential attractor A U A ( x ) = inf φ,T { S T ( φ ) | φ (0) ∈ A, φ ( T ) = x } basin of attraction B ( A ) Expected escape time 1 β log E [ τ exit ] = 1 lim β →∞ 2 inf x ∈ ∂B ( A ) U A ( x ) basin of attraction of A

Computing the quasipotential Geometric action (Heymann and Vanden-Eijnden, 2008) � � � L � ψ ′ �� b ( ψ ) � − ψ ′ · b ( ψ ) S ( ψ ) = ds 0 Hamilton-Jacobi equation Nonlinear PDE for the quasipotential �∇ U A � 2 + 2 b ( x ) · ∇ U A = 0 , U A ( x ) = 0 , x ∈ A Minimum action paths Orthogonal decomposition of b b ( x ) = −∇ U ( x ) / 2 + ℓ ( x ) d ψ ∗ d s � ∇ U ( ψ ∗ ) + b ( ψ ∗ ) ℓ ( x ) = ∇ U ( x ) / 2 + b ( x ) ℓ ( x ) ⊥ ∇ U ( x )

Computing the quasipotential on a mesh: goal Develop numerical methods for computing: ◮ Quasipotential: U A ( x ) ◮ Minimum action paths: ψ ∗ What can we do with U A ? ◮ More complete information about the dynamics in a region of interest ◮ A useful visualization tool

Computing the quasipotential on a mesh: key ideas The ordered upwind method for solving HJ PDEs (Sethian and Vladimirsky, 2001, 2003) � � ∇ U �∇ U ( x ) � = 1 F x, �∇ U � U ( x, y ) where 0 < F min ≤ F ≤ F max < ∞ y Key idea and motivation Dynamic programming principle Motivated by Sethian’s fast marching method (1996) for solving the eikonal equation: F ( x ) �∇ U � = 1 x

Computing the quasipotential on a mesh: previous methods OUM for the quasipotential (Cameron, 2012) �∇ U � 2 + 2 b ( x ) · ∇ U = 0 = 1 ⇒ �∇ U � �∇ U � = 1 − 2 b ( x ) · ∇ U Ordered line integral methods unbounded (Dahiya and Cameron, 2017) faster convergence and Solve the minimization problem directly: 100–1000 times more accurate than OUM � L � ψ ′ �� b ( ψ ) � − ψ ′ · b ( ψ ) � � S ( ψ ) = ds 0 � � � 1.5–4 times faster U A ( x ) = inf ψ S ( ψ ) � ψ (0) ∈ A, ψ ( L ) = x than OUM Apply a time saving hierarchical update strategy

play video

Ordered line integral methods in 3D Developed a 3D solver that: ◮ is 1–2 orders of magnitude faster than the OUM extended to 3D, ◮ improves the 2D hierarchical update strategy, ◮ skips higher-dimensional updates using the KKT conditions for constrained optimization, ◮ and searches for updates intelligently, ◮ without compromising accuracy.

Ordered line integral methods in 3D: simplex updates Three types of updates needed: line , triangle , and tetrahedron updates x x x m 0 x m 2 x m x m 1 1 x m x m x m 0 x 0 λ λ x 2 x 1 x λ x 0 x 1

Ordered line integral methods in 3D: local minimization problem E.g., the minimization problem solved to compute U ( x ): minimize U λ + Q ( x λ , x ) subject to λ 1 ≥ 0 , λ 2 ≥ 0 , λ 1 + λ 2 ≤ 1 Where: U λ = U ( x 0 ) + λ 1 ( U ( x 1 ) − U ( x 0 )) + λ 2 ( U ( x 2 − x 0 )) x λ = x 0 + λ 1 ( x 1 − x 0 ) + λ 2 ( x 2 − x 0 ) Q ( x λ , x ) = � b m λ �� x − x λ � − b m λ · ( x − x λ ) b m λ = b m 0 + λ 1 ( b m 1 − b m 0 ) + λ 2 ( b m 2 − b m 0 ) � x + x i � b m i = b i = 0 , 1 , 2 . , 2

Ordered line integral methods in 3D: simplex search

Ordered line integral methods in 3D: hierarchical update strategy The idea: ◮ Line updates are incident in adjacent triangle updates ◮ Triangle updates incident in adjacent tetrahedron updates ◮ For h small enough, strictly convex minimization problems ◮ So: do line updates first, check KKT conditions for adjacent triangle updates, ditto adjacent tetrahedron updates

Ordered line integral methods in 3D: complexity For N × N × N grid with N 3 grid points: ◮ Heap sort: O ( N 3 log N ) ◮ Line updates: O ( f ( N ) N 3 ) ◮ Triangle updates: O ( N 3 ) ◮ Tetrahedron updates: O ( N 3 ) The function f ( N ) is not O (1)—more likely, f ( N ) = O ( N ) or O ( N 2 ) with a very small constant

Numerical results: overview Five examples to test accuracy ◮ Three examples using linear SDEs ◮ Two nonlinear examples with analytic U A Two examples by Tao (2018)—two point attractors with a hyperbolic periodic orbit as a transition state A genetic switch model with asymptotically stable equilibria

Numerical results: errors and runtimes Extensive numerical tests in our JCP paper, brief synposis here Example #1 #2 #3 31 . 4 K 1 . 79 8 . 57 K 1 . 80 93 . 8 K 1 . 78 T vs K ( N = 513) [s] 0 . 181 h 0 . 750 0 . 0808 h 0 . 658 12 . 3 h 1 . 49 E ∞ vs h ( K opt.) 27 . 4 N 4 . 05 61 . 0 N 4 . 11 131 . 0 N 3 . 98 T vs N ( K opt.) [ns] Example #4 #5 66 . 3 K 1 . 80 1 . 94 K 1 . 80 T vs K ( N = 513) [s] 0 . 635 h 0 . 989 19 . 4 h 1 . 54 E ∞ vs h ( K opt.) 49 . 2 N 4 . 10 92 . 6 N 4 . 00 T vs N ( K opt.) [ns]

Numerical results: example E

Numerical results: a system with hyperbolic periodic orbits

Numerical results: a system with hyperbolic periodic orbits

Numerical results: genetic switch model

Numerical results: genetic switch model

Conclusion ◮ Extended the 2D OLIM to 3D, preserving the improved accuracy compared with the OUM ◮ Developed a hierarchical update strategy which reduces runtime in 3D on large meshes from days to hours ◮ Presented guidelines for choosing the neighborhood truncation parameter K

Future work ◮ Compute the quasipotential on (or around) the lower dimensional manifold occupied by the dynamics ◮ Optimal time complexity: O ( N ) ◮ Can the solver be made higher-order accurate?

References Shuo Yang, Samuel F. Potter, and Maria K. Cameron. “Computing the quasipotential for nongradient SDEs in 3D.” Journal of Computational Physics 379 (2019): 325–350. Samuel F. Potter and Maria K. Cameron. “Ordered line integral methods for solving the eikonal equation.” Journal of Scientific Computing (2019): 1–41. Cameron, Maria, and Shuo Yang. ”Computing the quasipotential for highly dissipative and chaotic SDEs an application to stochastic Lorenz’63.” Communications in Applied Mathematics and Computational Science 14.2 (2019): 207-246.