Ordered Line Integral Methods for Solving the Eikonal Equation - - PowerPoint PPT Presentation

Ordered Line Integral Methods for Solving the Eikonal Equation

Samuel F. Potter Maria K. Cameron July 17, 2019

The eikonal equation

Solve ∇u(x) = s(x) for u, where s is the slowness function (c = 1/s is speed of sound/light) Where does it come from? Reduced wave (Helmholtz) equation: ∆ + ω2

c2 w = 0

Ansatz (WKB): w(x) ∼ A(x)eiωu(x) → yields eikonal (and transport equation)

1/25

Ordered line integral methods (OLIMs)

Use local characteristic information to skip updates (important in 3D) Higher accuracy through use of midpoint quadrature rules and increased directional coverage Semi-Lagrangian = ⇒ parametrize local characteristics More tricks that work for the eikonal equation but not the quasipotential concurrent: quasipotential OLIM in 3D [YPC ‘19] OLIMs for computing the quasipotential [DC ‘18]

Related work

2/25

History of related algorithms

eikonal solvers label-setting algorithms Dijkstra-like Dial-like iterative solvers Bou´ e & Dupuis [‘99] fast sweeping method [Zhao] Tsitsiklis [‘95] Sethian [‘96] OLIMs for eikonal [PC ‘19]

3/25

Semi-Lagrangian updates

∂Ω x y ϕ Ω u(y) = minx,ϕ

• u(x) +
• ϕ s
• Approximate and minimize

Approximate using simple quadrature rules over straight line segments Minimize using linear algebra + geometry or numerical optimization 4/25

Stencils, neighborhoods, other methods

FMM (2D) Tsitsiklis (2D) Tsitsiklis (3D) FMM (3D)

5/25

A generic Dijkstra-like eikonal solver

Until done: → Remove minimum node from priority queue → Set node state to valid → Set far node neighbors’ states to trial → Update neighbors

• ur focus

far, trial, valid Grid of N d nodes (d = 2, 3) Compute numerical eikonal U Boundary data: U ≡ 0 on a subset of grid points States: Use a priority queue to sort trial nodes by U value trial

6/25

Contributions

Approximate with simple quadrature rules: mp0, mp1, rhr Update using efficient algorithms: top-down and bottom-up Exact semi-Lagrangian solution for mp0/rhr in all dimensions Numerical results + code profiling Theorem relating mp0 and mp1

7/25 Skip redundant updates

OLIM Classification

top-down bottom-up mp0 mp1 rhr mp0 mp1 rhr

6 18 26 6 18 26 6 18 26 26 26 26

• l

i m 6 m p

• l

i m 1 8 m p

• l

i m 2 6 m p

• l

i m 6 m p 1

• l

i m 1 8 m p 1

• l

i m 2 6 m p 1

• l

i m 6 r h r

• l

i m 1 8 r h r

• l

i m 2 6 r h r

• l

i m 3 d m p

• l

i m 3 d m p 1

• l

i m 3 d r h r

eikonal OLIMs

update alg.

• quad. rule
• num. neighb.

OLIM name 8/25

OLIM Classification

top-down bottom-up mp0 mp1 rhr mp0 mp1 rhr

6 18 26 6 18 26 6 18 26 26 26 26

• l

i m 6 m p

• l

i m 1 8 m p

• l

i m 2 6 m p

• l

i m 6 m p 1

• l

i m 1 8 m p 1

• l

i m 2 6 m p 1

• l

i m 6 r h r

• l

i m 1 8 r h r

• l

i m 2 6 r h r

• l

i m 3 d m p

• l

i m 3 d m p 1

• l

i m 3 d r h r

FMM 3D Tsitsiklis 3D

eikonal OLIMs

update alg.

• quad. rule
• num. neighb.

OLIM name 8/25

OLIM Classification

top-down bottom-up mp0 mp1 rhr mp0 mp1 rhr

6 18 26 6 18 26 6 18 26 26 26 26

• l

i m 6 m p

• l

i m 1 8 m p

• l

i m 2 6 m p

• l

i m 6 m p 1

• l

i m 1 8 m p 1

• l

i m 2 6 m p 1

• l

i m 6 r h r

• l

i m 1 8 r h r

• l

i m 2 6 r h r

• l

i m 3 d m p

• l

i m 3 d m p 1

• l

i m 3 d r h r

eikonal OLIMs

update alg.

• quad. rule
• num. neighb.

OLIM name best 8/25

Quadrature Rules: Setup

U(ˆ x) = min0≤λ≤1

• (1 − λ)U(x0) + λU(x1) + Q

L

0 s(ϕ)dl

• Uλ

ˆ U

ˆ x, ˆ s, ˆ U x0, s0, U0 x1, s1, U1 xλ, sλ, Uλ ϕ Approximate integral with quadrature rule Q:

rhr mp0 mp1 ˆ sˆ x − xλ ˆ

s+sλ 2

• ˆ

x − xλ

• ˆ

s+ s0+s1

2

