Topological Data Analysis for Imaging and Machine Learning

Contents

• Introduction and overview
• 0. Review: topology, linear algebra over F2 and union-find
• 1. Simplicial complexes and filtrations
• 2. Simplicial homology
• 3. Smooth and discrete Morse theory
• 4. Persistent homology and persistence diagrams
• 5. Reeb graphs, contour trees and Mapper
• 6. Morse–Smale complexes
• 7. Topological inference in high dimensions
• 8. Topological signatures and features for learning
• 9. Implementations in Python and C++
• 10. Complements for self-study
• 11. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with an M1 in mathematics, computer science or physics who knows linear algebra and some Python, and who wants to measure the shape of data — its components, loops and cavities — in images, scalar fields and point clouds, robustly to noise.

Empirical first: a 64 × 64 image of three Gaussian bumps plus Gaussian noise of standard deviation 0.05 has 721 local maxima. Persistence gives the two secondary bumps the values 0.734 and 0.64, and every noise peak at most 0.172 (Code 4.2). One threshold in that gap recovers the three objects, and the segmentation that follows finds the three regions of the clean image (Code 6.1). In the same way, 5 % of outliers cut the persistence of the loop of a sampled circle from 0.80 to 0.44, while the distance to a measure keeps 0.63 (Code 7.1).

Topological data analysis (TDA) describes data by the features that survive continuous deformation: connected components, loops, cavities, and the way they appear and disappear across scales. The course builds these descriptors from scratch — simplicial complexes, homology, Morse theory, persistence — and uses them for scientific visualization in 2D and 3D and for machine learning in high dimension. It shares tools with Geometric Data Analysis, Geometry Processing and Geometric Deep Learning and Optimal Transport for Machine Learning, whose Wasserstein distances reappear between persistence diagrams.

What you will be able to do

• Build Čech, Vietoris–Rips, alpha and lower-star filtrations from point clouds, images and scalar fields, and bound how far the Rips and Čech complexes differ.
• Compute Betti numbers by Gaussian elimination over the two-element field, and explain when the choice of coefficients matters.
• Classify the critical points of a piecewise-linear scalar field and check the result against the Euler characteristic.
• Compute persistence diagrams with the reduction algorithm and with union-find, and compare them with the bottleneck and Wasserstein distances.
• State and use the stability theorem to separate topological signal from noise and to simplify a scalar field.
• Compute merge trees, contour trees and Mapper graphs, and relate them to the Reeb graph of a function.
• Segment an image with the cells of a Morse–Smale complex, simplified by persistence.
• Infer the homology of a sampled shape with guarantees, and make the inference robust to outliers with the distance to a measure.
• Cluster data with persistence-based clustering (ToMATo) and choose its parameter from a diagram.
• Turn persistence diagrams into stable features (landscapes, images, kernels) for classification, and differentiate a topological loss.
• Choose data structures and algorithms (simplex trees, clearing, cohomology) that make these computations fast, in Python or C++.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com presents topological data analysis [1–6] as a sound family of techniques for the interactive analysis and visualization of data in imaging and machine learning. As acquired or simulated data sets in 2D, 3D and nD grow in size and complexity, these approaches help users understand their data through its topological and geometric structure.

In low dimensions, the description continues, TDA lets users extract, interact with, measure and classify geometric features defined by level sets or integral lines. Simplification based on persistent homology builds multi-scale topological representations that support robust analyses and comparisons despite noise. In higher dimensions, these techniques have been adapted into new clustering algorithms.

The stated purpose is to introduce the main concepts of this recent field and to illustrate their use in imaging (scientific visualization) and machine learning, from a mathematical and a practical point of view, with the Gudhi [7] and TTK [8] libraries. Each item of the official page is matched below.

Beyond the syllabus, this page adds a review chapter (ch. 0), a chapter of complements with the longer proofs (ch. 10), a research chapter (ch. 11), graded exercises, project subjects and a written-exam proposal.

Where this course fits in the programme

Official prerequisites: admission to the MVA; programming skills in Python and/or C++ are recommended (master-mva.com). The instructors listed on master-mva.com are Frédéric Chazal (INRIA) and Julien Tierny (CNRS). The course is taught in period 1, so only the refreshers come before it.

• Refreshers — Linear Algebra, Analysis and Optimization for ranks, kernels and quotients (ch. 0, 2, 4); Probability and Statistics for the sampling and bootstrap arguments of ch. 7; Scientific Python, NumPy and PyTorch for the labs and the differentiable losses of ch. 8.
• Related courses of semester 1 — Geometric Data Analysis (intrinsic dimension, neighbourhood graphs), Geometry Processing and Geometric Deep Learning (meshes and Laplacians), Optimal Transport for Machine Learning (Wasserstein distances), Medical Image Analysis, Computational Statistics (bootstrap) and Introduction to Statistical Learning (the classifiers of ch. 8).
• Courses of semester 2 that use its tools — Point Clouds and 3D Modelling, Graphs in Machine Learning, Kernel Methods for Machine Learning (kernels on diagrams), Geometry and Shape Spaces, Deformable Models and Geodesic Methods (watersheds), Deep Learning for Medical Imaging (topological losses for segmentation) and Algorithms and Learning for Protein Science.
• Afterwards — the Research Internship; ch. 11 lists directions that a first research project can start from.

Course outline

• Review: topology, linear algebra over F2 and union-find — homotopy equivalence, the Hausdorff distance, ranks over the two-element field, disjoint-set forests.
• Simplicial complexes and filtrations — abstract complexes, nerves, Čech, Rips and alpha complexes, triangulated images, filtrations.
• Simplicial homology — chains, boundaries, Betti numbers by rank, the Euler–Poincaré formula, the role of coefficients.
• Smooth and discrete Morse theory — critical points and cell attachments, piecewise-linear critical points, Banchoff’s formula, Forman’s discrete gradients.
• Persistent homology and persistence diagrams — persistence modules, the elder rule, the reduction algorithm, bottleneck distance, stability, simplification of images.
• Reeb graphs, contour trees and Mapper — level-set topology, merge trees by union-find, Mapper as a discrete Reeb graph.
• Morse–Smale complexes — integral lines, ascending and descending manifolds, quadrangle cells, watershed segmentation and its simplification.
• Topological inference in high dimensions — reach and sampling theorems, Rips inference and stability, distance to a measure, ToMATo.
• Topological signatures and features for learning — landscapes, persistence images, kernels, learned layers and differentiable persistence.
• Implementations in Python and C++ — simplex trees, clearing and cohomology, the Gudhi and TTK libraries, a union-find in C++.
• Complements for self-study — the pairing lemma, interleavings and the isometry theorem, collapses, the nerve theorem, statistics of diagrams.
• Research directions and open problems — multiparameter persistence, statistics and learning, large-scale topological visualization.

Session plan

24 h in class = 12 sessions of 2 h.

The volume is the official one of the Université Paris-Saclay programme table (12 h of lectures and 12 h of labs). master-mva.com organises the course in six lectures of two hours, six practical sessions of two hours and a final research seminar of three hours, which the programme table does not count. The plan alternates one lecture and one lab; the content of each session is a proposal. The review chapter is not an official lecture: it is read before session 1.

Level 1 = foundations, 2 = core results and algorithms, 3 = advanced material close to research; the level never decreases from one session to the next.

Self-study, outside the 24 h (about 8.5 h in total): ch. 0 Review, to read before session 1 (≈ 1.5 h) · ch. 9 Implementations in Python and C++, read along the labs, whose two codes are run in Lab 6 (≈ 2 h) · ch. 10 Complements — the pairing lemma, interleavings and the isometry theorem, collapses, the nerve theorem and the statistics of diagrams, each read after the session that states the result (≈ 3 h) · ch. 11 Research directions and open problems, which prepares the final research seminar announced on master-mva.com (≈ 2 h).

How to study this module

• Read ch. 0 before session 1. The Hausdorff distance (Proposition 0.3) and union-find (Algorithm 0.6) return in ch. 4, 5, 6 and 7.
• Follow one thread through the page: the union-find sweep with the elder rule computes connected components (ch. 0), 0-dimensional persistence (ch. 4), merge trees (ch. 5), watershed regions (ch. 6) and ToMATo clusters (ch. 7).
• After each lab, change one parameter of the code — the noise level, the radius, the threshold, the number of neighbours — and predict the new diagram before running it.
• Do the exercises listed in the session plan before the next session; read ch. 10 after the session that states each result. The ★★★ exercises and the paper studies of Projects and assessment prepare the research seminar.
• Python environment: NumPy, SciPy, scikit-learn and PyTorch on CPU; every block runs in a few seconds. The official labs use Gudhi and TTK; the code here rewrites the needed algorithms in plain Python, so that each step is visible and nothing else must be installed.

0. Review: topology, linear algebra over F2 and union-find

Topological data analysis counts features that survive continuous deformation, and it computes them with linear algebra and graph algorithms. This chapter recalls the three M1 tools the course uses everywhere: equivalences of spaces, linear algebra over the two-element field, and the union-find data structure.

0.1 Spaces up to deformation

An annulus deformation retracts onto its middle circle, and a disk onto its centre. A union of small balls around points sampled on a circle is never homeomorphic to the circle, but it can be homotopy equivalent to it. Homotopy equivalence is therefore the right notion for sampled data, and every invariant of this course is a homotopy invariant.

0.2 Distances between sets

Every stability result of ch. 4 and ch. 7 starts from this identity: two samples that are close in Hausdorff distance have close distance functions, hence close offsets.

0.3 Linear algebra over the two-element field

0.4 Union-find

In this course, union-find runs on vertices sorted by a function value. The order in which components merge then carries the information: it gives 0-dimensional persistence (ch. 4), merge trees (ch. 5), watershed regions (ch. 6) and clusters (ch. 7).

1. Simplicial complexes and filtrations

A computer cannot store a circle, but it can store points, edges and triangles. This chapter defines these combinatorial models, builds them from point clouds and images, and nests them along a scale.

1.1 Simplicial complexes

1.2 Nerves

1.3 Čech and Vietoris–Rips complexes

By the nerve theorem, since balls are convex, \check C(X, r) \simeq X^r . The Čech complex models the offset faithfully, but testing whether k + 1 balls intersect needs a smallest enclosing ball; the Rips complex only needs pairwise distances.

The constant of Jung’s theorem is reached: for the vertices of an equilateral triangle of side 1, the Rips triangle appears at r = 1/2 and the Čech triangle at the circumradius 1/\sqrt 3, a ratio of 1.155 (Exercise 2).

import numpy as np
from itertools import combinations

rng = np.random.default_rng(0)
theta = rng.uniform(0, 2 * np.pi, 30)                     # 30 points near the unit circle
P = np.c_[np.cos(theta), np.sin(theta)] + 0.05 * rng.normal(size=(30, 2))
D = np.linalg.norm(P[:, None] - P[None], axis=2)

def meb_radius(a, b, c):
    """Radius of the smallest disk containing three points of the plane."""
    ab, bc, ca = np.linalg.norm(a - b), np.linalg.norm(b - c), np.linalg.norm(c - a)
    s = sorted([ab, bc, ca])
    if s[2] ** 2 >= s[0] ** 2 + s[1] ** 2:                 # right or obtuse: half the longest side
        return s[2] / 2
    u, v = b - a, c - a
    area = abs(u[0] * v[1] - u[1] * v[0]) / 2
    return ab * bc * ca / (4 * area)                       # acute: circumradius

triples = list(combinations(range(30), 3))
meb = np.array([meb_radius(*P[list(t)]) for t in triples])
diam = np.array([max(D[i, j], D[i, k], D[j, k]) for i, j, k in triples])
jung = np.sqrt(4 / 3)                                      # sqrt(2d/(d+1)) for d = 2
for r in [0.3, 0.6, 0.8, 0.9, 1.0]:
    edges = int((D[np.triu_indices(30, 1)] <= 2 * r).sum())
    cech, rips = int((meb <= r).sum()), int((diam <= 2 * r).sum())
    ok = np.all(meb[diam <= 2 * r] <= jung * r + 1e-12)   # Rips(r) inside Cech(1.155 r)
    print(f"r = {r:.1f}: edges {edges:3d}, Cech triangles {cech:4d}, "
          f"Rips triangles {rips:4d}, Rips(r) in Cech({jung:.3f} r): {ok}")

The 48 extra triangles at r = 0.9 are acute triangles spanning the circle. An inscribed equilateral triangle has side \sqrt 3: the Rips complex fills it from r = 0.866, the Čech complex only from its circumradius r = 1. The Rips complex closes the hole earlier, which ch. 4 measures as the death of a loop.

1.4 Alpha complexes and triangulated images

