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

1/29 Setting System evolving according to SDE (Itˆ o): dx = b ( x ) dt + √ ǫdw, x ∈ R 3 standard Brownian motion C 1 vector field small noise parameter For attractor A ⊆ R 3 , compute quasipotential: � L 0 ( � b ( ψ ) �� ψ ′ � − b ( ψ ) · ψ ′ ) ds U ( x ) = min ψ,L geometric action (obtained from Freidlin-Wentzell action) Only makes sense to compute U if b is nongradient

2/29 Computing the Quasipotential on a Grid attractor A U | A ≡ 0 goal : march U forward causally in semi-Lagrangian fashion (approximate locally minimum action paths)

3/29 Dijkstra-like Solvers for HJ Equations Computing quasipotential involves solving Hamilton-Jacobi equation of the form: � � ∇ u ( x ) �∇ u ( x ) � = 1 F x, �∇ u ( x ) � Previous approach: ordered upwind method (OUM)—requires bounded anisotropy (Sethian & Vladimirsky ‘01): Υ = F max /F min May be unbounded —OUM for quasipotential handles this (Cameron ‘12). Additional error is O ( h 2 ) OLIMs (Dahiya and Cameron ‘17) reduce error constant by 10 2 – 10 3 by higher-order quadrature rules

4/29 Description of the Algorithm For grid points x , four states: unknown considered accepted attractor A accepted front

4/29 Description of the Algorithm Update considered neighbors of new accepted front point unknown considered Kh accepted attractor A accepted front

4/29 Description of the Algorithm Update considered neighbors of new accepted front point unknown considered accepted attractor A accepted front

5/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 612 3 was taking days. . .

6/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 overall number of updates significantly

