Persistent Homology

Overview & Motivation

Consider two point clouds with identical first- and second-order statistics — same mean, same variance, same covariance. One is sampled from an annulus (a ring), the other from a Gaussian cluster. Every classical statistical test says they’re the same. But one has a hole and the other doesn’t, and that difference can be the signal you need — whether you’re detecting anomalies in sensor networks, classifying dynamical regimes in time series, or finding structure in high-dimensional embeddings.

The problem is: what does “has a hole” even mean for a finite point cloud? Points are discrete; they don’t have holes in any obvious sense. The Vietoris-Rips construction (see Simplicial Complexes) gives us a way to build a continuous-ish structure from discrete data, but it introduces a new problem: the result depends on a scale parameter $\varepsilon$. Choose $\varepsilon$ too small and you see only isolated points. Too large and everything collapses into a single blob.

Persistent homology is the answer: instead of choosing a single scale, track the topology at every scale simultaneously. Features that persist across a wide range of $\varepsilon$ are real structure; features that flicker in and out are noise. This insight — that persistence is a proxy for significance — is what makes TDA useful for machine learning.

Not all complexes are created equal: the Vietoris-Rips complex can introduce phantom cycles — topological artifacts that do not correspond to real features of the data. The Čech complex avoids this via the Nerve Theorem.

---

Formal Framework

Filtrations

For the Vietoris-Rips complex, the filtration is parameterized by $\varepsilon$: as $\varepsilon$ increases from 0 to $\infty$, simplices enter the complex at their birth time (the minimum $\varepsilon$ at which they appear). An edge $[v_i, v_j]$ is born at $\varepsilon = d(v_i, v_j)$. A triangle $[v_i, v_j, v_k]$ is born at $\varepsilon = \max(d(v_i, v_j), d(v_i, v_k), d(v_j, v_k))$.

Chain Groups and the Boundary Operator

To detect “holes,” we need algebra. The key idea: a loop is a collection of edges that forms a cycle, and a hole is a loop that isn’t the boundary of something filled in.

The following property is the engine of homology:

Homology Groups

The $\partial \circ \partial = 0$ property gives us two important subgroups:

Homology counts holes that are not filled in:

• $\beta_0$ = number of connected components
• $\beta_1$ = number of loops (independent 1-cycles not bounding a filled region)
• $\beta_2$ = number of voids (enclosed cavities)

Persistent Homology

Now we combine filtrations with homology. As $\varepsilon$ increases, new simplices enter the complex. Each new simplex either:

1. Creates a topological feature (a new cycle appears that isn’t a boundary), or
2. Destroys a topological feature (a cycle that was a hole gets filled in)

Points far from the diagonal have high persistence — they represent features that survive across a wide range of scales. Points near the diagonal have low persistence — they are topological “noise,” features that appear and immediately disappear. The diagonal points are a technical convenience that makes the space of persistence diagrams well-behaved.

---

The Matrix Reduction Algorithm

Persistent homology is defined algebraically, but it is computed via matrix reduction — a variant of Gaussian elimination applied to the boundary matrix. Understanding this algorithm demystifies what ripser and gudhi are actually doing.

The key operation is column reduction: for each column $j$ (processed left to right), find the lowest nonzero entry (the “pivot” or “low” of column $j$). If another column $i < j$ has the same pivot row, add column $i$ to column $j$ (over $\mathbb{Z}_2$, addition is XOR). Repeat until no two columns share the same pivot.

Matrix Reduction Animator

Interactive visualization coming soon — step through the column reduction algorithm and watch birth-death pairs emerge.

---

Visual Intuition

The visualization below is the core interactive experience. The top panel shows a Vietoris-Rips complex built from points sampled near a circle. The bottom panels show the persistence diagram and barcode.

Drag the $\varepsilon$ slider and watch all three panels update together:

ε = 0.30

Persistence Diagram

Persistence Barcode

What to notice:

1. $H_0$ features (teal): At $\varepsilon = 0$, every point is its own connected component — 15 features are born simultaneously. As $\varepsilon$ increases, components merge. Each merge event kills a component: you see teal points appear in the persistence diagram with increasing death values. One component survives to $\infty$ (the essential $H_0$ feature).