These sets are convex with union X^r, so A(X, r) \simeq X^r by the nerve theorem. In general position, the alpha complex is a subcomplex of the Delaunay triangulation: O(n) simplices in the plane, against up to 2^n for the Čech complex.

In 3D, each cube is cut into six tetrahedra; cubical complexes, which keep squares and cubes, are the other common choice.

1.5 Filtrations

The Čech, Rips and alpha complexes form filtrations as r grows; for the Rips complex, f(\sigma) is half the diameter of \sigma. For scalar fields, the central example is the following.

Superlevel sets use the upper-star filtration, with \min instead of \max and the order reversed.

2. Simplicial homology

Homology counts the holes of a complex in each dimension — components, loops, cavities — by linear algebra on its simplices. This chapter defines it over \mathbb{F}_2 and computes it by Gaussian elimination.

2.1 Chains, cycles and boundaries

For a triangle, \partial [a, b, c] = ab + bc + ca, and the boundary of this loop is (a + b) + (b + c) + (c + a) = 0. This is general.

A k-cycle that is not a boundary encloses a k-dimensional hole; two cycles are homologous when their difference bounds.

2.2 First computations

The torus of Example 1.3 has \chi = 0, and Code 2.1 finds \beta = (1, 2, 1): two independent loops and one cavity.

2.3 Maps and invariance

A simplex is sent to zero in C_k(L) when its image has fewer vertices. Commutation with \partial makes \varphi_* well defined: cycles go to cycles and boundaries to boundaries.

With the nerve theorem, the homology of the Čech complex is the homology of the union of balls. Contiguity is the tool of the stability proof for Rips filtrations in ch. 7.

2.4 Other coefficients

Over \mathbb{Q} or \mathbb{Z}, the vertices of each simplex are ordered, and faces receive signs: \partial [v_0, \dots, v_k] = \sum_{i=0}^{k} (-1)^i\, [v_0, \dots, \widehat{v_i}, \dots, v_k]. Lemma 2.2 still holds, because the two copies of each (k-2)-face now carry opposite signs (Exercise 9).

The Betti numbers can depend on the field. On the projective plane, a loop that goes around twice bounds, but the loop itself does not: over \mathbb{Z}, H_1 is \mathbb{Z}/2. Over \mathbb{F}_2 this torsion is visible as \beta_1 = 1, and over \mathbb{Q} it disappears.

Over \mathbb{F}_2 no orientation is needed and each column of a boundary matrix is a set of indices. TDA software therefore works over \mathbb{F}_2 or a small prime field.

import numpy as np
from itertools import combinations

def closure(top):
    """All faces of a list of top simplices, sorted by dimension."""
    faces = {tuple(sorted(f)) for s in top for k in range(1, len(s) + 1)
             for f in combinations(s, k)}
    return [sorted(f for f in faces if len(f) == k + 1) for k in range(len(top[0]))]

def boundary(S, k, signed):
    """Matrix of the boundary map C_k -> C_{k-1} (rows: (k-1)-faces)."""
    index = {f: i for i, f in enumerate(S[k - 1])}
    B = np.zeros((len(S[k - 1]), len(S[k])), dtype=np.int64)
    for j, s in enumerate(S[k]):
        for i in range(len(s)):
            B[index[s[:i] + s[i + 1:]], j] = (-1) ** i if signed else 1
    return B

def rank_f2(B):
    """Rank over Z/2 by Gaussian elimination."""
    A, r = B.copy() % 2, 0
    for c in range(A.shape[1]):
        piv = np.nonzero(A[r:, c])[0]
        if len(piv) == 0:
            continue
        A[[r, r + piv[0]]] = A[[r + piv[0], r]]
        rows = np.nonzero(A[:, c])[0]
        A[rows[rows != r]] ^= A[r]
        r += 1
        if r == A.shape[0]:
            break
    return r

def betti(top, field):
    S = closure(top)
    rk = [0] + [rank_f2(boundary(S, k, False)) if field == "F2"
                else np.linalg.matrix_rank(boundary(S, k, True).astype(float))
                for k in range(1, len(S))] + [0]
    return [int(len(S[k]) - rk[k] - rk[k + 1]) for k in range(len(S))]

n = 3                                    # 3 x 3 periodic grid: a 9-vertex torus
v = lambda i, j: (i % n) * n + (j % n)
torus = [t for i in range(n) for j in range(n) for t in
         [(v(i, j), v(i + 1, j), v(i + 1, j + 1)), (v(i, j), v(i, j + 1), v(i + 1, j + 1))]]
rp2 = [(0, 1, 3), (0, 1, 5), (0, 2, 4), (0, 2, 5), (0, 3, 4),
       (1, 2, 3), (1, 2, 4), (1, 4, 5), (2, 3, 5), (3, 4, 5)]   # 6-vertex projective plane
for name, top in [("torus", torus), ("projective plane", rp2)]:
    S = closure(top)
    chi = sum((-1) ** k * len(S[k]) for k in range(3))
    print(f"{name:16s} f-vector {[len(s) for s in S]}, chi = {chi}, "
          f"Betti over F2 {betti(top, 'F2')}, over Q {betti(top, 'Q')}")

Gaussian elimination on n simplices costs O(n^3) operations in the worst case. Ch. 4 reorganises this elimination so that one pass over a filtration gives the homology of every subcomplex at once.

3. Smooth and discrete Morse theory

A scalar field changes the topology of its sublevel sets only at critical points. Morse theory makes this precise; its piecewise-linear and discrete versions turn it into algorithms.

3.1 Smooth Morse theory

Minima have index 0 and maxima index d; on a surface, saddles have index 1. Generic Morse functions are dense among smooth functions.

For the torus, 1 - 2 + 1 = 0 = \chi, and m = (1, 2, 1) equals the Betti numbers: the height function is perfect.

3.2 Piecewise-linear Morse theory

Data come as values on the vertices of a triangulation, extended linearly. Ties are broken by the vertex index (simulation of simplicity, Edelsbrunner and Mücke, 1990), so all values are distinct.

This classification follows Banchoff (1967): entering v adds a cone over \mathrm{Lk}^-(v), which merges k fronts when the lower link has k arcs.

On a surface, the term 1 - \chi(\mathrm{Lk}^-(v)) of a vertex is:

• 1 for a minimum (empty lower link) and for a maximum (a cycle, of Euler characteristic 0);
• 0 for a regular vertex (an arc, of Euler characteristic 1);
• 1 - k for a saddle with k lower arcs.

Hence:

import numpy as np

def critical_points(f):
    """PL critical vertices of f on the Freudenthal triangulation of a periodic grid.
    Ties are broken by the vertex index (simulation of simplicity)."""
    n, m = f.shape
    key = f.ravel().argsort(kind="stable").argsort().reshape(n, m)  # rank of each vertex
    ring = [(1, 0), (1, 1), (0, 1), (-1, 0), (-1, -1), (0, -1)]  # cyclic link of a vertex
    low = np.stack([np.roll(key, (-di, -dj), axis=(0, 1)) < key for di, dj in ring])
    arcs = (low & ~np.roll(low, 1, axis=0)).sum(axis=0)     # number of lower arcs in the hexagon
    n_low = low.sum(axis=0)
    minima, maxima = n_low == 0, n_low == 6
    saddle_mult = np.where(~minima & ~maxima, arcs - 1, 0)  # a k-fold saddle has k lower arcs
    return minima.sum(), saddle_mult.sum(), maxima.sum()

n = 64
x, y = np.meshgrid(np.arange(n) / n, np.arange(n) / n, indexing="ij")
rng = np.random.default_rng(0)
for noise in [0.0, 0.05, 0.2]:
    f = np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y) + noise * rng.normal(size=(n, n))
    m0, m1, m2 = critical_points(f)
    print(f"noise {noise:.2f}: minima {m0:4d}, saddles {m1:4d}, maxima {m2:4d}, "
          f"m0 - m1 + m2 = {m0 - m1 + m2}")

Noise multiplies the critical points by hundreds, while their alternating sum cannot move. Deciding which of them matter is the job of persistence (ch. 4).

3.3 Discrete Morse theory

Forman’s theory replaces the gradient by a matching between simplices, which suits exact combinatorial algorithms.

Robins, Wood and Sheppard (2011) extend this pairing to all dimensions, with critical simplices matching the topology changes of the lower-star filtration; their V-paths replace integral lines in ch. 6.

4. Persistent homology and persistence diagrams

Data come with a scale or a threshold that nobody knows how to choose. Persistence follows every hole across all scales at once, records when it appears and disappears, and summarizes the result stably.

4.1 Persistence modules and barcodes

Let \emptyset = K_0 \subset K_1 \subset \dots \subset K_m = K be a filtration in which K_i adds one simplex \sigma_i to K_{i-1} (ties broken arbitrarily, faces before cofaces).

\beta_k^{i,j} counts the k-dimensional holes of K_i that are still holes in K_j.

An interval [b, d) of indices becomes the interval [f(\sigma_b), f(\sigma_d)) of scales; intervals of length zero are discarded.

4.2 The elder rule

For the Rips filtration, the finite H_0 bars die at half the edge lengths of a minimum spanning tree: the sweep is Kruskal’s algorithm (Exercise 14).

4.3 The reduction algorithm

Order the simplices as in the filtration, and let \partial be the m \times m boundary matrix over \mathbb{F}_2: column j holds the faces of \sigma_j. For a non-zero column, \mathrm{low}(j) is the index of its lowest 1.

A simplex with R_i = 0 is positive: it creates a class. A simplex with R_j \neq 0 is negative: it kills the class created by \sigma_{\mathrm{low}(j)}, the youngest simplex of the cycle R_j.

import numpy as np
from itertools import combinations

def rips_filtration(P, rmax):
    """Simplices of dimension <= 2 of the Rips filtration, value r = diameter / 2."""
    D = np.linalg.norm(P[:, None] - P[None], axis=2) / 2
    n = len(P)
    S = [((i,), 0.0) for i in range(n)]
    S += [((i, j), D[i, j]) for i, j in combinations(range(n), 2) if D[i, j] <= rmax]
    S += [((i, j, k), max(D[i, j], D[i, k], D[j, k])) for i, j, k in combinations(range(n), 3)
          if max(D[i, j], D[i, k], D[j, k]) <= rmax]
    return sorted(S, key=lambda s: (s[1], len(s[0])))       # faces before cofaces

def persistence(S):
    """Standard column reduction over Z/2; returns (dim, birth, death) triples."""
    index = {s: i for i, (s, _) in enumerate(S)}
    low_of, cols, pairs, paired = {}, [], [], set()
    for j, (s, val) in enumerate(S):
        col = {index[s[:i] + s[i + 1:]] for i in range(len(s))} if len(s) > 1 else set()
        while col and max(col) in low_of:                  # add the column with the same pivot
            col ^= cols[low_of[max(col)]]
        cols.append(col)
        if col:
            i = max(col)
            low_of[i] = j
            paired |= {i, j}
            if val > S[i][1]:                              # skip zero-length intervals
                pairs.append((len(S[i][0]) - 1, S[i][1], val))
    pairs += [(len(s) - 1, val, np.inf) for j, (s, val) in enumerate(S)
              if j not in paired and len(s) < 3]           # essential classes (dim < top)
    return pairs

rng = np.random.default_rng(0)
theta = rng.uniform(0, 2 * np.pi, 30)
P = np.c_[np.cos(theta), np.sin(theta)] + 0.05 * rng.normal(size=(30, 2))
S = rips_filtration(P, rmax=1.2)
dgm = persistence(S)
h1 = sorted([(b, d) for k, b, d in dgm if k == 1], key=lambda p: p[0] - p[1])
print(f"{len(S)} simplices; H0: {sum(k == 0 for k, _, _ in dgm)} intervals, "
      f"H1: {len(h1)} intervals")
for b, d in h1[:3]:
    print(f"H1 interval [{b:.3f}, {d:.3f}), persistence {d - b:.3f}")

The loop is born when the largest gap between consecutive points closes, and dies near \sqrt 3 / 2 = 0.866, when the Rips complex fills an inscribed triangle (Code 1.1).

4.4 Persistence diagrams and their distances

Points far from the diagonal are long-lived features; the diagonal lets diagrams of different sizes be matched.

4.5 Stability

Cohen-Steiner, Edelsbrunner and Harer (2007) proved it for tame functions on triangulable spaces; Section 10.2 gives the algebraic version.

4.6 Images, noise and simplification

For an image, bright blobs are maxima, tracked by superlevel sets. In degree 0, each local maximum dies at the saddle level where it joins a higher peak: its persistence is its prominence.

import numpy as np

