Stochastic Heat Kernel Estimation on Sampled Manifolds Symposium on - - PowerPoint PPT Presentation

Stochastic Heat Kernel Estimation on Sampled Manifolds

Symposium on Geometry Processing 2017

Tristan Aumentado-Armstrong & Kaleem Siddiqi July 4, 2017

Centre for Intelligent Machines & School of Computer Science McGill University, Canada

Table of contents

• 1. Introduction
• 2. Theoretical Background
• 3. Algorithm
• 4. Empirical Results
• 5. Discussion

1

Introduction

Sampled Manifolds and Heat Kernels

Many datasets can be viewed as random samples from a Riemannian manifold. (e.g. 3D point clouds, natural images). The underlying manifold has an associated heat kernel, with many applications in computer vision and graphics:

• Feature extraction (Sun et al, 2009; Gebal et al, 2009)
• Shape matching (Ovsjanikov et al, 2010)
• Shape retrieval (Bronstein et al, 2011)

Introduction 2

Heat Kernel Applications for Shapes I

Feature point detection & description (Sun et al, 2009) Invariance properties (Bronstein et al, 2010.)

Introduction 3

Heat Kernel Applications for Shapes II

Symmetry (Sun et al, 2009; Ovsjanikov et al, 2010) Segmentation (Gebal et al, 2009)

Introduction 4

Heat Kernel Applications for Shapes III

Point matching (Ovsjanikov et al, 2010)

Introduction 5

Heat Kernel Applications for Shapes IV

Shape retrieval & global description (Ovsjanikov et al, 2009)

Introduction 6

Computing the Heat Kernel

The most straightforward spectral approach involves a global eigenvalue problem. Others have devised more sophisticated approaches:

• Multiresolution prolongation (Vaxman et al, 2010)
• Rational Chebyshev approximation (Patan´

e and Spagnuolo, 2014)

Introduction 7

A Stochastic Perspective

The heat kernel on a manifold has been studied with stochastic methods since Ito (1950). Basic idea: the heat kernel describes the transition density function of Brownian motion on the manifold. Question Can we use this outlook to compute the heat kernel? What properties would such an algorithm have?

Introduction 8

Answer

We show how to use a Monte Carlo algorithm, via trajectories simulated on the sampled manifold.

Introduction 9

Theoretical Background

The Laplace-Beltrami Operator

Let (M, g) be a Riemannian manifold in RD with dim(M) = d. The Laplace-Beltrami operator (LBO) on it is given by: ∆g = divg gradg = ∇g · ∇g In local coordinates, this can also be written as: ∆g =

• j
• k

gjk ∂2 ∂xj∂xk

• Diffusion

−

• ℓ

gjkΓℓ

jk

∂ ∂xℓ

• Convection

where Γℓ

jk are the Christoffel symbols. Theory 10

The Riemannian Heat Equation

The heat kernel is the fundamental solution to the heat equation on the manifold: ∆gu = ∂ ∂tu In spectral form: Kt(x, y) =

∞

• i=0

exp(−λit)φi(x)φi(y) where ∆gφi = −λiφi. The heat kernel Kt(x, y) is our quantity of interest.

Theory 11

Stochastic Calculus on Manifolds (I)

An Ito Diffusion is a stochastic differential equation (SDE): dXt = µ(Xt) dt

• Drift

+ σ(Xt) dBt

• Diffusion

The drift (convection) and diffusion coefs define the SDE. Changing µ Changing σ

Theory 12

Stochastic Calculus on Manifolds (II)

Every Ito diffusion has an infinitesimal generator: L =

• ℓ

µℓ ∂ ∂xℓ + 1 2

• i
• j

[σ σT]ij ∂2 ∂xi∂xj that encodes its short-time behaviour. Definition: Brownian motion (BM) on a manifold BM on (M, g) is an Ito Diffusion with L = ∆g/2. For BM, the infinitesimal generator can be rewritten as: L = 1 2∆g = 1 2

• j
• k

gjk

• [σ σT ]jk

∂2 ∂xj∂xk +

• ℓ

−1 2 gjkΓℓ

jk

• µℓ

∂ ∂xℓ

Theory 13

Stochastic Calculus on Manifolds (III)

Hence, given local coordinates on a manifold, we can simulate Brownian motion on it via: dXi

t = µi(Xt) dt +

• k

σi

k(Xt) dBk t

with convection (drift) µ and diffusion σ coefficients given by: µk = −1 2

• i
• j