2. $H_1$ features (purple): At a critical $\varepsilon$ around 0.7-0.9, the edges form a cycle around the point cloud but the interior hasn’t filled in yet — a loop is born. As $\varepsilon$ increases further, 2-simplices fill the interior, and the loop dies. The persistence (death minus birth) of this $H_1$ feature is the key topological signal: it tells you “this point cloud has a circular structure.”

3. Signal vs. noise: The long-lived $H_1$ bar stands out from the short bars near the diagonal. That gap between persistent and ephemeral features is what makes TDA work in practice — it’s a built-in notion of statistical significance, grounded in geometry rather than distributional assumptions.

---

The Stability Theorem

The practical utility of persistent homology rests on a single theorem: small perturbations of the input data produce small changes in the persistence diagram. Without this, the method would be useless — any noise would scramble the output.

For point cloud data, the relevant corollary is:

In plain language: if you jiggle each point by at most $\delta$, the persistence diagram moves by at most $2\delta$. This is what makes persistent homology safe to use on noisy real-world data.

---

Extended Persistence

Ordinary persistent homology tracks features born and killed in a sublevel set filtration: as $\varepsilon$ increases from 0 to $\infty$, simplices enter and topological features appear and disappear. But some features are essential — they are born and never die. The most important essential feature is the single surviving connected component ($H_0$ with death $= \infty$), but for spaces with nontrivial higher homology, essential features carry real topological information that ordinary persistence cannot capture.

Extended persistence (Cohen-Steiner, Edelsbrunner & Harer, 2009) solves this by appending a superlevel set filtration to the sublevel set filtration, creating a complete picture where every feature has a finite death time.

Extended Persistence Diagram

Interactive visualization coming soon — explore ordinary, relative, and extended persistence pairs on a torus.

---

Persistence Modules

The algebraic structure underlying persistent homology is the persistence module — a one-parameter family of vector spaces connected by linear maps. Understanding this algebraic framework clarifies why barcodes exist, why stability holds, and why multi-parameter persistence is fundamentally harder.

The persistence module arising from a filtration $K_0 \subseteq K_1 \subseteq \cdots$ is the family of homology groups $V_t = H_k(K_t)$, with the maps $v_{s,t}$ induced by the inclusions $K_s \hookrightarrow K_t$. The maps tell us which features at scale $s$ “survive to” scale $t$.

This is the theorem that makes barcodes well-defined. Without it, the barcode would be an arbitrary visualization choice; with it, the barcode is the unique decomposition of the persistence module into indecomposable summands. The proof relies on the classification of finitely generated modules over a principal ideal domain — the same algebraic machinery that gives the Jordan normal form for matrices.

For the metric structure on the space of persistence modules — including the interleaving distance, the Isometry Theorem, and the stability framework — see the Barcodes & Bottleneck Distance page.

---

Working Code

Computing Persistence with Ripser

import numpy as np
from ripser import ripser
from persim import plot_diagrams

# Generate an annulus (has a hole) and a Gaussian cluster (no hole)
np.random.seed(42)
n = 200

theta = np.random.uniform(0, 2 * np.pi, n)
r = 1.0 + np.random.normal(0, 0.1, n)
annulus = np.column_stack([r * np.cos(theta), r * np.sin(theta)])

cluster = np.random.randn(n, 2) * np.std(annulus, axis=0) + np.mean(annulus, axis=0)

# Compute persistence (H₀ and H₁)
result_annulus = ripser(annulus, maxdim=1)
result_cluster = ripser(cluster, maxdim=1)

# The annulus has a high-persistence H₁ feature (the loop)
dgm1 = result_annulus['dgms'][1]
max_pers = np.max(dgm1[:, 1] - dgm1[:, 0])
print(f"Annulus max H₁ persistence: {max_pers:.4f}")

# The cluster has only low-persistence H₁ noise
dgm1_c = result_cluster['dgms'][1]
max_pers_c = np.max(dgm1_c[:, 1] - dgm1_c[:, 0])
print(f"Cluster max H₁ persistence: {max_pers_c:.4f}")

Verifying the Stability Theorem Empirically

from persim import bottleneck