2

• ˆ

x − xλ

approximate U linearly 9/25

Quadrature Rules: Cost Functions

Frhr(λ) = uλ + ˆ sˆ x − xλ Fmp0(λ) = uλ +

• ˆ

s+

1 d+1

d

i=0 si

2

• ˆ

x − xλ Fmp1(λ) = uλ + ˆ

s+sλ 2

• ˆ

x − xλ U(ˆ x) = minλ∈∆d

• Uλ + Q

L

0 s(ϕ)dl

• ˆ

U = minλ∈∆d F(λ)

With ∆d = {λ ∈ Rd : λi ≥ 0 for i = 1, . . . , d and d

i=1 λi ≤ 1}:

=: F(λ) For each quadrature rule (rhr, mp0, mp1):

10/25

Quadrature Rules: θ-rules

F0(λ) = uλ +

• (1 − θ)ˆ

s +

θ d+1

d

i=0 si

• ˆ

x − xλ F1(λ) = uλ + [(1 − θ)ˆ s + θsλ] ˆ x − xλ Frhr = F0, θ = 0 Fmp0 = F0, θ = 1/2 Fmp1 = F1, θ = 1/2

λ s0 sλ

s0+s1 2

s1 sθ ˆ s θ = 0 θ = 1/2

= ⇒ Cost functions Frhr, Fmp0, Fmp1 generalized by: Specifically:

11/25

E.g., setup for tetrahedron update

ˆ p

p1 p2 p0 δp1 δp2 pλ = p0 + δPλ sθ sθ

λ

ˆ s = s(ˆ p) s1 s2 s0 sλ = s(pλ) 12/25 si = s(pi) δP = p1 − p0 p2 − p0

• =

δp1 δp2

Exact solution for F0 (Frhr and Fmp0)

Theorem: Let δP = p1 − p0 · · · pd − p0

• , and let

QR = qr(δP, 0) and λ∗ = argminλ∈RdF0(λ). Then: pλ∗ =

• p⊤

0 (I − QQ⊤)p0

1 −

• R−⊤ δU

sθh

• 2

λ∗ = −R−1

• Q⊤p0 + pλ∗R−⊤ δU

sθh

• ˆ

U = F(λ∗)

Works in any dimension for any simplex, parametrizes local characteristic in addition to computing ˆ

• U. Proof: geometry & linear

algebra. 13/25

Problem with continuity for Fmp0

ˆ p p2 p2 p1 p0

14/25

Problem with continuity for Fmp0

ˆ p p2 p2 p1 p0

14/25

Set ˆ U = Fmp1(λ∗

mp0)

Validating mp0

Proof: Kantorovich’s theorem + 3 lemmas for verification of conditions. Theorem: Let h be small enough that F1 is strictly convex, let λ∗

0 and λ∗ 1 be (unconstrained) minimizers of θ-rules F0 and

• F1. Then |F1(λ∗

1) − F1(λ∗ 0)| = O(h3).

Compute λ∗

0 using QR

decomposition, then set ˆ U ← Fmp1(λ∗

0) (i.e., use

Fmp1 to evaluate).

∆2 15/25

Solution for mp1

Nothing fancy: − → For tri updates, use hybrid method − → For tetra updates, use sequential quadratic programming

∆2 Our implementation is carefully tested and reasonably optimized 16/25

A few other things...

Some relevant things in the paper which I won’t go into here Causality (of individual updates) for F0 and F1 Convexity of θ-rules F0 and F1 General formula for the update gap of a simplex in d dimensions (for Dial-like algorithm) Simple idea: ∇2F1 = ∇2F0 + “small indefinite perturbation” Hence, perturbation = O(h2) lets proofs for F0 apply to F1 for h small enough

17/25

Local Factoring

At point sources and rarefaction fans, order of convergence can degrade from O(h) to O(h log 1

h) — see also Qi’s talk

Very easy to modify our cost functions to work with the additively factored eikonal equation

point source

rarefaction fan

• bstacle

Qi, Dongping, and Alexander Vladimirsky. ”Corner cases, singularities, and dynamic factoring.” Journal of Scientific Computing 79.3 (2019): 1456-1476. 18/25

Enumerating Updates (top-down)

ˆ p Manually enumerate update tetrahedra and triangles = olim6 neighbors + + + = olim18 neighbors = olim26 neighbors Consider only one octant at a time (by symmetry) Can then enumerate small number of tetrahedra (hence triangle) updates 19/25

Enumerating Updates (top-down)

ˆ p I II III IVa IVb V VIa VIb ˆ p 19/25

Enumerating Updates (top-down)

ˆ p I II III IVa IVb V VIa VIb ˆ p

• lim6

19/25

Enumerating Updates (top-down)

ˆ p I II III IVa IVb V VIa VIb ˆ p

• lim18

19/25

Enumerating Updates (top-down)

ˆ p I II III IVa IVb V VIa VIb ˆ p

• lim26

19/25

Skipping Updates (top-down)