def superlevel_h0(f):
    """0-dim persistence of the superlevel sets of an image (4-neighbour grid),
    by union-find and the elder rule. Returns (birth, death) pairs, birth >= death."""
    n, m = f.shape
    flat = f.ravel()
    parent = -np.ones(n * m, dtype=int)                # -1 = pixel not yet in the filtration
    def find(a):
        while parent[a] != a:
            parent[a] = parent[parent[a]]
            a = parent[a]
        return a
    pairs = []
    for v in np.argsort(-flat, kind="stable"):         # decreasing values
        parent[v] = v
        i, j = divmod(v, m)
        for a, b in [(i - 1, j), (i + 1, j), (i, j - 1), (i, j + 1)]:
            if 0 <= a < n and 0 <= b < m and parent[a * m + b] >= 0:
                ru, rv = find(a * m + b), find(v)
                if ru != rv:                           # elder rule: the lower peak dies
                    young, old = (ru, rv) if flat[ru] < flat[rv] else (rv, ru)
                    if flat[young] > flat[v]:          # skip zero-length pairs
                        pairs.append((flat[young], flat[v]))
                    parent[young] = old
    pairs.append((flat.max(), -np.inf))                # the essential class
    return np.array(pairs)

n = 64
x, y = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")
bump = lambda cx, cy, h: h * np.exp(-((x - cx) ** 2 + (y - cy) ** 2) / (2 * 6.0 ** 2))
clean = bump(16, 16, 1.0) + bump(44, 20, 0.8) + bump(30, 48, 0.6)
img = clean + 0.05 * np.random.default_rng(0).normal(size=(n, n))
dgm = superlevel_h0(img)
pers = np.sort(dgm[:, 0] - dgm[:, 1])[::-1]
print(f"{len(dgm)} local maxima; five largest persistences: {np.round(pers[:5], 3)}")
print(f"points with persistence > 0.25: {(pers > 0.25).sum()}")

The gap between 0.64 and 0.172 separates objects from noise. Theorem 4.10 explains it: noise of sup-norm \delta only creates points of persistence at most 2\delta.

This bound is a worst case. Here \delta = 0.195, so it allows noise points up to 0.39, while the highest one reaches 0.172 (Exercise 6).

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import maximum_bipartite_matching

def sublevel_h0(f):
    """Finite 0-dim intervals (birth, death) of the sublevel sets of a 1D signal."""
    parent, pairs = -np.ones(len(f), dtype=int), []
    def find(a):
        while parent[a] != a:
            a = parent[a]
        return a
    for v in np.argsort(f, kind="stable"):
        parent[v] = v
        for u in (v - 1, v + 1):
            if 0 <= u < len(f) and parent[u] >= 0:
                ru, rv = find(u), find(v)
                if ru != rv:
                    young, old = (ru, rv) if f[ru] > f[rv] else (rv, ru)
                    if f[young] < f[v]:
                        pairs.append((f[young], f[v]))
                    parent[young] = old
    return np.array(pairs).reshape(-1, 2)

def bottleneck(A, B):
    """Bottleneck distance between finite diagrams, by thresholded bipartite matching.
    Rows: points of A, then diagonal copies of B; columns: points of B, then copies of A."""
    a, b = len(A), len(B)
    C = np.full((a + b, b + a), np.inf)
    C[:a, :b] = np.abs(A[:, None, :] - B[None, :, :]).max(axis=2)
    C[np.arange(a), b + np.arange(a)] = (A[:, 1] - A[:, 0]) / 2   # A_i to the diagonal
    C[a + np.arange(b), np.arange(b)] = (B[:, 1] - B[:, 0]) / 2   # diagonal to B_j
    C[a:, b:] = 0.0                                                # diagonal to diagonal
    values = np.unique(C[np.isfinite(C)])
    lo, hi = 0, len(values) - 1
    while lo < hi:                                     # smallest threshold with a perfect matching
        mid = (lo + hi) // 2
        match = maximum_bipartite_matching(csr_matrix(C <= values[mid]), perm_type="column")
        lo, hi = (mid + 1, hi) if (match < 0).any() else (lo, mid)
    return values[lo]

rng = np.random.default_rng(1)
t = np.linspace(0, 1, 400)
f = np.sin(6 * np.pi * t) + 0.1 * rng.normal(size=t.size)
for delta in [0.01, 0.05, 0.2]:
    g = f + rng.uniform(-delta, delta, size=t.size)
    Df, Dg = sublevel_h0(f), sublevel_h0(g)
    big = [(D[:, 1] - D[:, 0] > 1.5).sum() for D in (Df, Dg)]
    print(f"|f - g|_inf = {np.abs(f - g).max():.4f}: bottleneck = {bottleneck(Df, Dg):.4f}; "
          f"{len(Df)} and {len(Dg)} points, persistence > 1.5: {big[0]} and {big[1]}")

Removing a point of persistence p needs a perturbation of sup-norm at least p/2, its distance to the diagonal. On surfaces this bound is reached.

5. Reeb graphs, contour trees and Mapper

The level sets of a scalar field split and merge as the level varies. The Reeb graph records this history; Mapper computes a discrete version of it for point clouds.

5.1 Reeb graphs

For the height function of the torus (Example 3.3), the Reeb graph is a minimum, a saddle where the circle splits into two, two parallel arcs, a saddle where they merge, and a maximum: one loop, as g = 1.

5.2 Merge trees and contour trees

On a simply connected domain, such as a rectangular image, \beta_1 = 0, so the Reeb graph is a tree by Theorem 5.2: the contour tree.

Both are computed by the union-find sweep of Proposition 4.3, recording a node at every merge. A persistence pair is then a leaf together with the node where its branch meets an older branch: the diagram in degree 0 is the branch decomposition of the merge tree (Pascucci, Cole-McLaughlin and Scorzelli, 2004).

In scientific visualization, an arc of the contour tree is a family of contours without topological change: selecting it segments a feature. Cutting the branch decomposition at a persistence threshold simplifies the tree.

5.3 Mapper

When the clusters are connected components and the intervals shrink, Mapper converges to the Reeb graph (Munch and Wang, 2016). Carrière and Oudot (2018) describe which features a given cover keeps, and Carrière, Michel and Oudot (2018) choose the cover from data.

import numpy as np
from scipy.cluster.hierarchy import fcluster, linkage

def mapper(X, filt, n_int=8, overlap=0.3, eps=0.25):
    """Mapper graph: cover the range of the filter by overlapping intervals,
    cluster each preimage by single linkage at scale eps, link clusters that share points."""
    lo, hi = filt.min(), filt.max()
    length = (hi - lo) / (n_int - (n_int - 1) * overlap)
    nodes = []
    for k in range(n_int):
        a = lo + k * length * (1 - overlap)
        idx = np.nonzero((filt >= a) & (filt <= a + length))[0]
        if len(idx) == 1:
            nodes.append(set(idx))
        elif len(idx) > 1:
            labels = fcluster(linkage(X[idx], "single"), eps, criterion="distance")
            nodes += [set(idx[labels == c]) for c in np.unique(labels)]
    edges = {(i, j) for i in range(len(nodes)) for j in range(i + 1, len(nodes))
             if nodes[i] & nodes[j]}
    return nodes, edges

def components(n, edges):
    parent = list(range(n))
    def find(a):
        while parent[a] != a:
            a = parent[a]
        return a
    for i, j in edges:
        parent[find(i)] = find(j)
    return len({find(i) for i in range(n)})

rng = np.random.default_rng(0)
t = rng.uniform(0, 2 * np.pi, 600)
data = {"circle": np.c_[np.cos(t), np.sin(t)],
        "figure eight": np.c_[np.sin(t), np.sin(t) * np.cos(t)]}
for name, X in data.items():
    X = X + 0.03 * rng.normal(size=X.shape)
    nodes, edges = mapper(X, filt=X[:, 0])
    c = components(len(nodes), edges)
    print(f"{name:12s}: {len(nodes)} nodes, {len(edges)} edges, {c} component(s), "
          f"independent loops = {len(edges) - len(nodes) + c}")

Too few intervals merge branches and too little overlap disconnects the graph (Exercise 13). Nicolau, Levine and Carlsson (2011) used Mapper on breast-cancer expression profiles to find a subgroup with excellent survival.

6. Morse–Smale complexes

Following the gradient from every point partitions the domain into cells, each flowing from one minimum to one maximum. These cells segment images along ridges and valleys, and persistence decides which to merge.

6.1 Integral lines and cells

By the stable manifold theorem, if p has index \lambda, then D(p) is an open cell of dimension \lambda and A(p) an open cell of dimension d - \lambda.

The descending manifolds partition M, and so do the ascending ones.

For \cos 2\pi x + \cos 2\pi y on the torus, m = (1, 2, 1): the complex has 4 nodes, 8 arcs and 4 quadrangles.

6.2 Watersheds

The ascending manifolds A(m) of the minima are the catchment basins of f: the points whose steepest descent ends at m. Their boundaries form the watershed used in image segmentation (Vincent and Soille, 1991). Symmetrically, the descending manifolds D(M) of the maxima of an image are the regions that climb to each bright peak.

On a grid, steepest ascent replaces the gradient: each pixel points to its highest neighbour, and the regions are the trees of these pointers. Each local maximum then has a region, and noise makes hundreds of them.

6.3 Discrete Morse–Smale complexes

Integrating integral lines numerically is fragile at saddles. The discrete gradients of Section 3.3 replace them: V-paths from critical edges to critical vertices are the descending separatrices, and dual paths to critical triangles the ascending ones. Gyulassy, Bremer, Hamann and Pascucci (2008) and Shivashankar and Natarajan (2012) compute such complexes in 3D and in parallel; TTK implements them.

6.4 Persistence-driven simplification

A saddle and an adjacent extremum connected by a single separatrix can be cancelled: the two critical points disappear, and the cells on both sides merge. Cancelling pairs in increasing order of persistence gives a hierarchy of Morse–Smale complexes (Edelsbrunner, Harer and Zomorodian, 2003).

In terms of regions, cancelling a maximum with the saddle where it meets an older peak merges its region into the older one. This is exactly the elder rule of Proposition 4.3, applied with a threshold \tau.

import numpy as np
from sklearn.metrics import adjusted_rand_score

def segment(f, tau):
    """Descending manifolds of the maxima of f on the 4-neighbour grid, by steepest
    ascent, with the regions of maxima of persistence < tau merged (elder rule)."""
    n, m = f.shape
    flat = f.ravel()
    parent = -np.ones(n * m, dtype=int)
    def find(a):
        while parent[a] != a:
            parent[a] = parent[parent[a]]
            a = parent[a]
        return a
    for v in np.argsort(-flat, kind="stable"):          # decreasing values
        i, j = divmod(v, m)
        nb = [a * m + b for a, b in [(i - 1, j), (i + 1, j), (i, j - 1), (i, j + 1)]
              if 0 <= a < n and 0 <= b < m and parent[a * m + b] >= 0]
        if not nb:
            parent[v] = v                               # a new maximum
            continue
        parent[v] = find(max(nb, key=lambda u: flat[u]))  # steepest ascent
        for u in nb:
            ru, rv = find(u), find(v)
            if ru != rv and flat[v] > min(flat[ru], flat[rv]) - tau:   # low persistence
                young, old = (ru, rv) if flat[ru] < flat[rv] else (rv, ru)
                parent[young] = old
    return np.array([find(v) for v in range(n * m)]).reshape(n, m)

n = 64
x, y = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")
bump = lambda cx, cy, h: h * np.exp(-((x - cx) ** 2 + (y - cy) ** 2) / (2 * 6.0 ** 2))
clean = bump(16, 16, 1.0) + bump(44, 20, 0.8) + bump(30, 48, 0.6)
img = clean + 0.05 * np.random.default_rng(0).normal(size=(n, n))
truth = segment(clean, 0.0)
for tau in [0.0, 0.1, 0.25]:
    lab = segment(img, tau)
    ari = adjusted_rand_score(truth.ravel(), lab.ravel())
    print(f"tau = {tau:.2f}: {len(np.unique(lab)):3d} regions, "
          f"adjusted Rand index with the clean segmentation {ari:.3f}")
print("clean image:", len(np.unique(truth)), "regions")

The remaining disagreement lies in the flat background, where noise moves the boundaries.

7. Topological inference in high dimensions

A point cloud in \mathbb{R}^{100} cannot be looked at, but its homology can be estimated. This chapter gives conditions under which the homology of a sampled shape is recovered, bounds the effect of noise and outliers, and turns persistence into a clustering algorithm.

7.1 Offsets, reach and sampling

A circle of radius R has reach R: its medial axis is the centre. The boundary of a square has reach 0, since its medial axis, the diagonals, reaches the corners. Closed smooth submanifolds have positive reach.

The retraction is the straight segment from each point to its projection onto M; the nerve theorem transfers the result to the Čech complex. Chazal, Cohen-Steiner and Lieutier (2009) extend it to compact sets of positive \mu-reach, which include sets of zero reach such as polygons.

