Topological Data Analysis: Finding Clusters, Loops, and Voids Across Scales

Clustering output depends on a scale hyperparameter (k, eps, dendrogram cut height). Topological data analysis (TDA) indexes every scale at once via a filtration and records which features persist across many thresholds.

1 Clustering scale parameters

Common clustering methods require one scale before analysis:

K-means: number of groups k.
DBSCAN: neighbourhood radius eps.
Hierarchical clustering: cut height on the merge tree.

Wrong scale collapses real structure or splits one group into noise.

TDA replaces a single threshold with a nested sequence of shapes indexed by radius $\varepsilon \geq 0$. Features with long lifetimes are treated as signal; features with birth $\approx$ death are treated as noise.

2 Vietoris-Rips complex

A point cloud has no intrinsic shape until pairs are joined by distance.

Simplices and simplicial complex:

0-simplex: point.
1-simplex: edge (2 points).
2-simplex: filled triangle (3 points).
3-simplex: filled tetrahedron (4 points).

Given $X = \{x_1, \dots, x_n\}$ and $\varepsilon \geq 0$, the Vietoris-Rips complex $R_\varepsilon(X)$ includes a simplex on every subset whose pairwise distances are $\leq \varepsilon$:

\[ \{x_{i_0}, \dots, x_{i_k}\} \in R_\varepsilon(X) \iff \|x_{i_p} - x_{i_q}\| \leq \varepsilon \text{ for every pair } p, q. \]

The rule uses only the pairwise distance matrix.

As $\varepsilon$ increases, simplices only accumulate: $R_0(X) \subseteq R_{\varepsilon_1}(X) \subseteq R_{\varepsilon_2}(X) \subseteq \cdots$. This nested family is the filtration.

3 Filtration on four points

Four corners of a unit square illustrate components and a loop at distinct scales. Neighbouring corners are distance 1 apart; opposite corners are $\sqrt{2} \approx 1.414$ apart.

Code

import numpy as np
import matplotlib.pyplot as plt
from itertools import combinations

square = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
side, diagonal = 1.0, np.sqrt(2)


def plot_rips(points, eps, ax):
    n = len(points)
    dists = np.linalg.norm(points[:, None] - points[None, :], axis=-1)
    for i, j, k in combinations(range(n), 3):
        if dists[i, j] <= eps and dists[j, k] <= eps and dists[i, k] <= eps:
            ax.add_patch(plt.Polygon(points[[i, j, k]], alpha=0.3, color="tab:blue", zorder=1))
    for i, j in combinations(range(n), 2):
        if dists[i, j] <= eps:
            ax.plot(*zip(points[i], points[j]), color="tab:blue", zorder=2)
    ax.scatter(points[:, 0], points[:, 1], s=80, zorder=3, color="black")
    ax.set_title(f"$\\varepsilon$ = {eps:.3f}")
    ax.set_xlim(-0.4, 1.4)
    ax.set_ylim(-0.4, 1.4)
    ax.set_aspect("equal")
    ax.set_xticks([])
    ax.set_yticks([])


fig, axes = plt.subplots(1, 4, figsize=(13, 3.5))
for ax, eps in zip(axes, [0.5, side, 1.2, diagonal]):
    plot_rips(square, eps, ax)
plt.tight_layout()
plt.show()

Panel sequence:

$\varepsilon = 0.5$: four isolated points ($\beta_0 = 4$).
$\varepsilon = 1$: boundary edges form a hollow square ($\beta_0 = 1$, $\beta_1 = 1$).
$\varepsilon = 1.2$: loop unchanged.
$\varepsilon = \sqrt{2}$: diagonals connect; triangles fill the interior; loop dies.

3.1 Betti numbers

Betti number $\beta_k$ counts $k$-cycles modulo $(k{+}1)$-boundaries:

\[ \beta_k = \dim \ker \partial_k - \dim \operatorname{im} \partial_{k+1}. \]

Each term is a matrix rank or nullity.

Code

def rips_simplices(points, eps, max_dim):
    n = len(points)
    dists = np.linalg.norm(points[:, None] - points[None, :], axis=-1)
    simplices = {0: [(i,) for i in range(n)]}
    for dim in range(1, max_dim + 1):
        simplices[dim] = [
            combo
            for combo in combinations(range(n), dim + 1)
            if all(dists[i, j] <= eps for i, j in combinations(combo, 2))
        ]
    return simplices


