GPU-Accelerated Computation of Persistent Homology for Topological Analysis of Image Data
Authors: Fan Wang, Hubert Wagner, Rezaul Chowdhury, Chao Chen
Organizations: Department of Computer Science, Stony Brook University, Stony Brook, NY 11794 USA · Department of Mathematics, University of Florida, Gainesville, FL 32611 USA · Department of Biomedical Informatics, Stony Brook University, Stony Brook, NY 11794 USA
In recent years, persistent homology has seen rapid adoption in deep learning, yet its computation remains a major bottleneck in network training. This paper introduces TopoGPU, a GPU streaming pipeline that computes persistence diagrams of cubical complexes induced by 2D and 3D images. TopoGPU streams the input image chunk by chunk, processing each chunk with massively parallel GPU kernels on a grid of GPU blocks; the resulting boundary relations are accumulated in host memory, where the CPU performs the boundary matrix reduction. TopoGPU introduces a stratification-aware discrete Morse matching that provably preserves persistent homology under streaming, together with a parallel topological sorting algorithm and a parallel V-path parity algorithm for deriving Morse boundaries on the GPU. TopoGPU outperforms Cubical Ripser, a state-of-the-art method for persistent homology computation, on every benchmark evaluated, achieving an average end-to-end speedup of 53.24x and a maximum of 198.01x. We further integrate TopoGPU into a topology-preserving deep network, demonstrating that it substantially reduces the cost of persistent homology computation during network training. TopoGPU is open source, with pre-built binaries, Google Colab notebooks, and Docker images available at the project's GitHub page: https://github.com/seravee08/GPU-Computation-of-Persistent-Homology-for-Image-Data.
Figures & tables
Fig. 1: (a) Input image; (b) Filtered cubical complex of (a) with solid arrows initially pointing from cells to their proper faces; (c) Matching graph of (b), with dashed arrows indicating Morse matchings, critical cells in red, and red arrows indicating a V -path; (d) Stratification-aware Morse matching after splitting (a) into two blocks, with border region shaded; (e) Persistence diagram encoding boundary relations with a persistence pair ( 10,12 ).
Fig. 2: (a) A 3D input partitioned into chunks; (b) collar-padded chunk tiled by a grid of non-overlapping GPU blocks, with one GPU thread per voxel; (c) front 2 -dimensional border complex of the green GPU block in (b), surrounded by 26 neighboring blocks.
Fig. 3: Examples of 2D local-dependency configurations for vertex v within its 8-connected neighborhood; for clarity, only the right half is shown.
Fig. 4: Illustration of Algorithm 1 for parallel topological sorting. Step I : identify predecessors. Step II : trace from green vertices (those with identified predecessors) backward along red arrows ( pre(⋅) ) to red vertices (source vertices), deleting the outgoing edges of each traceable vertex; for clarity, traceable vertices are also removed from the illustration.
Fig. 5: DAG configurations illustrating three representative regimes for parallel topological sorting: large trace depth, constant depth, and large iteration count; legend as in Figure 4 .
Dataset
Avg. #Persist.
TopoGPU
Cubicle
CRipser
DIPHA
GUDHI
Topo/CRipser
Topo/best
TopoGPU
Pairs
(secs)
(secs)
(secs)
(secs)
(secs)
Speedup
Speedup a
MVox/s b
Silicium
479
0.02 ± 0.00
0.62 ± 0.00
0.11 ± 0.00
0.66 ± 0.00
0.56 ± 0.03
5.50
5.50
5.66
Fuel
126
0.03 ± 0.00
0.65 ± 0.04
0.28 ± 0.00
1.43 ± 0.03
1.33 ± 0.00
9.33
9.33
8.74
Neghip
346
0.03 ± 0.00
0.76 ± 0.02
0.31 ± 0.00
1.48 ± 0.02
1.36 ± 0.03
10.33
10.33
8.74
Tooth
204769
0.59 ± 0.01
2.29 ± 0.03
3.80 ± 0.11
9.90 ± 0.10
11.19 ± 0.17
6.44
3.88
2.64
Hydrogen
7
0.12 ± 0.00
0.89 ± 0.03
3.67 ± 0.07
12.04 ± 0.31
12.16 ± 0.04
30.58
7.42
17.48
TABLE I: End-to-end runtime for persistent homology computation: TopoGPU vs. state-of-the-art baselines.
Dataset
Size
Data Type
BoundMat Size a
Dev ↔ Host
Morse Matching
Topo Sort
Path Parity
Matrix Reduction
I/O
RAM b
VRAM c
(GB)
(GB)
Silicium
98×34×34
uint8
3.24⋅104
0.00%
1.58%
2.40%
22.20%
44.72%
10.17%
0.19
0.34
Fuel
64×64×64
uint8
5.47⋅104
0.00%
1.48%
3.12%
47.76%
28.02%
14.62%
0.19
0.31
Neghip
64×64×64
uint8
9.09⋅104
0.00%
1.39%
3.57%
48.08%
26.78%
10.18%
0.20
0.31
Tooth
103×94×161
uint8
2.80⋅106
0.19%
0.63%
1.25%
15.38%
77.85%
3.54%
0.42
0.53
Hydrogen
128×128×128
uint8
2.72⋅105
0.02%
2.64%
6.94%
55.33%
19.82%
14.05%
0.28
0.91
TABLE II: Stage-wise runtime breakdown and peak host RAM and device memory usage for TopoGPU.
Fig. 6: Per-stage runtime breakdown of TopoGPU. Stacked bars show stage-wise percentages (left y-axis); the black curve shows total runtime of the first 16 volumes (right y-axis). Percentages inside the green bars indicate the share of boundary matrix reduction.
Fig. 7: Runtime curves for CPUKahn, GPUKahn, and TopoSort (left y-axis) with TopoSort-over-GPUKahn speedup curve (right y-axis). Stacked bar charts show the percentage of vertices topologically sorted while d−(⋅)>0 (green; percentage displayed) versus when d−(⋅)=0 (red). Volumes are ordered by increasing vertex count in their matching graphs, with vertex count annotated above each bar whose height corresponds to the count.
Dataset
Resolution / #Images
Avg. time per image (ms)
Speedup
CuMPerLay
TopoGPU
ISIC 2018
512×512 / 1,000
21.30
6.43
3.31 ×
Glaucoma
1690×2004 / 5,335
1128.44
45.87
24.60 ×
TABLE III: End-to-end runtime versus CuMPerLay on 2D image datasets.
Fig. 8: PARSE dataset visualizations with ground truth segmentations in red; subjects 005, 082, 120, and 285 shown from left to right.
Fig. 9: Left y-axis: training loss versus training time. Right y-axis: cumulative persistent homology (PH) computation time versus training time. Both networks trained for 480 epochs.
Persistence-based topological optimization deforms a point cloud X⊂Rd by minimizing objectives of the form L(X)=ℓ(Dgm(X)), where Dgm(X) is a persistence diagram. In practice, optimization is limited by two coupled issues: persistent homology is typically computed on subsamples, and the resulting topological gradients are highly sparse, with only a few anchor points receiving nonzero updates. Motivated by diffeomorphic interpolation, which extends sparse gradients to smooth ambient vector fields via Reproducing Kernel Hilbert Space (RKHS) interpolation, we propose a more scalable pipeline that improves both subsampling and gradient extension. We introduce subsampling via random slicing, a lightweight scheme that promotes iteration-wise geometric coverage and mitigates density bias. We further replace the costly kernel solve with a fast Nadaraya-Watson (NW) Gaussian convolution, producing a globally defined smooth update field at a fraction of the computational cost, while being more suited for topological optimization tasks. We provide theoretical guarantees for NW smoothing, including anchor approximation bounds and global Lipschitz estimates. Experiments in 2D and 3D show that combining random slicing with NW smoothing yields consistent speedups and improved objective values over other baselines on common persistence losses.
Persistent homology (PH) is a frequently used tool for extracting and preserving topological information from image data, particularly in image segmentation, where preservation of topological structures is important. However, despite its general applicability across dimensionality, domains, and target structures, the runtime cost of PH-based methods often makes their practical use infeasible. In this work, we argue that this runtime cost is largely driven by processing information that is unimportant for downstream application (e.g. as optimization objective). We propose sparse cubical filtrations as an alternative foundation for PH computation, reducing subsequent computational costs by factors of up to 100 on real datasets. We show close agreement with the optimization signal of the dense counterpart and empirically evaluate our solution's effectiveness as an optimization objective in realistic training regimes where other PH-based objectives can practically not operate (i.e., 3D data with large patch sizes). We show how our solution improves topological accuracy by up to 80% across six diverse datasets while maintaining pixel- and region-based accuracy.
Alexander H. Berger, Marco Fontana, Daniel Rueckert +3
Weill Cornell Medicine, New York, USA · Technical University of Munich, Munich, Germany · Department of Computing, Imperial College London, UK +3
We present Flash Cubical, a highly efficient computation of cubical persistence on a V-filtration for 2D and 3D images over F2. The implementation is built around three core ideas. First, cubical complexes satisfy properties that allow for the computation of persistence of the highest dimension via union-find and duality. Second, pruning of certain edges allows for a fast and efficient implementation of union-find. Third, the use of a lookup table, which exploits the regularity of cubical complexes to pre-compute local information. This avoids the need to compute local information at run time. To the best of our knowledge, this is the most efficient implementation of cubical persistence with a V-filtration, both in terms of time and memory costs. Although the paper focuses on persistence for V-filtration cubical complexes, the underlying ideas generalise naturally to T-filtrations on cubical complexes and suggest promising directions for other complexes.
Titouan Le Breton, Karol Szustakowski, Marie Piraud
1Helmholtz AI, Helmholtz Munich, Neuherberg, Germany · 2École des Ponts ParisTech, ENS Paris-Saclay · Institute of AI for Health, Helmholtz Munich, Neuherberg, Germany