Case 1: F0 (i.e. Frhr and mp0) = ⇒ exact solution via QR Case 2: F1 (i.e. Fmp1) = ⇒ constrained optimization

Skip all incident lower-dimensional updates (i.e., after a tetra update, skip all incident tri and line updates (6 total))

∆2 2 2 2 1 1 1

λ1 λ2

skip

λ∗ p0

p1

p2 If λ∗

0 /

∈ ∆d, check its location For d = 2, easy to divide R2\∆2 into disjoint “skip zones” Skip zone tells us which lower-dimensional updates to skip 20/25

Enumerating Updates (bottom-up)

p0 = pnew

ˆ p p2 p1

Check if tri update minimizer λ∗ is minimizer of for adjacent tetrahedron updates (using KKT conditions) λ ∈ ∆d is a simple linear matrix inequality Check if tri update argmin is argmin

• f tetra update by computing only
• ne Lagrange multiplier

Start with p0 = pnew, do line update Check gradient to see if line update “minimizer” is minimizer of adjacent triangle updates (olim3d) 21/25

Enumerating Updates (bottom-up)

p0 = pnew

ˆ p p2 p1

Check if tri update minimizer λ∗ is minimizer of for adjacent tetrahedron updates (using KKT conditions) λ ∈ ∆d is a simple linear matrix inequality Check if tri update argmin is argmin

• f tetra update by computing only
• ne Lagrange multiplier

Start with p0 = pnew, do line update Check gradient to see if line update “minimizer” is minimizer of adjacent triangle updates (olim3d) 21/25 minimizer of triangle update?

Enumerating Updates (bottom-up)

p0 = pnew

ˆ p p2 p1

Check if tri update minimizer λ∗ is minimizer of for adjacent tetrahedron updates (using KKT conditions) λ ∈ ∆d is a simple linear matrix inequality Check if tri update argmin is argmin

• f tetra update by computing only
• ne Lagrange multiplier

Start with p0 = pnew, do line update Check gradient to see if line update “minimizer” is minimizer of adjacent triangle updates (olim3d) 21/25 minimizer of tetra update?

Enumerating Updates (bottom-up)

p0 = pnew

ˆ p p2 p1

Check if tri update minimizer λ∗ is minimizer of for adjacent tetrahedron updates (using KKT conditions) λ ∈ ∆d is a simple linear matrix inequality Check if tri update argmin is argmin

• f tetra update by computing only
• ne Lagrange multiplier

Start with p0 = pnew, do line update Check gradient to see if line update “minimizer” is minimizer of adjacent triangle updates (olim3d) 21/25 minimizer!

Enumerating Updates (bottom-up)

p0 = pnew

ˆ p p2 p1

Check if tri update minimizer λ∗ is minimizer of for adjacent tetrahedron updates (using KKT conditions) λ ∈ ∆d is a simple linear matrix inequality Check if tri update argmin is argmin

• f tetra update by computing only
• ne Lagrange multiplier

Start with p0 = pnew, do line update Check gradient to see if line update “minimizer” is minimizer of adjacent triangle updates Only try higher-dim updates if λ∗ in interior = ⇒ only one constraint active (olim3d) 21/25

Test Problem (Single Point Source)

u(x) = 1

2x⊤A1/2x

s(x) = xA = √ x⊤Ax A is spd 22/25

Test Problem (Single Point Source)

u(x) = 1

2x⊤A1/2x

s(x) = xA = √ x⊤Ax A is spd 22/25

Profiling Using Valgrind

Too little “heap” time to measure when running olim26 and olim3d FMM and olim6 rhr spend smallest amount of time updating = ⇒ can detect time spent in heap Very sensitive to how state is stored (valid, trial, far) Related: “Fast Methods for Eikonal Equations: an Experimental Survey”, G´

• mez, ´

Alvarez, Garrido, Moreno. 2015. 23/25

Summary

OLIMs are fast and accurate, especially in 3D

• lim3d mp0 best overall: rhr speed + mp1 accuracy

Improved directional coverage leads to improved accuracy—midpoint quadrature rule important here Bottom-up algorithm generally faster and a bit more accurate, but top-down may be useful for multiple arrivals. Data structures just as important as algorithms, but little emphasis on them in the literature

24/25

Future Work

Dahiya, Daisy, and Maria Cameron. “Ordered line integral methods for computing the quasi-potential.” Journal of Scientific Computing (2018). Potter, Samuel F., and Maria K. Cameron. “Ordered Line Integral Methods for Solving the Eikonal Equation.” arXiv:1902.06825 (2019). In revision. Yang, Shuo, Samuel F. Potter, and Maria K. Cameron. “Computing the quasipotential for nongradient SDEs in 3D.” Journal of Computational Physics 379 (2019). more details Higher order semi-Lagrangian Dijkstra-like algorithm Use two-level grid or quad-tree to store front states Detect multiple arrivals (caustics...) — HF acoustics application Use update gap to allow some parallelization (at least SIMD...)

References

memory savings 25/25