7.2 Inference with Rips complexes

The Rips complex is not homotopy equivalent to the union of balls, but two Rips complexes at different scales are enough.

The hypothesis holds for a submanifold sampled as in Theorem 7.2 when r \geq \varepsilon and 2r < \sqrt{3/5}\, \tau: the three offsets then retract onto M, compatibly with the inclusions. For samples of more general compact sets, Chazal and Oudot (2008) prove a variant with a longer chain of scales.

7.3 Stability of Rips persistence

This is due to Chazal, Cohen-Steiner, Guibas, Mémoli and Oudot (2009) and Chazal, de Silva and Oudot (2014); Section 10.2 gives the proof. Two \varepsilon-dense samples of the same shape thus have diagrams within \varepsilon of each other. But a single outlier can change d_H, hence the diagram, by an arbitrary amount.

7.4 Distance to a measure

Moving one of n sample points by a distance L changes W_2 by at most L / \sqrt n, hence the empirical DTM by at most L / \sqrt k, everywhere (Exercise 7). The sampling theorems of Section 7.1 extend to its sublevel sets.

import numpy as np
from scipy.spatial import cKDTree

def lower_star_h1(F):
    """Finite H1 intervals of the sublevel sets of F on the Freudenthal triangulation
    of a grid (lower-star filtration: a simplex enters with its highest vertex)."""
    n, m = F.shape
    v = lambda i, j: i * m + j
    S = [((v(i, j),), F[i, j]) for i in range(n) for j in range(m)]
    for i in range(n):
        for j in range(m):
            for a, b in [(i + 1, j), (i, j + 1), (i + 1, j + 1)]:
                if a < n and b < m:
                    S.append(((v(i, j), v(a, b)), max(F[i, j], F[a, b])))
            if i + 1 < n and j + 1 < m:
                for c in [(i + 1, j), (i, j + 1)]:
                    tri = (v(i, j), v(*c), v(i + 1, j + 1))
                    S.append((tuple(sorted(tri)), max(F[i, j], F[c], F[i + 1, j + 1])))
    S.sort(key=lambda s: (s[1], len(s[0])))
    index = {s: k for k, (s, _) in enumerate(S)}
    low_of, cols, out = {}, [], []
    for k, (s, val) in enumerate(S):
        col = {index[s[:i] + s[i + 1:]] for i in range(len(s))} if len(s) > 1 else set()
        while col and max(col) in low_of:
            col ^= cols[low_of[max(col)]]
        cols.append(col)
        if col:
            low_of[max(col)] = k
            if len(s) == 3 and val > S[max(col)][1]:
                out.append((S[max(col)][1], val))
    return sorted(out, key=lambda p: p[0] - p[1])

rng = np.random.default_rng(0)
theta = rng.uniform(0, 2 * np.pi, 500)
X = np.c_[np.cos(theta), np.sin(theta)] + 0.05 * rng.normal(size=(500, 2))
Y = np.r_[X, rng.uniform(-1.5, 1.5, size=(25, 2))]          # add 5 % of uniform outliers
g = np.linspace(-1.5, 1.5, 61)
Z = np.stack(np.meshgrid(g, g, indexing="ij"), axis=-1).reshape(-1, 2)
dist, _ = cKDTree(Y).query(Z, k=50)                        # 50 nearest points: mass m = 50 / 525
functions = {"distance, clean": cKDTree(X).query(Z)[0],
             "distance, outliers": dist[:, 0],
             "DTM, outliers": np.sqrt((dist ** 2).mean(axis=1))}
for name, f in functions.items():
    bars = lower_star_h1(f.reshape(61, 61))
    pers = [d - b for b, d in bars]
    print(f"{name:18s}: longest H1 interval [{bars[0][0]:.3f}, {bars[0][1]:.3f}), "
          f"persistence {pers[0]:.3f}; intervals longer than 0.05: {sum(p > 0.05 for p in pers)}")

Outliers inside the circle fill its hole early in the offsets of the distance function. For the DTM they carry too little mass: the hole dies at 0.885, close to the clean 0.857.

7.5 Persistence-based clustering (ToMATo)

An estimated density has many spurious local maxima, like the image of Code 4.2. ToMATo (Chazal, Guibas, Oudot and Skraba, 2013) keeps the prominent ones.

Under regularity assumptions, the diagram of \hat f is close to that of the true density, and the prominent peaks and their clusters are recovered (Chazal, Guibas, Oudot and Skraba, 2013).

import numpy as np
from scipy.spatial import cKDTree
from sklearn.cluster import KMeans
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score

def tomato(X, k):
    """ToMATo with tau = infinity: hill climbing on the k-NN graph of a k-NN density estimate.
    Returns the peak of each point, the merge history and the prominence of each merged peak."""
    dist, nbrs = cKDTree(X).query(X, k=k + 1)
    dens = 1.0 / dist[:, -1] ** X.shape[1]                     # k-NN density, up to a constant
    rank = np.empty(len(X), dtype=int)
    rank[np.argsort(-dens, kind="stable")] = np.arange(len(X))
    parent, peak, merged_into, prominence = -np.ones(len(X), dtype=int), {}, {}, {}
    def find(a):
        while parent[a] != a:
            parent[a] = parent[parent[a]]
            a = parent[a]
        return a
    for v in np.argsort(rank):                                 # decreasing density
        up = [u for u in nbrs[v, 1:] if rank[u] < rank[v]]
        parent[v] = v if not up else find(min(up, key=lambda u: rank[u]))
        peak[v] = parent[v]                                    # the hill v climbs to
        for u in up:
            ru, rv = find(u), find(v)
            if ru != rv:                                       # elder rule
                young, old = (ru, rv) if rank[ru] > rank[rv] else (rv, ru)
                merged_into[young], prominence[young] = old, dens[young] - dens[v]
                parent[young] = old
    return peak, merged_into, prominence

def labels_at(peak, merged_into, prominence, tau):
    """Merge every peak of prominence < tau into its elder; one label per point."""
    def resolve(p):
        while p in merged_into and prominence[p] < tau:
            p = merged_into[p]
        return p
    return np.array([resolve(peak[v]) for v in range(len(peak))])

X, y = make_moons(n_samples=600, noise=0.06, random_state=0)
km = KMeans(n_clusters=2, n_init=10, random_state=0).fit_predict(X)
print(f"k-means with 2 clusters: ARI = {adjusted_rand_score(y, km):.3f}")
for k in [15, 30]:
    peak, merged_into, prominence = tomato(X, k)
    top = np.sort(list(prominence.values()))[::-1]
    gap = np.argmax(top[:-1] - top[1:])                        # largest gap of the diagram
    labels = labels_at(peak, merged_into, prominence, tau=(top[gap] + top[gap + 1]) / 2)
    print(f"ToMATo, k = {k}: prominences / largest {np.round(top[:3] / top[0], 3)}, "
          f"{len(np.unique(labels))} clusters, ARI = {adjusted_rand_score(y, labels):.3f}")

The density estimate must be smooth enough. When no gap stands out, as at k = 15, the diagram warns that \tau is not determined by the data.

8. Topological signatures and features for learning

A persistence diagram is a multiset, not a vector: it cannot be averaged or fed to a linear model as it stands. This chapter maps diagrams to vectors and kernels without losing stability, and differentiates through them.

8.1 The space of diagrams

Diagrams with d_B or W_q form a metric space. Fréchet means exist (Mileyko, Mukherjee and Harer, 2011) but need not be unique, and the space does not embed bi-Lipschitz into a Hilbert space (Carrière and Bauer, 2019). Learning therefore goes through stable vectorizations or kernels.

8.2 Betti curves and landscapes

The Betti curve counts every short bar and is not stable. Landscapes are.

Landscapes can be averaged; mean landscapes satisfy a law of large numbers and a central limit theorem (Bubenik, 2015).

8.3 Persistence images

The weight makes the image continuous in the diagram, and persistence images are stable for the 1-Wasserstein distance (Adams and coauthors, 2017).

8.4 Kernels on diagrams

The first kernel is an inner product, hence positive definite. For the second, transport on a line is an L^1 distance between cumulative distribution functions, which is conditionally negative definite; by Schoenberg’s theorem k(D, E) = \exp\big( -SW(D, E) / 2\sigma^2 \big) is positive definite, ready for a support vector machine.

import numpy as np
from itertools import combinations
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

def h1_diagram(P, rmax=1.2):
    """Finite H1 intervals of the Rips filtration (r = diameter / 2), as in Code 4.1."""
    D = np.linalg.norm(P[:, None] - P[None], axis=2) / 2
    n = len(P)
    S = [((i,), 0.0) for i in range(n)]
    S += [((i, j), D[i, j]) for i, j in combinations(range(n), 2) if D[i, j] <= rmax]
    S += [((i, j, k), max(D[i, j], D[i, k], D[j, k])) for i, j, k in combinations(range(n), 3)
          if max(D[i, j], D[i, k], D[j, k]) <= rmax]
    S.sort(key=lambda s: (s[1], len(s[0])))
    index = {s: i for i, (s, _) in enumerate(S)}
    low_of, cols, out = {}, [], []
    for j, (s, val) in enumerate(S):
        col = {index[s[:i] + s[i + 1:]] for i in range(len(s))} if len(s) > 1 else set()
        while col and max(col) in low_of:
            col ^= cols[low_of[max(col)]]
        cols.append(col)
        if col:
            low_of[max(col)] = j
            if len(s) == 3 and val > S[max(col)][1]:
                out.append((S[max(col)][1], val))
    return np.array(out).reshape(-1, 2)

def landscape(dgm, grid, K=2):
    """First K persistence landscapes on a grid: k-th largest tent max(0, min(t - b, d - t))."""
    tents = np.maximum(0, np.minimum(grid[None] - dgm[:, :1], dgm[:, 1:] - grid[None]))
    tents = np.vstack([np.sort(tents, axis=0)[::-1], np.zeros((K, len(grid)))])
    return tents[:K].ravel()

rng = np.random.default_rng(0)
def sample(label):
    """30 noisy points on a circle (label 0) or on an arc missing a quarter (label 1)."""
    span = 2 * np.pi if label == 0 else 1.5 * np.pi
    t = rng.uniform(0, span, 30)
    R = rng.uniform(0.8, 1.2)
    return R * np.c_[np.cos(t), np.sin(t)] + 0.05 * rng.normal(size=(30, 2))

labels = np.repeat([0, 1], 40)
clouds = [sample(c) for c in labels]
grid = np.linspace(0, 1.2, 60)
topo = np.array([landscape(h1_diagram(P), grid) for P in clouds])
hist = np.array([np.histogram(np.linalg.norm(P[:, None] - P[None], axis=2)[np.triu_indices(30, 1)],
                              bins=30, range=(0, 2.6), density=True)[0] for P in clouds])
for name, F in [("H1 landscapes", topo), ("distance histogram", hist)]:
    acc = cross_val_score(LogisticRegression(C=10, max_iter=5000), F, labels, cv=5)
    print(f"{name:18s}: 5-fold accuracy {acc.mean():.3f} +/- {acc.std():.3f}")

The histogram sees the missing quarter only through a few long distances, blurred by the random radius; the landscape sees whether a long H_1 bar exists.

8.5 Learned layers and differentiable persistence

PersLay (Carrière, Chazal, Ike, Lacombe, Royer and Umeda, 2020) is a neural layer of the form D \mapsto \operatorname{op}_{p \in D}\; w(p)\, \varphi(p), where \operatorname{op} is a sum, a maximum or a k-th largest value, and w and \varphi are learned; landscapes and persistence images are special cases.

Persistence can also sit inside a model, when a filtration depends on parameters \theta such as pixel values or point positions.

Stochastic subgradient descent converges for such losses (Carrière, Chazal, Glisse, Ike, Kannan and Umeda, 2021). They enforce connectivity in segmentation (Hu, Li, Samaras and Chen, 2019) and regularize models (Brüel-Gabrielsson, Nelson, Dwaraknath and Skraba, 2020).

import numpy as np
import torch

torch.set_num_threads(1)

def h0_pairs(f):
    """Index pairs (birth vertex, death vertex) of the finite 0-dim sublevel intervals."""
    parent, pairs = -np.ones(len(f), dtype=int), []
    def find(a):
        while parent[a] != a:
            a = parent[a]
        return a
    for v in np.argsort(f, kind="stable"):
        parent[v] = v
        for u in (v - 1, v + 1):
            if 0 <= u < len(f) and parent[u] >= 0:
                ru, rv = find(u), find(v)
                if ru != rv:
                    young, old = (ru, rv) if f[ru] > f[rv] else (rv, ru)
                    if f[young] < f[v]:
                        pairs.append((young, v))
                    parent[young] = old
    return np.array(pairs, dtype=int).reshape(-1, 2)