# Clean circle
n_pts = 100
theta_base = np.linspace(0, 2 * np.pi, n_pts, endpoint=False)
base_cloud = np.column_stack([np.cos(theta_base), np.sin(theta_base)])
base_dgm = ripser(base_cloud, maxdim=1)['dgms'][1]

# Add increasing noise and measure diagram movement
for noise_std in [0.0, 0.1, 0.2, 0.3, 0.5]:
    noise = noise_std * np.random.randn(n_pts, 2)
    noisy_cloud = base_cloud + noise
    noisy_dgm = ripser(noisy_cloud, maxdim=1)['dgms'][1]

    bn = bottleneck(base_dgm, noisy_dgm)
    linf = np.max(np.abs(noise)) if noise_std > 0 else 0
    print(f"noise σ={noise_std:.1f}: ||perturbation||∞={linf:.3f}, d_B={bn:.3f}, "
          f"bound holds: {bn <= linf + 1e-10 or noise_std == 0}")

Vectorization: Persistence Images for ML

Persistence diagrams aren’t fixed-size vectors, so we can’t feed them directly to most classifiers. Persistence images (Adams et al., 2017) solve this by placing a Gaussian at each point $(b, d-b)$ and evaluating on a grid:

from persim import PersistenceImager

# Create a persistence image from H₁
pimgr = PersistenceImager(pixel_size=0.05, birth_range=(0, 1.5), pers_range=(0, 1.0))
pimgr.fit([result_annulus['dgms'][1]])

img = pimgr.transform(result_annulus['dgms'][1])
print(f"Feature vector shape: {img.flatten().shape}")
# This fixed-size vector goes into any sklearn classifier

Full Pipeline: Time Series Classification

The Takens delay embedding theorem connects time series analysis to topology:

$x(t) \mapsto \bigl(x(t),\; x(t+\tau),\; x(t+2\tau),\; \ldots,\; x(t+(d-1)\tau)\bigr)$

Periodic signals produce loops in the embedding; chaotic signals do not.

1. 1D time series x(t)
2. Takens delay embedding dim d, delay τ
3. Point cloud in ℝᵈ
4. Persistent homology Ripser
5. Feature vector summary stats of persistence
6. Classifier Random Forest, SVM, …
7. Regime label

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score

def takens_embedding(ts, dim=2, delay=10):
    """Takens delay embedding of a 1D time series."""
    n = len(ts) - (dim - 1) * delay
    return np.array([ts[i:i + dim * delay:delay] for i in range(n)])

def topological_features(ts):
    """Extract topological feature vector from a time series."""
    embedded = takens_embedding(ts, dim=2, delay=10)
    result = ripser(embedded, maxdim=1)
    features = []
    for d in range(2):  # H₀ and H₁
        dgm = result['dgms'][d]
        finite = dgm[np.isfinite(dgm[:, 1])]
        if len(finite) > 0:
            pers = finite[:, 1] - finite[:, 0]
            features.extend([
                np.max(pers),        # Max persistence
                np.mean(pers),       # Mean persistence
                np.std(pers),        # Std of persistence
                np.sum(pers > 0.1),  # Count of significant features
                np.sum(pers ** 2),   # Total persistence (L² norm)
            ])
        else:
            features.extend([0, 0, 0, 0, 0])
    return np.array(features)

# Generate labeled time series from 3 dynamical regimes
X, y = [], []
for regime in range(3):
    for _ in range(50):
        t = np.linspace(0, 6 * np.pi, 300)
        noise = 0.05 * np.random.randn(300)
        if regime == 0:    # Periodic
            ts = np.sin(t) + noise
        elif regime == 1:  # Quasi-periodic
            ts = np.sin(t) + 0.5 * np.sin(np.sqrt(2) * t) + noise
        else:              # Chaotic (logistic map)
            ts = np.zeros(300)
            ts[0] = 0.1
            for i in range(1, 300):
                ts[i] = 3.9 * ts[i-1] * (1 - ts[i-1])
        X.append(topological_features(ts))
        y.append(regime)

X, y = np.array(X), np.array(y)
scores = cross_val_score(RandomForestClassifier(n_estimators=100, random_state=42),
                         X, y, cv=5)