gijΓk

ij

σ = √g via the metric tensor g.

Theory 14

Relation to the Manifold Heat Kernel

The BM Xt is a Markov process with a transition density function p(x, t|y). There is an intuitive connection between p(x, t|y) & Kt(x, y):

• p(x, t|y): probability density of a random walk (BM) on

(M, g) reaching x from y in time t.

• Kt(x, y): amount of some substance that diffuses over

(M, g) from x to y in time t.

• Prop. (Hsu, 2002): Equivalence of p(x, t|y) & Kt(x, y)

p(x, t|y) = Kt(x, y) is the heat kernel of ∂tu = Lu.

Theory 15

Transition Density Estimation (I)

We may estimate p(t, x|y) instead of Kt(x, y). But how? A Monte Carlo approach. Use kernel density estimation (KDE) over a set

• f nT sample trajectories.

x y High p(x, t|y) x y Low p(x, t|y)

Theory 16

Transition Density Estimation (II)

Kernel density estimation (KDE) performs distance weighting. p(x, t|y) = 1 nT

nT

• j=1

Kδ

• X(j)

x,t − y

• X(j)

x,t: trajectory j at time t starting from x

• Kδ: a kernel function (e.g. Gaussian density)
• δ: the kernel bandwidth (controls the smoothing level)

Includes some theoretical guarantees. Problem: assumes process is in Euclidean space, not on a curved manifold

Theory 17

Transition Density Estimation (III)

KDE on manifolds is possible, but requires geodesic distances. However, we can approximate the manifold KDE.

• Prop. (Ozakin & Gray, 2009): Manifold KDE Approx.
• Kt(x, y) =

1 nTδd

h nT

• j=1

Ψ

• ||X(j)

x,t − y||

δh

• with X(j)

x,t trajectory j at time t starting from x, Ψ a kernel,

|| · || the ambient norm, and δh the bandwidth. Intuition: at small scales, (M, g) has intrinsic distances like RD, but has surface area like Rd.

Theory 18

Algorithm

Outline of Algorithm

The algorithm has three steps:

• 1. Local surface construction
• 2. Stochastic trajectory generation
• 3. Transition density estimation

Algorithm 19

• 1. Moving Least Squares Surface Construction

Use PCA on neighbours Nk(p) to get local coordinates per point p. Fit zp with weighted moving least squares (MLS), where: zp(x, y) =

• i,j

γi,j xiyj with deg(zp) = 2. Local surface: Λp = {x, y, zp(x, y)}

Global Coordinates Local Coordinates

Algorithm 20

• 2. Stochastic Trajectory Generation (I)

Can use standard stochastic numerical integration on a Λp: Xt = Xs + t

s

µ(Xr) dr + t

s

σ(Xr) dBr Problem: how to transition from Λp to Λq? Repeat until Tmax:

• Simulate for time Tδ on Λp
• Find closest point q to Xt
• Project Xt to Λq & p ← q

Algorithm 21

• 2. Stochastic Trajectory Generation (II)

Problem: evaluating µ and σ is expensive. Assumption: minor geometric effects on Xt are less important than the variance between trajectories at high t. Solution: if t > τ, linearize the curvature terms in the SDE. Similar to the multiresolution idea (Vaxman et al, 2010).

Algorithm 22

• 3. Kernel Density Estimation (I)

Parameters:

• P: input points (|P| = N)
• S: set of feature points
• nT: trajectories per point
• HT: times of interest

Approach:

• Generate nT trajectories from each q ∈ S
• ∀ p ∈ P, q ∈ S, t ∈ HT:
• Kt(p, q) =

1 nTδd

h

√ πd

nT

• j=1

exp

• −||X(j)

q,t − p||2

δ2

h

• which is manifold KDE with a Gaussian kernel.

This fills Kt(x, y) as an |S| × N × |HT| array.

Algorithm 23

• 3. Kernel Density Estimation (II)

Problem: how to choose the kernel bandwidth δh? (a) Theory (Milstein et al, 2004) shows optimal δh ≈ n

−1 (4+d)

T

. (b) Use the short-time autodiffusion expansion:

• Kt(p, p) ≈

1 (4πt)d/2

• 1 + R(p)t

6

• where R(p) is the Ricci curvature (for d = 2, R = 2K).

Calibrate δh with Kt(p, p) around the theory-based value using a few small times.

Algorithm 24

Algorithmic Complexity

• 1. Moving least squares surface construction

O(N log(N))

