Delaunay Core Bifiltration

Odin Hoff Gardå, https://odinhg.github.io

In this notebook, we demonstrate how one can use the Delaunay core bifiltration in multipers to compute persistent homology of a noisy point cloud.

import multipers as mp
import matplotlib.pyplot as plt
import numpy as np
from multipers.data import three_annulus
from multipers.filtrations import CoreDelaunay

np.random.seed(0)

[KeOps] Warning : CUDA libraries not found or could not be loaded; Switching to CPU only.

Background

The core bifiltration was introduced in the paper “Core Bifiltration” and is a generalization of the degree-Čech bifiltration (setting \(\beta=1/2\)). The core bifiltration is interleaved with the well-known multicover bifiltration and is a parameter-free construction that is robust to noise (in the sense of stability with respect to the Prokhorov metric). Intersecting with Voronoi cells in the usual way, one obtains the Delaunay core bifiltration which is a bifiltration of the Delaunay complex. The Delaunay core bifiltration is interleaved with the core bifiltration and also enjoys similar stability properties.

Consider a finite point cloud \(A\) in a metric space. Let \(B^\beta_{r,k}(a)\) be the closed metric ball \(B(a,r)\) whenever \(r/\beta\geq d_k^A(a)\) and the empty set otherwise. Here, \(d_k^A(a)\) denotes the distance from \(a\) to its \(k\)-th nearest neighbor in \(A\). Now, let \(W^\beta_{r,k}(a)\) denote \(B^\beta_{r,k}(a)\) intersected with the Voronoi cell of \(a\). The Delaunay core bifiltration is then defined in filtration degree \((r,-k)\) as the union of \(W^\beta_{r,k}(a)\) for all \(a\in A\). Taking the nerve of the covering \(\{W^\beta_{r,k}(a)\}\) yields a bifiltered simplicial complex, which is the one implemented in multipers.

Example Usage

Point Cloud Sample

We start by sampling \(2000\) points from three annuli in the plane with \(2000\) outliers added.

num_pts = 2000
num_outliers = 2000

X = three_annulus(num_pts, num_outliers)

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
ax.scatter(*X.T, s=5, alpha=0.7)
ax.set_aspect(1)
plt.show()

Constructing the Delaunay Core Bifiltration

The (Delaunay) core bifiltration takes as input a point cloud in \(\mathbb{R}^d\) and a parameter \(\beta\in\mathbb{R}_{\geq 0}\) balancing the effect of the radius of balls and \(k\)-nearest neighbor distances. The keyword argument ks specifies which values of \(k\) to consider when constructing the bifiltration.

k_max = len(X)/10
num_steps = 100
ks = np.unique(np.linspace(1, k_max, num=num_steps, dtype=int))
beta = 0.5

st = CoreDelaunay(points=X, beta=beta, ks=ks)

Computing and Visualizing Persistent Homology

Using module_approximation from multipers, we compute persistence and visualize it.

Notice the [1,0] direction, with the max_error set to the ks precision. This is motivated by the fact that the filtration values are all of the form \((x,k)\) for some \(x\in\mathbb{R}\) and \(k\in \mathrm{ks}\). These parameters thus ensure that the line slices drawn by MMA hit each filtration value of the filtration, leading to a much more efficient computation.

pers = mp.module_approximation(
    st,
    direction=[1,0], 
    max_error=.5*k_max/num_steps
)
box = mp.grids.compute_bounding_box(pers)
pers.plot(degree=1, alpha=0.9, xlabel="radius", ylabel="k", min_persistence=10e-3, box=box)
plt.show()

For completeness’ sake, we now try some different values of the input parameters and include \(H_0\) in the plot.

k_max = len(X)//10
num_steps = 200
ks = np.unique(np.linspace(1, k_max, num=num_steps, dtype=int))
beta = 1.5

st = CoreDelaunay(points=X, beta=beta, ks=ks)

pers = mp.module_approximation(
    st,
    direction=[1,0],
    max_error=.5*k_max/num_steps,
)
box = mp.grids.compute_bounding_box(pers)
pers.plot(alpha=0.9, dpi=400, xlabel="radius", ylabel="k", min_persistence=10e-3, box=box)
plt.tight_layout()
plt.show()