print(f"5-fold CV accuracy: {scores.mean():.3f} ± {scores.std():.3f}")

This pipeline — time series → Takens embedding → persistent homology → feature vector → classifier — is one of the most successful applications of TDA in practice. It works because the persistence features are invariant to reparameterization of the time series (you can stretch or compress time and the topology doesn’t change) and stable under noise (by the Stability Theorem).

Matrix Reduction from Scratch

The standard algorithm for computing persistence, implemented directly from Theorem 3:

import numpy as np

def reduce_boundary_matrix(D):
    """Column-reduce a boundary matrix over Z₂.

    Parameters
    ----------
    D : ndarray of shape (m, m) — the filtered boundary matrix (Z₂ entries)

    Returns
    -------
    R : ndarray — the reduced matrix
    pairs : list of (creator_idx, destroyer_idx) or (creator_idx, None)
    """
    R = D.copy() % 2
    m = R.shape[1]
    low = {}     # row_index -> col_index of column that has this row as its pivot
    pairs = []

    for j in range(m):
        # Find the lowest 1 in column j
        nonzero = np.where(R[:, j] == 1)[0]
        while len(nonzero) > 0 and nonzero[-1] in low:
            # Another column i has the same low — add it (Z₂: XOR)
            i = low[nonzero[-1]]
            R[:, j] = (R[:, j] + R[:, i]) % 2
            nonzero = np.where(R[:, j] == 1)[0]

        if len(nonzero) > 0:
            pivot = nonzero[-1]
            low[pivot] = j
            pairs.append((pivot, j))  # pivot creates, j destroys
        # else: column j is zero — creates an essential feature

    # Identify essential features (creators not paired as destroyers)
    all_destroyers = set(j for _, j in pairs)
    all_creators = set(i for i, _ in pairs)
    for j in range(m):
        if j not in all_destroyers and j not in all_creators:
            pairs.append((j, None))  # essential feature

    return R, pairs

# Example: three points with edges at distances 1, 2, 3
# Filtration order: v0, v1, v2, [v0,v1], [v0,v2], [v1,v2]
D = np.array([
    # v0 v1 v2 e01 e02 e12
    [0, 0, 0, 1,  1,  0],  # v0
    [0, 0, 0, 1,  0,  1],  # v1
    [0, 0, 0, 0,  1,  1],  # v2
    [0, 0, 0, 0,  0,  0],  # e01
    [0, 0, 0, 0,  0,  0],  # e02
    [0, 0, 0, 0,  0,  0],  # e12
])
R, pairs = reduce_boundary_matrix(D)
print("Persistence pairs (creator_idx, destroyer_idx):")
for c, d in pairs:
    print(f"  σ_{c} paired with {'σ_' + str(d) if d is not None else '∞'}")

Extended Persistence with GUDHI

import gudhi
import numpy as np

# Height function on a sampled torus
n = 500
theta = np.random.uniform(0, 2 * np.pi, n)
phi = np.random.uniform(0, 2 * np.pi, n)
R, r = 2.0, 1.0
torus = np.column_stack([
    (R + r * np.cos(phi)) * np.cos(theta),
    (R + r * np.cos(phi)) * np.sin(theta),
    r * np.sin(phi)
])

# Build alpha complex and compute extended persistence
alpha = gudhi.AlphaComplex(points=torus)
stree = alpha.create_simplex_tree()

# Extended persistence of the height function (z-coordinate)
# gudhi computes this via the extended filtration
stree.extend_filtration()
dgms = stree.extended_persistence()

print("Extended persistence pairs by type:")
for i, (name, dgm) in enumerate(zip(
    ['Ordinary', 'Relative', 'Extended+', 'Extended-'], dgms)):
    print(f"  {name}: {len(dgm)} pairs")
    for dim, (b, d) in dgm[:3]:
        print(f"    H_{dim}: [{b:.2f}, {d:.2f})")

---

Connections & Cross-Track Bridges

Persistent homology sits at the center of a growing web of mathematical and applied techniques, connecting deeply to several other areas on formalML.

Within the Topology & TDA track:

• Simplicial Complexes provides the combinatorial substrate: the VR filtration is a nested sequence of simplicial complexes, and the boundary operator defined there is the algebraic engine of homology.
• Barcodes & Bottleneck Distance develops the metric structure on persistence diagrams. The Structure Theorem (Theorem 4 above) guarantees that barcodes are well-defined; the bottleneck and Wasserstein distances turn the space of diagrams into a metric space where stability makes sense.
• Statistical TDA treats persistence diagrams as random objects — statistical estimators of the true homology of an underlying distribution. Bootstrap confidence sets separate topological signal from sampling noise.
• The Mapper Algorithm provides a complementary geometric/visual summary. Where persistent homology produces algebraic invariants (barcodes), Mapper produces an interpretable graph. PH on the Mapper graph quantifies its topological features; multiscale Mapper uses PH to select optimal parameters.

Across tracks:

• Graph Laplacians & Spectrum: The 0th persistent homology of a VR filtration is equivalent to single-linkage clustering — the threshold $\varepsilon$ at which two points first become connected in $H_0$ equals the single-linkage distance. The number of connected components at scale $\varepsilon$ equals the number of zero eigenvalues of the graph Laplacian of the $\varepsilon$-neighborhood graph.
• Measure-Theoretic Probability: The convergence of empirical persistence diagrams to the population diagram relies on measure-theoretic concentration bounds. The convergence rate $O((\log n / n)^{1/d})$ exhibits the curse of dimensionality — the same dimension-dependent rate that appears in density estimation and minimax theory.
• Concentration Inequalities: The stability of persistence diagrams under noise is a Lipschitz condition, and the convergence rate of empirical diagrams involves Hausdorff concentration bounds that parallel McDiarmid-type arguments.

Applications: Persistence features have been used for protein classification (Cang & Wei, 2017), materials science (Hiraoka et al., 2016), cosmology (Xu et al., 2019), neuroscience (Giusti et al., 2015), and financial market analysis. The Takens delay embedding connects persistent homology to time series classification — see the working code above for the full pipeline.

---

Exercises

Foundational

Exercise 1. Consider four points at the corners of a unit square in $\mathbb{R}^2$: $(0,0)$, $(1,0)$, $(1,1)$, $(0,1)$. Write out the filtered boundary matrix for the VR filtration (vertices born at 0, four side-edges born at 1, two diagonal edges born at $\sqrt{2}$). Perform column reduction by hand and read off all birth-death pairs for $H_0$ and $H_1$.

Exercise 2. Verify $\partial_1 \circ \partial_2 = 0$ for the boundary of a tetrahedron $[v_0, v_1, v_2, v_3]$ over $\mathbb{Z}_2$. Write out the $6 \times 4$ matrix $\partial_2$ (rows = edges, columns = triangular faces) and the $4 \times 6$ matrix $\partial_1$ (rows = vertices, columns = edges). Compute the product and verify every entry is 0 mod 2.

Intermediate

Exercise 3. Prove that the longest $H_0$ bar in any VR persistence barcode always has death time $= \infty$. Hint: The last surviving connected component is never merged into anything.

Exercise 4. The total persistence of a persistence diagram is $\sum_i (d_i - b_i)^p$ for some power $p > 0$. Show that for $H_0$ persistence of $n$ points, the total persistence (with $p = 1$) equals the total weight of the minimum spanning tree of the point cloud. Hint: Each $H_0$ death corresponds to an MST edge.

Advanced

Exercise 5. Implement the full Takens embedding pipeline from the Working Code section on a real dataset: download daily closing prices for a stock index (e.g., S&P 500), compute the delay embedding, run persistent homology, and classify 60-day rolling windows into “high-volatility” vs “low-volatility” regimes using persistence features. Report cross-validation accuracy.

Exercise 6. Prove that persistent homology is invariant under bijective reparameterization of the filtration. That is, if $\phi: \mathbb{R} \to \mathbb{R}$ is a monotonically increasing bijection and $K_t$ is a filtration, then the persistence diagram of $K_{\phi(t)}$ is related to that of $K_t$ by the coordinate change $(b, d) \mapsto (\phi^{-1}(b), \phi^{-1}(d))$. Hint: Reparameterization preserves the inclusion order of simplices, hence the matrix reduction is identical.