def boundary_matrix(simplices_k, simplices_km1):
    index = {s: i for i, s in enumerate(simplices_km1)}
    d = np.zeros((len(simplices_km1), len(simplices_k)))
    for col, simplex in enumerate(simplices_k):
        for i in range(len(simplex)):
            face = simplex[:i] + simplex[i + 1 :]  # drop vertex i -> a face
            d[index[face], col] = (-1) ** i
    return d


def betti_numbers(points, eps, max_dim=2):
    simplices = rips_simplices(points, eps, max_dim)
    ranks = {0: 0}
    for dim in range(1, max_dim + 1):
        if simplices[dim] and simplices[dim - 1]:
            ranks[dim] = np.linalg.matrix_rank(boundary_matrix(simplices[dim], simplices[dim - 1]))
        else:
            ranks[dim] = 0
    betti = {}
    for dim in range(max_dim):
        cycles = len(simplices[dim]) - ranks[dim]      # dim ker(d_dim)
        boundaries = ranks.get(dim + 1, 0)              # dim im(d_{dim+1})
        betti[dim] = cycles - boundaries
    return betti


for eps in [0.5, side, 1.2, diagonal]:
    print(f"eps = {eps:.3f}:  betti_0 = {betti_numbers(square, eps)[0]},  betti_1 = {betti_numbers(square, eps)[1]}")

eps = 0.500:  betti_0 = 4,  betti_1 = 0
eps = 1.000:  betti_0 = 1,  betti_1 = 1
eps = 1.200:  betti_0 = 1,  betti_1 = 1
eps = 1.414:  betti_0 = 1,  betti_1 = 0

ripser reports each feature as a (birth, death) pair instead of a per-scale Betti count.

Code

from ripser import ripser

dgms = ripser(square, maxdim=1)["dgms"]
print("H0 (birth, death):\n", dgms[0])
print("H1 (birth, death):\n", dgms[1])

H0 (birth, death):
 [[ 0.  1.]
 [ 0.  1.]
 [ 0.  1.]
 [ 0. inf]]
H1 (birth, death):
 [[1.         1.41421354]]

One $H_1$ interval [1.0, 1.41421354]: born at $\varepsilon = 1$, dead at $\varepsilon = \sqrt{2}$. Three finite $H_0$ bars die at $\varepsilon = 1$; one component persists ($\infty$).

4 Homology dimensions

$\beta_0$ equals connected components: threshold the distance graph at $\varepsilon$ and count components. That is single-linkage clustering.

Higher dimensions use the same quotient:

\[ \beta_k = \underbrace{\text{(all $k$-dimensional cycles)}}_{\text{the ``things''}} \Big/ \underbrace{\text{(the ones that bound a $(k{+}1)$-dimensional patch)}}_{\text{``trivial'' cycles}} \]

$k=0$: points modulo paths → clusters.
$k=1$: loops modulo filled disks → holes.
$k=2$: shells modulo filled solids → voids.

Persistent homology runs this at every $k$ and every $\varepsilon$ simultaneously.

5 Noisy ring vs blobs

Data provenance (synthetic):

make_circles (scikit-learn): points on two concentric circles + Gaussian noise.
make_blobs: three separated Gaussian clouds.
Ground truth is known: one dataset has a loop; one does not.
Real analogues: periodic trajectories (ring) vs mixture populations (blobs).

Code

from sklearn.datasets import make_circles, make_blobs
from persim import plot_diagrams

rng = np.random.default_rng(42)

ring, _ = make_circles(n_samples=200, noise=0.05, factor=0.6, random_state=42)
blobs, _ = make_blobs(n_samples=200, centers=3, cluster_std=0.4, random_state=42)

fig, axes = plt.subplots(2, 2, figsize=(10, 8))

for ax, data, title in zip(axes[0], [ring, blobs], ["Noisy ring", "Three blobs"]):
    ax.scatter(data[:, 0], data[:, 1], s=10)
    ax.set_title(title)
    ax.set_aspect("equal")

for ax, data, title in zip(
    axes[1], [ring, blobs], ["Ring: persistence diagram", "Blobs: persistence diagram"]
):
    dgms = ripser(data, maxdim=1)["dgms"]
    plot_diagrams(dgms, ax=ax, show=False)
    ax.set_title(title)

plt.tight_layout()
plt.show()

Persistence diagram reading:

Each point: (birth, death) of one feature.
Points lie above the diagonal; vertical distance = lifetime.
Near-diagonal points: short-lived noise.
Far from diagonal: persistent structure.

Ring: one $H_1$ point far above the diagonal (the loop). Blobs: $H_1$ near the diagonal; $H_0$ shows three long-lived components before merge.

5.1 Single-linkage equivalence

$H_0$ death times equal single-linkage merge heights on the same data.

Code

from scipy.cluster.hierarchy import linkage

dgm0 = ripser(blobs, maxdim=0)["dgms"][0]
death_times = np.sort(dgm0[np.isfinite(dgm0[:, 1]), 1])

Z = linkage(blobs, method="single")
merge_heights = np.sort(Z[:, 2])

print("H0 death times match single-linkage merge heights:", np.allclose(death_times, merge_heights))

H0 death times match single-linkage merge heights: True

6 TDA vs clustering

Aspect	Traditional clustering	TDA
Scale	One hyperparameter (`k`, `eps`, cut height)	All scales in one filtration
Output	Single partition	Persistence diagram (birth/death per feature)
Shape assumptions	Often convex blobs or one density scale	Components, loops, voids without convexity
Noise handling	Depends on hyperparameter	Short lifetime ≈ noise
Sensitivity	Different `k`/`eps` → different partitions	Deterministic given distance + filtration
Cost	$O(n)$–$O(n^2)$	$O(n^3)$ or worse for full Rips; subsample at scale

Clustering answers “what is the grouping at this scale?” TDA answers “at which scales does structure exist, and with what lifetime?” Loops and voids are not expressible as point partitions.

7 Diagram vectorization

A persistence diagram is a variable-size multiset of points. Models need fixed-length vectors.

Vectorization families:

Summary statistics: max lifetime, total persistence, count above threshold, persistence entropy.
Persistence images: rasterize (birth, lifetime) into a fixed grid, weighted by persistence.
Persistence landscapes: functional summary in a vector space; supports averaging and kernel methods.

Applications: material microstructures, molecular conformations, time series via Takens delay embedding.

7.1 Ring vs blob classifier

Data provenance (synthetic):

60 Gaussian blobs and 60 jittered unit circles, 100 points each.
Labels known by construction.
Task: can $H_1$ summary statistics separate the two shape classes?

Each cloud → three numbers from its $H_1$ diagram; random forest for classification.

Code

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split


def sample_blob(n, rng):
    return rng.normal(size=(n, 2))


def sample_ring(n, rng, radius=1.0, noise=0.05):
    theta = rng.uniform(0, 2 * np.pi, n)
    r = radius + rng.normal(0, noise, n)
    return np.column_stack([r * np.cos(theta), r * np.sin(theta)])


def persistence_features(point_cloud):
    """Vectorize a point cloud's H1 persistence diagram into 3 summary numbers."""
    h1 = ripser(point_cloud, maxdim=1)["dgms"][1]
    if len(h1) == 0:
        return np.zeros(3)
    lifetimes = h1[:, 1] - h1[:, 0]
    return np.array([lifetimes.max(), lifetimes.sum(), (lifetimes > 0.1).sum()])


rng = np.random.default_rng(0)
X, y = [], []
for _ in range(60):
    X.append(persistence_features(sample_blob(100, rng)))
    y.append(0)  # "blob" shape class
    X.append(persistence_features(sample_ring(100, rng)))
    y.append(1)  # "ring" shape class

X, y = np.array(X), np.array(y)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)

clf = RandomForestClassifier(n_estimators=100, random_state=0)
clf.fit(X_train, y_train)
print(f"Test accuracy: {clf.score(X_test, y_test):.2f}")

Test accuracy: 1.00

Perfect hold-out accuracy on this toy split demonstrates the pipeline (diagram → features → classifier), not a general accuracy claim.

8 Practical constraints

Cost: full Rips persistence is roughly $O(n^3)$. Subsample, use landmark/witness complexes, or cap maximum simplex dimension beyond a few thousand points.
Metric choice: TDA removes the scale hyperparameter, not the metric (Euclidean vs cosine vs learned distance).
Complement to clustering: use $H_0$ persistence to read cluster count and merge confidence; then run k-means or DBSCAN at an informed scale.
Mapper: cover-based graph summary (separate from persistence diagrams).
Libraries: ripser/persim for fast low-dimensional Rips; GUDHI and giotto-tda for broader complexes and scikit-learn transformers.