A parallel algorithm for computing the flow complex - Friedrich-Schiller ...

A parallel algorithm for computing the flow complex Joachim Giesen

Lars K¨ uhne

Fakult¨ at f¨ ur Mathematik und Informatik Friedrich-Schiller Universit¨ at, Jena Germany [email protected] [email protected]

December 14, 2012

The point set

Let P ⊂ Rd be the point set in general position, with at least d + 1 points.

The distance function

x N(x)

I

distance function h : Rd 3 x 7→ minp∈P kx − pk

I

Neighbors of x N(x) = {p ∈ P|h(x) = kx − pk}

Gradient of the distance function

∂h

x

N(x) = d(x)

I

driver d(x) is the center of the smallest enclosing ball of N(x)

I

gradient of h(x), ∂h (x) =

x−d(x) h(x)

The flow

∂h x d(x)

I

flow φ : [0, ∞) × Rd → Rd

I

defined by φ(0, x) = x, and lim

t↓t0

φ(t, x) − φ(t0 , x) = ∂h (φ(t0 , x)) t − t0

The Flow line

d(x)

I

x

∂h

orbit/flow-line of a point x: φ(x) = {φ(t, x)|t ≥ 0}

Critical point

0 = ∂h x = d(x)

I

x ∈ Rd : ∂h (x) = 0

I

x ∈ conv (N(x))

I

index of a critical point is the dimension of the affine hull of N(x)

I

Minima with i(x) = 0, and Maxima with i(x) = d

I

the remaining critical points are saddle points

Flow complex

Flow complex

Flow complex

d X

(−1)i · ni = 1

i=0 I

ni is the number of critical points of index i

Stable and unstable manifolds

I

stable manifold of x contains the points, that flow into x S(x) = {y ∈ Rd | lim φ(t, y ) = x} t↓∞

I

given a neighborhood U around x V (U) = {y ∈ Rd |∃z ∈ U, t ≥ 0 : φ(t, z) = y }

I

, unstable manifold U(x) =

\ NeighborhoodUofx

V (U)

The Hasse diagram

I

for two critical points x and y , S(y ) is incident to S(x), if U(y ) ∩ S(x) 6= ∅, i.e. there is a point in U(y ) that flows into x

The algorithm

The algorithm

Two basic operations 1. Ascent: increasing h(x) along the gradient 2. Descent: decreasing h(x) along the gradient

The primitives Nearest neighbor along a ray: Let v be the ray-vector, p ∈ V and q ∈ P \ V . Solve k(x + t · v ) − qk2 = k(x + t · v ) − pk2 , and pick the minimum tq > 0, along with q Projection onto an affine hull: In order to check for x ∈ CONV (N(x)), we compute the coefficients λi for which X X c= λi · pi , pi ∈ V , λi = 1 i

Here c is the orthogonal-projection of x onto V . We use a dynamic QR-decomposition to solve the least-squares problem of computing the projection of x onto V .

Example

d(x) x

Example

x

Example

x

Example

x

Example

Bounding simplex

I

low-dimensional(sub-)structures serialize the exploration

I

a bounding simplex allows for a scalable parallel exploration

Experimental results Critical point statistics data set random within sphere random within sphere curve spiral-3d spiral-3d spiral-3d beethoven-3d-model cup-3d-model

dimension 3 4 3 3 4 5 3 3

distribution 50-137-122-34 50-199-284-178-44 50-73-56-32 30-64-61-26 30-41-12-0-0 30-41-12-0-0-0 500-1076-649-72 300-766-599-132

time (s) 0.03 2.66 0.7 0.03 0.97 182.4 11.81 0.92

Experimental results Scalability

I

On the left: 200 random points in 4 dimensions

I

On the right: 100 random points in 5 dimensions

Conclusions

I

a new algorithm to compute the flow complex, without explicit computation of the Delaunay triangulation

I

numerical robustness

I

it is inherently parallelizable, and scales well with the number of cores in a shared-memory architecture

I

future work: high-level scalability by subdivision

Thank you!

Questions?