rng = np.random.default_rng(0)
t = np.linspace(0, 1, 300)
f0 = torch.tensor(np.sin(6 * np.pi * t) + 0.15 * rng.normal(size=t.size))
x = f0.clone().requires_grad_(True)
opt = torch.optim.Adam([x], lr=0.01)
keep = 2                                          # keep the two most persistent finite pairs
print("local minima before:", len(h0_pairs(f0.numpy())) + 1)
for step in range(300):
    pairs = h0_pairs(x.detach().numpy())          # the pairing is piecewise constant in x
    pers = x[pairs[:, 1]] - x[pairs[:, 0]]        # death - birth, differentiable in x
    small = torch.sort(pers, descending=True).values[keep:]
    loss = (small ** 2).sum() + 0.1 * ((x - f0) ** 2).mean()
    opt.zero_grad()
    loss.backward()
    opt.step()
pairs = h0_pairs(x.detach().numpy())
pers = np.sort(x.detach().numpy()[pairs[:, 1]] - x.detach().numpy()[pairs[:, 0]])[::-1]
print(f"after 300 steps: largest persistences {', '.join(f'{p:.3f}' for p in pers[:4])}, "
      f"pairs above 0.05: {(pers > 0.05).sum()}, |x - f0|_inf = {(x - f0).abs().max():.3f}")

9. Implementations in Python and C++

The algorithms of this course are simple, but naive implementations do not scale: a Rips complex has a cubic number of triangles, and the reduction is cubic in the number of simplices. This chapter, for self-study alongside the labs, presents the data structures and the algebraic shortcuts used by the Gudhi and TTK libraries, and runs two of them.

9.1 Data structures

A simplex of dimension k is found in O(k) steps along its path, and its faces are reached by removing one vertex. A Rips complex is built by computing the neighbourhood graph and then expanding it: the cofaces of a node are the common neighbours of its vertices.

Three remedies keep the size manageable: stop the filtration at a radius r_{\max}; use the alpha complex, of size O(n) in the plane and O(n^2) in the worst case in 3D; or simplify the complex while preserving its persistence, for instance by edge collapses (Boissonnat and Pritam, 2020).

9.2 Faster reductions

Processing the dimensions from the top down lets every pivot found in dimension k + 1 clear a column of dimension k: this is the twist algorithm (Chen and Kerber, 2011).

For a Rips filtration, cohomology runs the clearing in the favourable direction: dimension 0 first, then dimension 1, so that the many triangles that create no class are never touched. Ripser (Bauer, 2021) combines cohomology, clearing, an implicit coboundary matrix computed on the fly, and apparent pairs detected without reduction.

import time
import numpy as np
from itertools import combinations

rng = np.random.default_rng(0)
theta = rng.uniform(0, 2 * np.pi, 30)
P = np.c_[np.cos(theta), np.sin(theta)] + 0.05 * rng.normal(size=(30, 2))
D = np.linalg.norm(P[:, None] - P[None], axis=2) / 2
S = [((i,), 0.0) for i in range(30)]
S += [((i, j), D[i, j]) for i, j in combinations(range(30), 2)]
S += [((i, j, k), max(D[i, j], D[i, k], D[j, k])) for i, j, k in combinations(range(30), 3)]
S.sort(key=lambda s: (s[1], len(s[0])))                 # the full 2-skeleton, 4525 simplices
index = {s: i for i, (s, _) in enumerate(S)}
bd = [{index[s[:i] + s[i + 1:]] for i in range(len(s))} if len(s) > 1 else set() for s, _ in S]
cobd = [set() for _ in S]
for j, faces in enumerate(bd):
    for i in faces:
        cobd[i].add(j)
by_dim = [[j for j, (s, _) in enumerate(S) if len(s) == d + 1] for d in range(3)]

def homology(order, clearing):
    """Reduce boundary columns (pivot = largest row) in the given order."""
    low_of, cols, adds, t0 = {}, {}, 0, time.perf_counter()
    for j in order:
        if clearing and j in low_of:
            continue                                     # positive simplex: column is zero
        col = set(bd[j])
        while col and max(col) in low_of:
            col ^= cols[low_of[max(col)]]
            adds += 1
        cols[j] = col
        if col:
            low_of[max(col)] = j
    return {(i, j) for i, j in low_of.items()}, adds, time.perf_counter() - t0

def cohomology():
    """Reduce coboundary columns (pivot = smallest row), dimension 0 then 1, with clearing."""
    piv_of, cols, adds, t0 = {}, {}, 0, time.perf_counter()
    for d in (0, 1):
        for j in reversed(by_dim[d]):
            if j in piv_of:
                continue                                 # cleared: j already died
            col = set(cobd[j])
            while col and min(col) in piv_of:
                col ^= cols[piv_of[min(col)]]
                adds += 1
            cols[j] = col
            if col:
                piv_of[min(col)] = j
    return {(j, k) for k, j in piv_of.items()}, adds, time.perf_counter() - t0

ref, _, _ = homology(range(len(S)), False)
runs = [("homology, standard", homology(range(len(S)), False)),
        ("homology, twist", homology(by_dim[2] + by_dim[1], True)),
        ("cohomology, clearing", cohomology())]
for name, (pairs, adds, sec) in runs:
    print(f"{name:21s}: same pairs {pairs == ref}, {adds:5d} column additions, {sec:.3f} s")

The twist saves little here because the expensive columns are the 4060 triangles, which nothing clears in homology. In cohomology, the triangles are never reduced at all: only the 406 edges that create loops are processed, and most of their columns are reduced already.

9.3 Libraries: Gudhi and TTK

Gudhi [7] is a C++ library with a Python interface, developed at Inria. It covers the high-dimensional side of the course: simplex trees, Rips, Čech, alpha, witness and cubical complexes, persistence computation, bottleneck and Wasserstein distances, the representations of ch. 8 in a scikit-learn style, ToMATo and cover complexes such as Mapper.

TTK [8] (Tierny, Favelier, Levine, Gueunet and Michaux, 2018) is a C++ library for scalar data on 2D and 3D meshes and grids, integrated with ParaView and VTK, with Python bindings. It covers the low-dimensional side: critical points, persistence diagrams and curves, topological simplification, merge and contour trees, Reeb graphs and Morse–Smale complexes, with parallel implementations.

The labs of the official course use these libraries; the code of this page rewrites the needed algorithms in a few dozen lines each, which makes every step visible and runs without installation.

9.4 A union-find in C++

The sweep of Proposition 4.3 is the inner loop of merge trees, watersheds and ToMATo. In C++ it runs on millions of vertices per second.

// 0-dimensional sublevel-set persistence of a sampled signal: union-find + elder rule.
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <numeric>
#include <vector>

struct UnionFind {
  std::vector<int> parent;
  explicit UnionFind(int n) : parent(n, -1) {}
  int find(int a) {                       // path halving
    while (parent[a] != a) a = parent[a] = parent[parent[a]];
    return a;
  }
};

int main() {
  const int n = 400;
  const double pi = std::acos(-1.0);
  std::vector<double> f(n);
  for (int i = 0; i < n; ++i) {           // a deterministic "noisy" signal
    double t = i / (n - 1.0);
    f[i] = std::sin(6 * pi * t) + 0.1 * std::sin(173 * t) + 0.05 * std::sin(517 * t);
  }
  std::vector<int> order(n);
  std::iota(order.begin(), order.end(), 0);
  std::stable_sort(order.begin(), order.end(), [&](int a, int b) { return f[a] < f[b]; });

  UnionFind uf(n);
  std::vector<std::pair<double, double>> pairs;   // (birth, death)
  for (int v : order) {
    uf.parent[v] = v;
    for (int u : {v - 1, v + 1}) {
      if (u < 0 || u >= n || uf.parent[u] < 0) continue;
      int ru = uf.find(u), rv = uf.find(v);
      if (ru == rv) continue;
      int young = f[ru] > f[rv] ? ru : rv, old = young == ru ? rv : ru;
      if (f[young] < f[v]) pairs.emplace_back(f[young], f[v]);
      uf.parent[young] = old;               // elder rule
    }
  }
  std::sort(pairs.begin(), pairs.end(), [](auto& a, auto& b) {
    return a.second - a.first > b.second - b.first;
  });
  std::printf("%zu finite pairs; the three most persistent:\n", pairs.size());
  for (int k = 0; k < 3; ++k)
    std::printf("  [%.3f, %.3f)  persistence %.3f\n", pairs[k].first, pairs[k].second,
                pairs[k].second - pairs[k].first);
}

10. Complements for self-study

This chapter collects the proofs that the sessions only sketch: the pairing lemma behind the reduction algorithm, the algebraic stability that covers both images and point clouds, collapses for discrete Morse theory, the nerve theorem, and the statistics of persistence diagrams. Read each section after the session that states the result.

10.1 The pairing lemma

Let \partial be the m \times m boundary matrix of a filtration (Section 4.3).

For indices i, j, let r_\partial(i, j) be the rank of its lower-left corner: the submatrix made of the rows i to m and the columns 1 to j.

To read the barcode from the pairs (Theorem 4.5), recall that V is upper triangular.

• If R_i = 0, the column V_i is a cycle whose last simplex is \sigma_i: a class is born at i.
• If moreover i = \mathrm{low}(j), the boundary R_j also ends with \sigma_i, so V_i - R_j is a cycle of K_{i-1}.
• In K_j, the new class therefore equals an older one: it dies at j.

Counting these classes gives the persistent Betti numbers of Theorem 4.2 (Edelsbrunner and Harer, 2010, Chapter VII).

10.2 Interleavings and the isometry theorem

10.3 Collapses and discrete Morse theory

To prove Theorem 3.10, order the simplices so that faces come before cofaces and each gradient pair (\sigma, \tau) is consecutive; such an order exists exactly when there is no closed V-path. Build K in this order.

• A pair (\sigma, \tau) enters while the other cofaces of \sigma are still absent, so \sigma is a free face of \tau: the step is an inverse elementary collapse and keeps the homotopy type.
• A critical simplex of dimension p enters along its whole boundary: the step attaches a p-cell.

Hence K is homotopy equivalent to a CW complex with one p-cell per critical p-simplex (Forman, 1998).

10.4 The nerve theorem for convex sets

Let \mathcal{U} = (U_i) be a finite family of open convex sets with union U, and N its nerve.

• From the nerve to the union. For each simplex J of N, choose a point p_J in the intersection of the sets U_j, j \in J.

  A simplex of the barycentric subdivision is a chain J_0 \subset \dots \subset J_k. Map its vertices to the points p_{J_i} and extend linearly. These points all lie in the convex set U_j for any j \in J_0, so the image lies in the union: this defines g.

• From the union to the nerve. Take a partition of unity (\phi_i) on U subordinate to the open cover (U_i), and set h(x) = \sum_i \phi_i(x)\, e_i, where e_i is the vertex i.

  The indices with \phi_i(x) > 0 label sets that contain x, so they form a simplex and h(x) \in |N|.

• Homotopies. The point h(x) lies in a simplex J_0 \subset \dots \subset J_k of the subdivision with \phi_j(x) > 0 for j \in J_0.

  So x and g(h(x)) both lie in the convex set U_j, and the straight segment between them stays in U: this gives g \circ h \simeq \mathrm{id}_U.

  A similar argument inside the simplices of N gives h \circ g \simeq \mathrm{id}.

Details are in Hatcher (2002, Section 4.G). Closed convex sets, such as the closed balls of the Čech complex, need one more step: Borsuk (1948) proves the nerve theorem for finite closed covers with contractible intersections, and Björner (1995) surveys both versions.

10.5 Statistics of persistence diagrams

Let X_n be n i.i.d. points from a measure \mu supported on a compact set M. Assume that, for some constants a, b > 0 and all x \in M, \mu(\bar B(x, r)) \geq \min(1, a r^b).

Fasy, Lecci, Rinaldo, Wasserman, Balakrishnan and Singh (2014) turn the same bound into confidence sets: a bootstrap estimate \hat c of the quantile of d_H(X_n, M) defines a band of width 2 \hat c along the diagonal, and the points outside the band are significant features. Chazal, Fasy, Lecci, Rinaldo and Wasserman (2015) build bootstrap confidence bands for mean landscapes (ch. 8).

11. Research directions and open problems

Topological data analysis is an active field, and the research seminar of the course is an entry point to it. This chapter, for self-study, points to directions that a project or an internship can start from, with references that the References list in full.

Multiparameter persistence