7/29 Simplex updates In 3D, three types of updates: line , triangle , tetrahedron E.g., for the tetrahedron update: Q 3 ( x 0 , x 1 , x 2 , x ) = min λ { U λ + Q ( x λ , x ) } , where U λ = U ( x 0 ) + λ 1 ( U ( x 1 ) − U ( x 0 )) + λ ( U ( x 2 ) − U ( 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 i + x � b mi = b , i = 0 , 1 , 2 , 2 subject to λ 1 ≥ 0 , λ 2 ≥ 0 , λ 1 + λ 2 ≤ 1 .

8/29 Simplex update minimization problem Cost function for simplex update: f ( λ ) = U λ + � x − x λ �� b mλ � − ( x − x λ ) · b mλ ∇ f = δU + � x − x λ � � b mλ � � b mλ � B ⊤ b mλ + � x − x λ � X ⊤ ( x − x λ ) − B ⊤ ( x − x λ ) − X ⊤ b mλ � ⊤ � ⊤ � ⊤ + � X ⊤ ( x − x λ ) B ⊤ b mλ X ⊤ ( x − x λ ) B ⊤ b mλ � � H = � x − x λ �� b mλ � + � x − x λ � � x − x λ � X ⊤ X − � x − x λ � � b mλ � � b mλ � B ⊤ B + � b mλ � 3 B ⊤ b mλ B ⊤ b mλ � � � b mλ � � ⊤ − B ⊤ X − � ⊤ � x − x λ � X ⊤ ( x − x λ ) X ⊤ ( x − x λ ) B ⊤ X � � − Where: � b ( x m 1 ) − b ( x m 0 ) b ( x m 2 ) − b ( x m 0 ) � B = � x 0 − x 1 � x 0 − x 2 X = δU = ( U ( x 1 ) − U ( x 0 ) , U ( x 2 ) − U ( x 0 ))

8/29 Simplex update minimization problem Cost function for simplex update: f ( λ ) = U λ + � x − x λ �� b mλ � − ( x − x λ ) · b mλ ∇ f = δU + � x − x λ � � b mλ � away from equilibria of b ( x ) , H is � b mλ � B ⊤ b mλ + � x − x λ � X ⊤ ( x − x λ ) positive definite for h small enough if − B ⊤ ( x − x λ ) − X ⊤ b mλ K ∼ h − α , where 0 < α < 1 � ⊤ � ⊤ � ⊤ + � X ⊤ ( x − x λ ) B ⊤ b mλ X ⊤ ( x − x λ ) B ⊤ b mλ � � H H = � x − x λ �� b mλ � + � x − x λ � � x − x λ � X ⊤ X − � x − x λ � � b mλ � � b mλ � B ⊤ B + � b mλ � 3 B ⊤ b mλ B ⊤ b mλ � � � b mλ � � ⊤ − B ⊤ X − � ⊤ � x − x λ � X ⊤ ( x − x λ ) X ⊤ ( x − x λ ) B ⊤ X � � − Where: � b ( x m 1 ) − b ( x m 0 ) b ( x m 2 ) − b ( x m 0 ) � B = � x 0 − x 1 � x 0 − x 2 X = δU = ( U ( x 1 ) − U ( x 0 ) , U ( x 2 ) − U ( x 0 ))

9/29 A fast search for admissible updates Neighborhoods used to organize different types of updates: x ∈ Z 3 : � x − x 0 � 1 = 1 � � N 1 ( x 0 ) = x ∈ Z 3 : � x − x 0 � 1 = 2 and � x − x 0 � ∞ = 1 � � N 2 ( x 0 ) = x ∈ Z 3 : � x − x 0 � 1 = 3 and � x − x 0 � ∞ = 1 � � N 3 ( x 0 ) = N ( x 0 ) = N 1 ( x 0 ) ∪ N 2 ( x 0 ) ∪ N 3 ( x 0 ) x ∈ Z 3 : � x − x 0 � 2 ≤ K N K � � far ( x 0 ) ≈

9/29 A fast search for admissible updates Neighborhoods used to organize different types of updates: x ∈ Z 3 : � x − x 0 � 1 = 1 � � N 1 ( x 0 ) = x ∈ Z 3 : � x − x 0 � 1 = 2 and � x − x 0 � ∞ = 1 � � N 2 ( x 0 ) = x ∈ Z 3 : � x − x 0 � 1 = 3 and � x − x 0 � ∞ = 1 � � N 3 ( x 0 ) = N ( x 0 ) = N 1 ( x 0 ) ∪ N 2 ( x 0 ) ∪ N 3 ( x 0 ) x ∈ Z 3 : � x − x 0 � 2 ≤ K N K � � far ( x 0 ) ≈ considered point being updated

10/29 A fast search for admissible updates Motivation for using these simplexes: interpolation error E.g., assume U ( x ) = U 0 + g ⊤ x + 1 2 x ⊤ Hx In 1D, on interval of size h , max interpolation error = 1 8 h 2 H 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

11/29 Skipping updates using the KKT conds. Constrained minimization problem in canonical form: f ( λ ) = U λ + � b mλ �� x − x λ � − b mλ · ( x − x λ ) minimize λ 1 ≥ 0 subject to λ 2 ≥ 0 1 − λ 1 − λ 2 ≥ 0 Lagrangian : L ( λ, µ ) = f ( λ ) − µ 1 λ 1 − µ 2 λ 2 − µ 3 (1 − λ 1 − λ 2 ) KKT conditions : � 1 � � 0 � � − 1 � ∇ λ L = ∇ f − µ 1 − µ 2 − µ 3 = 0 , 0 1 − 1 µ 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: f ( λ ) = U λ + � b mλ �� x − x λ � − b mλ · ( x − x λ ) minimize λ 1 ≥ 0 subject to λ 2 ≥ 0 to verify, ends up requiring only simple comparisons involving ∇ f 1 − λ 1 − λ 2 ≥ 0 Lagrangian : L ( λ, µ ) = f ( λ ) − µ 1 λ 1 − µ 2 λ 2 − µ 3 (1 − λ 1 − λ 2 ) KKT conditions : � 1 � � 0 � � − 1 � ∇ λ L = ∇ f − µ 1 − µ 2 − µ 3 = 0 , 0 1 − 1 µ 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 .

12/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 x x 0 Important part of getting a speed-up in 3D and verified x 2 numerically x λ ∗ x 1

13/29 The hierarchical update strategy A brute force update algorithm: 1. x 0 ← considered point with smallest U value 2. x 0 becomes accepted front 3. accepted front points in N ( x 0 ) with no considered neighbors become accepted 4. for all considered x ∈ N K far ( x 0 ) , set U ( x ) to minimum Q 3 ( x, x 0 , x 1 , x 2 ) value where x 1 , x 2 are accepted front points in N 1 ( x 0 ) ∪ N 2 ( x 0 ) 5. for all unknown x ∈ N ( x 0 ) (a) x becomes considered (b) for all accepted front y 0 ∈ N K far ( x ) , set U ( x ) to minimum Q 3 ( y 0 , y 1 , y 2 , x ) value where y 1 , y 2 are accepted front points in N 1 ( y 0 ) ∪ N 2 ( y 0 )