• 2. Stochastic trajectory generation

O (|S| nT Tmax[log(N)/Tδ + Cstep/h])

• 3. Kernel density estimation

O(nT |S| N |HT|)

Assumes: using a KD-tree in 3D

Algorithm 25

Empirical Results

The Sphere (I)

We sample N = 500 points from S2. Manifold fidelity: t ||Xt|| Tδ Green: ||Xt|| = 1. Blue: ||Xt|| = 1.005. The true Kt can be computed from the spherical harmonics.

Experiments 26

The Sphere (II)

Kt(x, x) Kt(x, y)

(Antipodal)

t t t nT=70 nT=150 nT=1000

Green: stochastic estimate; Blue: analytic value.

Experiments 27

TOSCA Meshes (I)

Only the surface construction step (1) is changed by having a mesh structure (i.e. choosing the neighbours Nk(p)). Manifold Fidelity:

(Problems if high curvature & Tδ)

Experiments 28

TOSCA Meshes (II)

Compare to linear FEM spectral approach (Reuter et al, 2006):

Stoch. Spec. Stochastic approach has sparser log arrival-time maps (nT = 500).

Experiments 29

Timing on Larger Point Clouds (I)

Using KD-trees. (No edge/face information used). Fix: Tδ = 0.25, hs = 0.05, nT = 200, Tmax = 500, |HT| = 5.

Experiments 30

Timing on Larger Point Clouds (II)

For densely sampled manifolds, interpolate the MLS surfaces (4) and estimated heat kernel (3) to neighbours.

Heat kernels at 3 points for t = 25

Experiments 31

Timing on Larger Point Clouds (III)

Filling Kt(x, y) as an N × |S| × 5 array Model N MLS (s) |S| STG (s) KDE (s) Chinese Lion 152807 30.5 1 159.31 38.75 5 811.26 192.93 Armadillo 172974 36.75 1 174.02 43.24 5 874.21 214.85 Buddha 719560 143.44 1 176.11 173.63 5 896.88 870.8 Linear dependence on N (log(N) for STG) and |S|.

Experiments 32

Discussion

Summary

We’ve devised a three-step, stochastic Monte Carlo algorithm for estimating the heat kernel on sampled manifolds, and tested it on shapes in 3D.

Discussion 33

Pros & Cons

Advantages

• Good asymptotics
• Stable
• Generalizable (e.g. other

metrics, higher dimensions)

• Parallelizable
• Time-accuracy tradeoff

parameters

• Some theory guarantees

Disadvantages

• Currently slower than
• ther spectrum-free

methods

• Hard to determine

convergence

• Time-accuracy tradeoff

parameters

• HKS-only extraction

unclear

Discussion 34

Future Work

• Algorithmic
• Manifold fidelity
• Speed
• Theoretical guarantees
• Error bounds
• Convergence rates
• Extracting other signatures (analogous to Feynman-Kac

formula)

Discussion 35

Acknowledgements

Thanks to my supervisor Prof. Siddiqi and the rest of the lab. Thanks to Natural Sciences and Engineering Research Council

• f Canada (NSERC) for funding.

36

Thank you for listening!

36

First-Order Approximation to Kt(x, y) = p(x, t|y)

Suppose x and y are close and t is small. Then: p(x, t|y) ≈ PGaussian(x; y + µ(y)t, σ(y)2t) where µ, σ are from the SDE for BM. In geodesic normal coordinates at y, this is: p(x, t|y) ≈ (2πt)−d/2 exp −||x − y||2

2

2t

• This is equivalent to the heat kernel discretization used by

Belkin & Niyogi (2003) for non-linear dimensionality reduction.

Transition Density Estimation Theory

We may estimate p(t, x|y) instead of Kt(x, y). But how? Let Xt approximate BM, with integration step-size h.

• Prop. (Hu et al., 1996): Convergence of p(x, t|y)

p(x, t|y) = lim

h→0 E [φh(Xx,t − y)]

where φh(a) = PGaussian(a; 0, h2) and Xx,t starts at x. This suggests a Monte Carlo approach is possible. In practice, we can use kernel density estimation (KDE)

• ver a set of nT sample trajectories instead.

The Heat Equation Describes a BM Ensemble

Intuitively, the heat equation describes the average movement

• f an ensemble of BMs. Formally:

u(t, a) = E[f(Xt)] is a solution to the diffusion equation ∂tu = Lu; u(0, a) = f(a) where X0 = a and L = ∆g/2.