Data often come with two parameters at once — a scale and a density, or a scale and a time. Carlsson and Zomorodian (2009) showed that modules indexed by a grid of dimension two or more admit no complete discrete invariant like the barcode. Lesnick (2015) extended the interleaving distance to them, and Lesnick and Wright (2015) developed the RIVET software to explore two-parameter modules interactively. Carrière and Blumberg (2020) proposed multiparameter persistence images for learning.

Open problems: invariants that are stable, computable at scale and informative enough for statistics and learning; algorithms for minimal presentations in more than two parameters.

Statistics and learning with topology

The statistical theory of diagrams — convergence rates, bootstrap confidence sets, central limit theorems for landscapes (Section 10.5, ch. 8) — is surveyed by Chazal and Michel (2021). Learned filtrations go one step further: Hofer, Graf, Rieck, Niethammer and Kwitt (2020) learn the filtration of a graph end to end, and Horn, De Brouwer, Moor, Moreau, Rieck and Borgwardt (2022) add persistence to graph neural networks. Topological losses and regularizers (ch. 8) are now used in segmentation and representation learning (Moor, Horn, Rieck and Borgwardt, 2020).

Open problems: generalization guarantees for models with topological layers; the optimization landscape of topological losses, which are non-smooth and non-convex; filtrations that remain informative in very high dimension.

Topology for scientific visualization at scale

Simulations produce ensembles of large fields, and comparing their topological summaries requires geometry on the space of summaries. Vidal, Budin and Tierny (2020) compute Wasserstein barycenters of persistence diagrams progressively; Pont, Vidal, Delon and Tierny (2022) define Wasserstein distances, geodesics and barycenters of merge trees. Soler, Plainchault, Conche and Tierny (2018) compress fields lossily while controlling the change of their diagram.

Open problems: in-situ topological analysis during simulations, distributed computation of merge trees and Morse–Smale complexes on supercomputers, and summaries of uncertain or time-varying fields.

Shapes and transforms

The persistent homology transform of Turner, Mukherjee and Boyer (2014) records the diagrams of the height functions of a shape in all directions; for triangulated surfaces in \mathbb{R}^3, this collection of diagrams determines the shape. Related Euler characteristic transforms are used in biology and medical imaging.

Open problems: the number of directions needed in practice, and statistical models on these transforms.

Applications in the sciences

Persistent homology describes the medium-range order of amorphous solids (Hiraoka and coauthors, 2016), the filaments of the cosmic web (Sousbie, 2011), and the topology of anatomical structures in segmentation (Hu, Li, Samaras and Chen, 2019). In each case, a stable summary of multi-scale structure answers a question that pointwise statistics miss.

Open problems: the interpretation of diagram points in terms of the underlying physics or biology, and the design of filtrations adapted to each domain.

Exercises

The exercises are sorted by level: ★ applies the course, ★★ combines results or proves a short statement, ★★★ is open-ended. The session plan says which ones to do after each session; exercises 6, 13, 15 and 16 are computational.

1. ★ Counting simplices. (a) Compute the f-vector, the Euler characteristic and the Betti numbers of the boundary of an octahedron. (b) The five triangles \{i, i+1, i+2\}, indices modulo 5, form a complex on 5 vertices. Compute its f-vector and its Euler characteristic, find its boundary edges, and compute its Betti numbers over \mathbb{F}_2.

   Solution: (a) (6, 12, 8), \chi = 2, Betti numbers (1, 0, 1): it is a sphere.

   2. The f-vector is (5, 10, 5) and \chi = 0. Each edge \{i, i+1\} lies in two triangles and each edge \{i, i+2\} in one, so the boundary is the single 5-cycle 0, 2, 4, 1, 3: this is a Möbius band, since an annulus has two boundary circles.

   No non-empty set of triangles has zero boundary, so \beta_2 = 0; the complex is connected, so \beta_0 = 1. Theorem 2.7 then gives \beta_1 = 1.

2. ★ Čech or Rips? (a) For the vertices of an equilateral triangle of side 1, find the radii at which the edges and the triangle enter the Čech and the Rips filtrations, and the range of radii over which the edges form a loop that no triangle fills. (b) Same question for the vertices of the unit square.

   Solution: (a) All edges enter at r = 1/2 in both filtrations. The Rips triangle enters at once: the loop is never unfilled. The Čech triangle waits for the circumradius 1/\sqrt 3 = 0.577, so the loop is unfilled on [0.5, 0.577): this is the H_1 interval of ch. 4. The ratio 1.155 = \sqrt{4/3} is the constant of Theorem 1.7 in the plane, which is therefore optimal.

   2. The sides enter at 0.5 and the diagonals at \sqrt 2 / 2 = 0.707. In the Rips filtration, the diagonals fill everything at once; in the Čech filtration, the four disks meet at the centre at 0.707, and the right triangles enter at half their hypotenuse, also 0.707. Both leave the loop of the four sides unfilled on [0.5, 0.707).

3. ★ Classifying a vertex. A vertex of a triangulated surface has a hexagonal link whose values, in cyclic order, are 3, 7, 4, 8, 2, 9. Classify the vertex and compute its term 1 - \chi(\mathrm{Lk}^-) in Theorem 3.6 when its own value is 1, 2.5, 5, 7.5 and 10.

   Solution: at 1, the lower link is empty: a minimum, term +1. At 2.5, the lower link is the single vertex of value 2: regular, term 0. At 5, the lower vertices 3, 4, 2 alternate with 7, 8, 9: three arcs, a monkey saddle of multiplicity 2, term -2.

   At 7.5, the lower vertices 3, 7, 4 form one arc and 2 another: a simple saddle, term -1. At 10, the lower link is the whole hexagon: a maximum, term +1.

4. ★ A sublevel-set diagram by hand. The path v_1, \dots, v_7 carries the values f = (3, 1, 4, 2, 6, 0, 5). (a) Compute \mathrm{Dgm}_0(f) with the elder rule. (b) Do the same for g, equal to f except g(v_4) = 3.5, and check Theorem 4.10.

   Solution: (a) The minima are v_6 (0), v_2 (1) and v_4 (2). At v_3 (4), the component born at 2 dies; at v_5 (6), the component born at 1 dies. Hence \mathrm{Dgm}_0(f) = \{ (2, 4),\ (1, 6),\ (0, \infty) \}.

   2. \mathrm{Dgm}_0(g) = \{ (3.5, 4), (1, 6), (0, \infty) \}. Matching (2, 4) with (3.5, 4) costs 1.5, while sending both to the diagonal costs \max(1, 0.25) = 1. So d_B = 1 \leq \| f - g \|_\infty = 1.5.

5. ★ Bottleneck and Wasserstein by hand. Compute d_B and W_1 between D = \{(0, 2), (1, 5)\} and D' = \{(0.5, 2), (1, 4)\}.

   Solution: matching (0, 2) with (0.5, 2) costs 0.5, and (1, 5) with (1, 4) costs 1. Sending (0, 2) and (0.5, 2) to the diagonal would cost 1 and 0.75, and crossing the pairs costs at least 2.5. Hence d_B(D, D') = 1, \qquad W_1(D, D') = 1.5 .

6. ★ Computational: noise and prominence. Rerun Code 4.2 with noise levels \sigma = 0.02, 0.05, 0.1, 0.2 (same seed). Record the number of local maxima and the three largest finite persistences, and compare the largest noise persistence with the bound of Theorem 4.10.

   Hints and expected result: the numbers of maxima are 521, 721, 800 and 819. The three largest finite persistences are (0.705, 0.583, 0.065), (0.734, 0.640, 0.172), (0.810, 0.747, 0.352) and (1.029, 0.974, 0.712).

   The noise has sup-norm 3.90\,\sigma for this seed, so Theorem 4.10, applied to the clean and noisy images, bounds the noise persistence by 7.8\,\sigma: 0.156, 0.39, 0.78 and 1.56, all respected. At \sigma = 0.2 the third bump is barely separated from the noise (ratio 1.37, against 3.7 at \sigma = 0.05).

7. ★ Distance to a measure by hand. On the line, take X = \{-1, 0, 1, 10\} and k = 2. (a) Compute d_X and the empirical DTM at the points 0, 5 and 10. (b) Move the outlier from 10 to 20, and compare the changes of d_X and of the DTM at 20 with the bound L / \sqrt k of Section 7.4.

   Solution: (a) d_X(0) = 0, d_X(5) = 4, d_X(10) = 0. The DTM is \sqrt{(0 + 1)/2} = 0.707 at 0, \sqrt{(16 + 25)/2} = 4.528 at 5, and \sqrt{(0 + 81)/2} = 6.364 at 10: the isolated point is far for the DTM.

   2. At x = 20, d_X drops from 10 to 0. The DTM goes from \sqrt{(100 + 361)/2} = 15.18 to \sqrt{361/2} = 13.43, a change of 1.75, below L / \sqrt k = 10 / \sqrt 2 = 7.07.

8. ★ Landscapes by hand. For D = \{(0, 4), (1, 3)\} and D' = \{(0, 4)\}, compute \lambda_1 and \lambda_2 of both diagrams and check Theorem 8.2.

   Solution: the tent of (1, 3) lies below the tent of (0, 4), so \lambda_1 = \Lambda_{(0, 4)} for both diagrams, with peak 2 at t = 2.

   Then \lambda_2(D) = \Lambda_{(1, 3)}, with peak 1, and \lambda_2(D') = 0. So \| \lambda_2(D) - \lambda_2(D') \|_\infty = 1 = d_B(D, D'): the bound of Theorem 8.2 is attained.

9. ★★ Signs. Prove that \partial \circ \partial = 0 for the signed boundary of Section 2.4.

   Solution: let i < j, and follow the face obtained by removing v_i and v_j.

   • Removing v_j first, with sign (-1)^j, and then v_i, still in position i, gives the sign (-1)^{i + j}.
   • Removing v_i first, with sign (-1)^i, moves v_j to position j - 1; removing it then gives (-1)^{i + j - 1} in total.

   The two terms have opposite signs and cancel.

10. ★★ Loops of a graph. Prove that a graph with V vertices, E edges and C connected components has \beta_1 = E - V + C over any field, and use it, after ch. 5, to check the loop counts of Code 5.1.

    Solution: a graph has no 2-simplices, so Theorem 2.7 gives V - E = \beta_0 - \beta_1, and \beta_0 = C by Proposition 2.6. Hence \beta_1 = E - V + C. The circle gives 14 - 14 + 1 = 1 and the figure eight 13 - 12 + 1 = 2.

11. ★★ The double torus. A closed orientable surface of genus 2 stands vertically, and f is its height: one minimum, four saddles, one maximum. (a) Check Theorem 3.4. (b) Count the nodes, arcs and loops of the Reeb graph. (c) After a small perturbation that makes f Morse–Smale, count the arcs and cells of its Morse–Smale complex and check Euler’s formula.

    Solution: (a) 1 - 4 + 1 = -2 = 2 - 2g. (b) Six nodes: the extrema have degree 1 and each saddle degree 3, since a level circle splits or two merge. The degrees sum to 2 + 12 = 14, so there are 7 arcs and 7 - 6 + 1 = 2 loops, as Theorem 5.2 predicts.

    3. By Proposition 6.3, 4 \times 4 = 16 arcs and 8 quadrangles, and 6 - 16 + 8 = -2.

12. ★★ Rips stability on the line. Let X = \{0, 1, 3\} and Y = \{0, 1.2, 3\}. Compute the Rips diagrams in degrees 0 and 1 (with r = diameter / 2), their bottleneck distance, and d_H(X, Y).

    Solution: the edges of X have lengths 1, 2 and 3, so \mathrm{Dgm}_0 R(X) = \{(0, 0.5), (0, 1), (0, \infty)\}. For Y, the lengths 1.2, 1.8 and 3 give \{(0, 0.6), (0, 0.9), (0, \infty)\}. In degree 1 both diagrams are empty: the triangle enters with its longest edge.

    Hence d_B = 0.1 \leq d_H(X, Y) = 0.2, as Theorem 7.5 requires.

13. ★★ Computational: Mapper and its parameters. On the figure eight of Code 5.1, (a) vary the number of intervals in \{4, 8, 16\} and the overlap in \{0.1, 0.3, 0.5\}; (b) with 8 intervals and overlap 0.3, vary the clustering scale in \{0.05, 0.1, 0.25, 0.5\}. Count nodes, components and loops.

    Hints and expected result: (a) With 4 intervals, the graph is a path of 4 nodes, with no loop, for every overlap: each interval is too wide to separate the branches. With 8 intervals, there are 2 loops for every overlap (12, 12 and 10 nodes).

    With 16 intervals, overlap 0.1 gives 26 nodes, 22 edges, 4 components and no loop: some overlaps are too thin to contain a point of each branch, and the one at the crossing is empty. Overlaps 0.3 and 0.5 give 26 nodes, 27 edges and 2 loops.

    2. Scale 0.05 shatters the preimages into 32 nodes in 14 components, with no loop. At scale 0.1, six points of the upper right branch split from the central cluster, which opens the right lobe: 13 nodes, one component, a single loop. Scales 0.25 and 0.5 give 12 and 10 nodes, one component and 2 loops.

14. ★★ H_0 of Rips filtrations and spanning trees. Prove that the finite intervals of \mathrm{Dgm}_0 R(X) are [0, \ell_e / 2), where \ell_e runs over the edge lengths of a minimum spanning tree of X.

    Solution: all vertices enter at r = 0, and the edges enter in increasing order of length. The sweep of Proposition 4.3 is then Kruskal’s algorithm: an edge kills a bar exactly when it joins two components, that is, when Kruskal keeps it in the tree. It dies at half its length. The merge tree of this sweep is the single-linkage dendrogram.

15. ★★★ Computational: the torus. Sample 400 points \big((2 + \cos u) \cos v, (2 + \cos u) \sin v, \sin u\big) with u, v uniform on [0, 2\pi) (seed 0). Build the Rips complex up to r = 1 in dimension 2, from the triangles of common neighbours, and compute \mathrm{Dgm}_1 with the cohomology and clearing of Code 9.1. Explain the result.

    Hints and expected result: the complex has 407,279 simplices, and the computation takes a few seconds in plain Python. Two intervals stand out, [0.270, 0.887) and [0.272, 0.874), of persistence 0.617 and 0.602; the next one has persistence 0.370.

    They are the two generators of H_1 of the torus, \beta_1 = 2: a meridian around the tube and a longitude around the hole. Both die near \sqrt 3 / 2 = 0.866: the tube and the inner equator are circles of radius 1, and the Rips complex fills an inscribed equilateral triangle of such a circle at that scale (Code 1.1).

16. ★★★ Open-ended: a topological loss for thin structures. Generate synthetic images of a thin curve (a random smooth curve drawn with a width of 2 pixels) with noise and blur. Predict a probability map with a small convolutional network, trained (a) with a pixelwise cross-entropy, (b) with an added term that penalizes the persistence of every component but the largest one in \mathrm{Dgm}_0 of the superlevel sets of the prediction (Proposition 8.5). Compare the Dice score and the error on the number of components.

    Hints and expected result: with cross-entropy alone, thin parts break, and the predicted mask has several components on many images. The topological term should bring the error on \beta_0 close to zero on most images, with a Dice score close to that of cross-entropy alone. Hu, Li, Samaras and Chen (2019) report this pattern on real data: much lower Betti errors at comparable pixel-wise accuracy. Compute the pairs with the union-find of Code 8.2 extended to the 4-neighbour grid.

Projects and assessment

master-mva.com announces a short lesson exam and a research seminar, with no weights displayed; the Université Paris-Saclay course detail describes the assessment as “Evaluation Continue non Intégrale” (partial continuous assessment). The subjects below are proposals for the research seminar: each ends with a talk, most of them on a research paper and its reimplementation. The written exam is a proposal in the spirit of the short lesson exam, not an official subject.

Solo projects

1. [Solo · ≈ 20 h] Paper study: the stability of persistence diagrams. Reproduce the stability theorem of Cohen-Steiner, Edelsbrunner and Harer (2007) through the vineyards of Cohen-Steiner, Edelsbrunner and Morozov (2006).
   • Work: lower-star persistence of 32 \times 32 images with Algorithm 4.4; follow the linear homotopy f_t at 200 values of t and draw the vineyard; test the bound of Theorem 4.10 on 100 random pairs of images.
   • Success criterion: the bound is never violated, the vineyard figure shows the transpositions, and the report completes the proof sketch of Theorem 4.10 at the swaps.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 4 and 10.
2. [Solo · ≈ 20 h] Paper study: persistence-based clustering. Reproduce the main claims of Chazal, Guibas, Oudot and Skraba (2013) on synthetic data.
   • Work: extend Code 7.2 to a \delta-neighbourhood graph and a DTM-based density; data from the scikit-learn generators (moons, anisotropic blobs, varied densities); baselines k-means, DBSCAN and mean shift.
   • Success criterion: an adjusted Rand index above 0.95 on the two moons with \tau read in the largest gap of the diagram, and a table of the sensitivity to k and to the noise level.
   • Deliverables: code, report, talk. Chapters 4 and 7.
3. [Solo · ≈ 20 h] Paper study: persistence images. Reproduce the classification experiment of Adams and coauthors (2017) on the linked twist map, a synthetic dynamical system on the unit square.
   • Work: orbits of x_{n+1} = x_n + r y_n (1 - y_n) and y_{n+1} = y_n + r x_{n+1} (1 - x_{n+1}) modulo 1, for five values of r; H_0 and H_1 Rips diagrams of subsampled orbits; persistence images, landscapes and the sliced Wasserstein kernel.
   • Success criterion: the accuracy reported in the paper for persistence images is reproduced within 5 points, or the gap is explained (chance is 20 %); a table compares the three representations and the effect of the resolution and of the weight function.
   • Deliverables: code, report, talk. Chapters 7 and 8.

Group projects

1. [Group of 2–3 · ≈ 30 h each] A small topological toolkit for 3D scalar fields. Build, in the spirit of TTK, a pipeline from a 3D grid to a segmentation.
   • Data: synthetic fields on grids from 32^3 to 128^3: sums of Gaussians with noise, and the Marschner–Lobb test signal (Marschner and Lobb, 1994).
   • Split: one member computes the critical points of the 3D Freudenthal triangulation and checks Theorem 3.6; one computes the lower-star persistence and the simplification; one computes merge trees and the segmentation by simplified descending manifolds (Code 6.1 in 3D).
   • Success criterion: the number of Gaussians is recovered at every noise level below a threshold to be measured; the adjusted Rand index against the ground truth; running times that scale nearly linearly with the grid size.
   • Deliverables: code, a 6–8 page report, a 15-minute talk. Chapters 3 to 6 and 9.
2. [Group of 2–3 · ≈ 30 h each] Paper study: optimizing persistent homology based functions. Reproduce the applications of Carrière, Chazal, Glisse, Ike, Kannan and Umeda (2021) with the differentiable persistence of Proposition 8.5.
   • Split: one member optimizes point clouds to create or remove loops; one denoises 2D images topologically (Code 8.2 on grids); one adds a topological regularizer to a small classifier on the scikit-learn digits.
   • Success criterion: the target diagrams are reached on the synthetic tasks, and the effect of the regularizer on test accuracy is measured over 5 seeds.
   • Deliverables: code, report, talk. Chapters 4, 8 and 9.
3. [Group of 2–3 · ≈ 30 h each] Mapper with statistical guarantees. Reproduce the parameter selection and the bootstrap of Carrière, Michel and Oudot (2018).
   • Data: synthetic shapes with known Reeb graphs (circle, figure eight, a Y with flares), and the scikit-learn digits with a PCA, an eccentricity and a density filter.
   • Split: one member implements Mapper with automatic resolution and gain; one the bootstrap stability of its features; one the comparison with the Reeb graphs and the interpretation on the digits.
   • Success criterion: the correct numbers of loops and branches in at least 90 % of the bootstrap samples of the synthetic shapes.
   • Deliverables: code, report, talk. Chapters 5 and 7.

Written exam (proposal)

Format: 1 hour, no documents, no electronic devices, four exercises.

1. Question (4 points). Define the persistence diagram of a filtration and the bottleneck distance. State the stability theorem, and explain why it justifies thresholding persistence to separate features from noise in an image.

   Expected: Definitions 4.6 and 4.7, Theorem 4.10. Noise of sup-norm \delta moves each point of the clean diagram by at most \delta, and every extra point lies within \delta of the diagonal, so it has persistence at most 2\delta.

   A clean feature of persistence above 4\delta keeps a persistence above 2\delta, so a threshold in a gap above 2\delta keeps exactly these features.

2. Question (6 points). (a) Compute the Betti numbers of the boundary of a tetrahedron from the ranks of its boundary matrices. (b) Compute \mathrm{Dgm}_0 of the sublevel sets of the values (4, 1, 3, 0, 5, 2, 6) on a path. (c) Compute its bottleneck distance to \{(1, 3.5), (0, \infty)\}.

   Expected: (a) n = (4, 6, 4), \operatorname{rank} \partial_1 = 3, \operatorname{rank} \partial_2 = 3, so \beta = (1, 0, 1). (b) \{(1, 3), (2, 5), (0, \infty)\}, by the elder rule at the values 3 and 5.

   3. Matching (1, 3) with (1, 3.5) and sending (2, 5) to the diagonal costs \max(0.5, 1.5); the other matchings also cost at least 1.5. So d_B = 1.5.

3. Question (5 points). Prove Banchoff’s formula (Theorem 3.6), and deduce that on a triangulated surface the numbers of minima, saddles counted with multiplicity and maxima satisfy m_0 - m_1 + m_2 = \chi.

   Expected: the grouping of simplices by highest vertex, with the bijection between the simplices of highest vertex v and \{v\} \cup \mathrm{Lk}^-(v); then the Euler characteristics of the lower links of Corollary 3.7.

4. Question (5 points). Prove that \check C(X, r) \subseteq R(X, r) \subseteq \check C(X, 2r) in any metric space. State and prove the sandwich lemma, and deduce how to read the homology of an offset from two Rips complexes.

   Expected: the proof of Theorem 1.7 (a common point gives pairwise distances at most 2r; any vertex is a common point at radius 2r), Lemma 7.3 and Proposition 7.4, with the condition that the inclusions of offsets induce isomorphisms.

Grading rubric (proposal)

References

Official references (master-mva.com)

• [1] H. Edelsbrunner and J. Harer, Computational Topology: An Introduction, American Mathematical Society, 2010.
• [2] J. Milnor, Morse Theory, Annals of Mathematics Studies 51, Princeton University Press, 1963.
• [3] A. T. Fomenko and T. L. Kunii, Topological Modeling for Visualization, Springer, 1997.
• [4] V. Pascucci, X. Tricoche, H. Hagen and J. Tierny (eds.), Topological Methods in Data Analysis and Visualization: Theory, Algorithms, and Applications, Springer, 2010.
• [5] S. Y. Oudot, Persistence Theory: From Quiver Representations to Data Analysis, Mathematical Surveys and Monographs 209, American Mathematical Society, 2015.
• [6] F. Chazal, “High-dimensional topological data analysis”, in J. E. Goodman, J. O’Rourke and C. D. Tóth (eds.), Handbook of Discrete and Computational Geometry, 3rd edition (listed as “to appear” on master-mva.com; published by CRC Press, 2017).
• [7] Gudhi — Geometry Understanding in Higher Dimensions, open-source C++ and Python library, gudhi.inria.fr; see C. Maria, J.-D. Boissonnat, M. Glisse and M. Yvinec, “The Gudhi library: simplicial complexes and persistent homology”, International Congress on Mathematical Software (ICMS), Lecture Notes in Computer Science 8592, Springer, 2014.
• [8] TTK — The Topology ToolKit, topology-tool-kit.github.io; see J. Tierny, G. Favelier, J. A. Levine, C. Gueunet and M. Michaux, “The Topology ToolKit”, IEEE Transactions on Visualization and Computer Graphics 24(1), 2018.

Textbooks and surveys

• A. Hatcher, Algebraic Topology, Cambridge University Press, 2002.
• J. R. Munkres, Elements of Algebraic Topology, Addison-Wesley, 1984.
• A. Zomorodian, Topology for Computing, Cambridge University Press, 2005.
• J.-D. Boissonnat, F. Chazal and M. Yvinec, Geometric and Topological Inference, Cambridge University Press, 2018.
• T. K. Dey and Y. Wang, Computational Topology for Data Analysis, Cambridge University Press, 2022.
• F. Chazal, V. de Silva, M. Glisse and S. Oudot, The Structure and Stability of Persistence Modules, SpringerBriefs in Mathematics, Springer, 2016.
• F. Chazal and B. Michel, “An introduction to topological data analysis: fundamental and practical aspects for data scientists”, Frontiers in Artificial Intelligence 4, 2021.
• G. Carlsson, “Topology and data”, Bulletin of the American Mathematical Society 46(2), 2009.
• R. Ghrist, “Barcodes: the persistent topology of data”, Bulletin of the American Mathematical Society 45(1), 2008.
• A. Björner, “Topological methods”, in R. Graham, M. Grötschel and L. Lovász (eds.), Handbook of Combinatorics, vol. 2, Elsevier, 1995.

Complexes, homology and Morse theory (ch. 0–3)

• R. E. Tarjan, “Efficiency of a good but not linear set union algorithm”, Journal of the ACM 22(2), 1975.
• K. Borsuk, “On the imbedding of systems of compacta in simplicial complexes”, Fundamenta Mathematicae 35, 1948.
• H. W. E. Jung, “Über die kleinste Kugel, die eine räumliche Figur einschließt”, Journal für die reine und angewandte Mathematik 123, 1901.
• H. Edelsbrunner and E. P. Mücke, “Three-dimensional alpha shapes”, ACM Transactions on Graphics 13(1), 1994.
• H. Edelsbrunner and E. P. Mücke, “Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms”, ACM Transactions on Graphics 9(1), 1990.
• T. F. Banchoff, “Critical points and curvature for embedded polyhedra”, Journal of Differential Geometry 1, 1967.
• R. Forman, “Morse theory for cell complexes”, Advances in Mathematics 134(1), 1998.
• R. Forman, “A user’s guide to discrete Morse theory”, Séminaire Lotharingien de Combinatoire 48, 2002.
• V. Robins, P. J. Wood and A. P. Sheppard, “Theory and algorithms for constructing discrete Morse complexes from grayscale digital images”, IEEE Transactions on Pattern Analysis and Machine Intelligence 33(8), 2011.

Persistence (ch. 4 and 10)

• H. Edelsbrunner, D. Letscher and A. Zomorodian, “Topological persistence and simplification”, Discrete & Computational Geometry 28, 2002.
• A. Zomorodian and G. Carlsson, “Computing persistent homology”, Discrete & Computational Geometry 33, 2005.
• W. Crawley-Boevey, “Decomposition of pointwise finite-dimensional persistence modules”, Journal of Algebra and its Applications 14(5), 2015.
• D. Cohen-Steiner, H. Edelsbrunner and J. Harer, “Stability of persistence diagrams”, Discrete & Computational Geometry 37, 2007.
• D. Cohen-Steiner, H. Edelsbrunner and D. Morozov, “Vines and vineyards by updating persistence in linear time”, Symposium on Computational Geometry (SoCG), 2006.
• F. Chazal, D. Cohen-Steiner, M. Glisse, L. J. Guibas and S. Y. Oudot, “Proximity of persistence modules and their diagrams”, Symposium on Computational Geometry (SoCG), 2009.
• M. Lesnick, “The theory of the interleaving distance on multidimensional persistence modules”, Foundations of Computational Mathematics 15, 2015.
• U. Bauer and M. Lesnick, “Induced matchings and the algebraic stability of persistence barcodes”, Journal of Computational Geometry 6(2), 2015.
• H. Edelsbrunner, D. Morozov and V. Pascucci, “Persistence-sensitive simplification of functions on 2-manifolds”, Symposium on Computational Geometry (SoCG), 2006.
• U. Bauer, C. Lange and M. Wardetzky, “Optimal topological simplification of discrete functions on surfaces”, Discrete & Computational Geometry 47, 2012.
• J. Tierny and V. Pascucci, “Generalized topological simplification of scalar fields on surfaces”, IEEE Transactions on Visualization and Computer Graphics 18(12), 2012.
• J. E. Hopcroft and R. M. Karp, “An n^{5/2} algorithm for maximum matchings in bipartite graphs”, SIAM Journal on Computing 2(4), 1973.

Reeb graphs, Mapper and Morse–Smale complexes (ch. 5–6)

• G. Reeb, “Sur les points singuliers d’une forme de Pfaff complètement intégrable ou d’une fonction numérique”, Comptes Rendus de l’Académie des Sciences 222, 1946.
• K. Cole-McLaughlin, H. Edelsbrunner, J. Harer, V. Natarajan and V. Pascucci, “Loops in Reeb graphs of 2-manifolds”, Discrete & Computational Geometry 32, 2004.
• H. Carr, J. Snoeyink and U. Axen, “Computing contour trees in all dimensions”, Computational Geometry 24(2), 2003.
• V. Pascucci, K. Cole-McLaughlin and G. Scorzelli, “Multi-resolution computation and presentation of contour trees”, IASTED Conference on Visualization, Imaging and Image Processing, 2004.
• G. Singh, F. Mémoli and G. Carlsson, “Topological methods for the analysis of high dimensional data sets and 3D object recognition”, Eurographics Symposium on Point-Based Graphics, 2007.
• M. Carrière and S. Oudot, “Structure and stability of the one-dimensional Mapper”, Foundations of Computational Mathematics 18, 2018.
• M. Carrière, B. Michel and S. Oudot, “Statistical analysis and parameter selection for Mapper”, Journal of Machine Learning Research 19, 2018.
• E. Munch and B. Wang, “Convergence between categorical representations of Reeb space and Mapper”, Symposium on Computational Geometry (SoCG), 2016.
• M. Nicolau, A. J. Levine and G. Carlsson, “Topology based data analysis identifies a subgroup of breast cancers with a unique mutational profile and excellent survival”, Proceedings of the National Academy of Sciences 108(17), 2011.
• H. Edelsbrunner, J. Harer and A. Zomorodian, “Hierarchical Morse–Smale complexes for piecewise linear 2-manifolds”, Discrete & Computational Geometry 30, 2003.
• A. Gyulassy, P.-T. Bremer, B. Hamann and V. Pascucci, “A practical approach to Morse–Smale complex computation: scalability and generality”, IEEE Transactions on Visualization and Computer Graphics 14(6), 2008.
• N. Shivashankar and V. Natarajan, “Parallel computation of 3D Morse–Smale complexes”, Computer Graphics Forum 31(3), 2012.
• L. Vincent and P. Soille, “Watersheds in digital spaces: an efficient algorithm based on immersion simulations”, IEEE Transactions on Pattern Analysis and Machine Intelligence 13(6), 1991.

Inference and clustering (ch. 7)

• H. Federer, “Curvature measures”, Transactions of the American Mathematical Society 93, 1959.
• P. Niyogi, S. Smale and S. Weinberger, “Finding the homology of submanifolds with high confidence from random samples”, Discrete & Computational Geometry 39, 2008.
• F. Chazal, D. Cohen-Steiner and A. Lieutier, “A sampling theory for compact sets in Euclidean space”, Discrete & Computational Geometry 41, 2009.
• F. Chazal and S. Oudot, “Towards persistence-based reconstruction in Euclidean spaces”, Symposium on Computational Geometry (SoCG), 2008.
• F. Chazal, D. Cohen-Steiner, L. J. Guibas, F. Mémoli and S. Oudot, “Gromov–Hausdorff stable signatures for shapes using persistence”, Computer Graphics Forum 28(5), 2009.
• F. Chazal, V. de Silva and S. Oudot, “Persistence stability for geometric complexes”, Geometriae Dedicata 173, 2014.
• F. Chazal, D. Cohen-Steiner and Q. Mérigot, “Geometric inference for probability measures”, Foundations of Computational Mathematics 11, 2011.
• F. Chazal, M. Glisse, C. Labruère and B. Michel, “Convergence rates for persistence diagram estimation in topological data analysis”, Journal of Machine Learning Research 16, 2015.
• B. T. Fasy, F. Lecci, A. Rinaldo, L. Wasserman, S. Balakrishnan and A. Singh, “Confidence sets for persistence diagrams”, Annals of Statistics 42(6), 2014.
• F. Chazal, L. J. Guibas, S. Y. Oudot and P. Skraba, “Persistence-based clustering in Riemannian manifolds”, Journal of the ACM 60(6), 2013.

Signatures and learning (ch. 8)

• Y. Mileyko, S. Mukherjee and J. Harer, “Probability measures on the space of persistence diagrams”, Inverse Problems 27(12), 2011.
• K. Turner, Y. Mileyko, S. Mukherjee and J. Harer, “Fréchet means for distributions of persistence diagrams”, Discrete & Computational Geometry 52, 2014.
• M. Carrière and U. Bauer, “On the metric distortion of embedding persistence diagrams into separable Hilbert spaces”, Symposium on Computational Geometry (SoCG), 2019.
• P. Bubenik, “Statistical topological data analysis using persistence landscapes”, Journal of Machine Learning Research 16, 2015.
• F. Chazal, B. T. Fasy, F. Lecci, A. Rinaldo and L. Wasserman, “Stochastic convergence of persistence landscapes and silhouettes”, Journal of Computational Geometry 6(2), 2015.
• H. Adams, T. Emerson, M. Kirby, R. Neville, C. Peterson, P. Shipman, S. Chepushtanova, E. Hanson, F. Motta and L. Ziegelmeier, “Persistence images: a stable vector representation of persistent homology”, Journal of Machine Learning Research 18, 2017.
• J. Reininghaus, S. Huber, U. Bauer and R. Kwitt, “A stable multi-scale kernel for topological machine learning”, IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2015.
• M. Carrière, M. Cuturi and S. Oudot, “Sliced Wasserstein kernel for persistence diagrams”, International Conference on Machine Learning (ICML), 2017.
• M. Carrière, F. Chazal, Y. Ike, T. Lacombe, M. Royer and Y. Umeda, “PersLay: a neural network layer for persistence diagrams and new graph topological signatures”, AISTATS, 2020.
• M. Carrière, F. Chazal, M. Glisse, Y. Ike, H. Kannan and Y. Umeda, “Optimizing persistent homology based functions”, International Conference on Machine Learning (ICML), 2021.
• R. Brüel-Gabrielsson, B. J. Nelson, A. Dwaraknath and P. Skraba, “A topology layer for machine learning”, AISTATS, 2020.
• X. Hu, F. Li, D. Samaras and C. Chen, “Topology-preserving deep image segmentation”, Advances in Neural Information Processing Systems (NeurIPS), 2019.
• M. Moor, M. Horn, B. Rieck and K. Borgwardt, “Topological autoencoders”, International Conference on Machine Learning (ICML), 2020.

Implementations (ch. 9)

• J.-D. Boissonnat and C. Maria, “The simplex tree: an efficient data structure for general simplicial complexes”, Algorithmica 70, 2014.
• J.-D. Boissonnat and S. Pritam, “Edge collapse and persistence of flag complexes”, Symposium on Computational Geometry (SoCG), 2020.
• C. Chen and M. Kerber, “Persistent homology computation with a twist”, European Workshop on Computational Geometry (EuroCG), 2011.
• V. de Silva, D. Morozov and M. Vejdemo-Johansson, “Dualities in persistent (co)homology”, Inverse Problems 27(12), 2011.
• U. Bauer, “Ripser: efficient computation of Vietoris–Rips persistence barcodes”, Journal of Applied and Computational Topology 5, 2021.
• U. Bauer, M. Kerber, J. Reininghaus and H. Wagner, “PHAT — persistent homology algorithms toolbox”, Journal of Symbolic Computation 78, 2017.

Research directions and projects (ch. 11)

• G. Carlsson and A. Zomorodian, “The theory of multidimensional persistence”, Discrete & Computational Geometry 42, 2009.
• M. Lesnick and M. Wright, “Interactive visualization of 2-D persistence modules”, arXiv preprint, 2015.
• M. Carrière and A. Blumberg, “Multiparameter persistence image for topological machine learning”, Advances in Neural Information Processing Systems (NeurIPS), 2020.
• C. Hofer, F. Graf, B. Rieck, M. Niethammer and R. Kwitt, “Graph filtration learning”, International Conference on Machine Learning (ICML), 2020.
• M. Horn, E. De Brouwer, M. Moor, Y. Moreau, B. Rieck and K. Borgwardt, “Topological graph neural networks”, International Conference on Learning Representations (ICLR), 2022.
• J. Vidal, J. Budin and J. Tierny, “Progressive Wasserstein barycenters of persistence diagrams”, IEEE Transactions on Visualization and Computer Graphics 26(1), 2020.
• M. Pont, J. Vidal, J. Delon and J. Tierny, “Wasserstein distances, geodesics and barycenters of merge trees”, IEEE Transactions on Visualization and Computer Graphics 28(1), 2022.
• M. Soler, M. Plainchault, B. Conche and J. Tierny, “Topologically controlled lossy compression”, IEEE Pacific Visualization Symposium (PacificVis), 2018.
• K. Turner, S. Mukherjee and D. M. Boyer, “Persistent homology transform for modeling shapes and surfaces”, Information and Inference 3(4), 2014.
• Y. Hiraoka, T. Nakamura, A. Hirata, E. G. Escolar, K. Matsue and Y. Nishiura, “Hierarchical structures of amorphous solids characterized by persistent homology”, Proceedings of the National Academy of Sciences 113(26), 2016.
• T. Sousbie, “The persistent cosmic web and its filamentary structure — I. Theory and implementation”, Monthly Notices of the Royal Astronomical Society 414(1), 2011.
• S. R. Marschner and R. J. Lobb, “An evaluation of reconstruction filters for volume rendering”, IEEE Visualization, 1994.