GPU-Accelerated Computation of Persistent Homology for Topological Analysis of Image Data

AI summary

Overview

Research area: Topological data analysis (TDA) and topological deep learning (TDL), specifically GPU-accelerated computation of persistent homology for 2D and 3D image data.

Technical level: Advanced. The paper assumes familiarity with algebraic topology concepts (homology groups, filtrations, persistence diagrams, discrete Morse theory), GPU/CUDA programming models, and parallel algorithm design.

Scope: The paper presents TopoGPU, a GPU streaming pipeline that computes persistence diagrams of cubical complexes induced by 2D and 3D grayscale images, leaving only the boundary matrix reduction on the CPU.

Accepted for publication in IEEE Transactions on Pattern Analysis and Machine Intelligence, DOI 10.1109/TPAMI.2026.3741473.

What This Paper Is About

Persistent homology is a cornerstone of topological data analysis that summarizes how topological features such as connected components, loops, and voids appear and disappear across scales, but its computation remains a major bottleneck when it is used inside deep network training. The paper's goal is to move the heavy parts of this computation onto the GPU, avoiding the expensive host-device transfers and synchronization that occur when intermediate tensors such as likelihood maps must be offloaded to the CPU. TopoGPU targets only cubical complexes induced by 2D and 3D images; it explicitly does not target point clouds, Vietoris-Rips complexes, or general simplicial complexes.

Key Contributions

1. TopoGPU, a massively parallel GPU algorithm for persistent homology computation of cubical complexes induced by 2D/3D grayscale images. It is reported as faster than Cubical Ripser (CRipser), a state-of-the-art method, on every evaluated benchmark, with average and maximum end-to-end speedups of 53.24 times and 198.01 times.

2. Stratification-aware discrete Morse matching that preserves the persistent homology of the input under a streaming scheme. It prevents border-crossing paths by introducing auxiliary critical cells, at the cost of additional critical cells compared to matching on the entire cubical complex.

3. A parallel topological sorting algorithm for bounded-degree Morse matching graphs, reported to achieve a 46.13 times average and up to 54.81 times speedup over classic Kahn's algorithm.

4. A parallel V-path parity algorithm on the Morse matching graphs, running in time linear in the graph size under a parallel model with sufficient thread parallelism. Together with the parallel topological sorting algorithm, it enables efficient derivation of boundary relations of the Morse complex.

5. Integration into a topology-preserving deep learning framework, demonstrating that persistent homology computation is no longer the dominant cost during network training.

Main Findings

• End-to-end speedup: 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.24 times and a maximum of 198.01 times.

• Parallel topological sorting speedup: The proposed parallel topological sorting algorithm for bounded-degree Morse matching graphs achieves a 46.13 times average and up to 54.81 times speedup over classic Kahn's algorithm.

• Correctness under streaming: Corollary 1 (Global Correctness) states that the global discrete Morse complex encodes the same persistent homology up to zero-persistence pairs. The paper's Border Blocking Lemma establishes that once a directed path reaches a border complex, it cannot leave, and within a border complex it can only remain in the same-dimensional border stratum or move to lower-dimensional border strata.

• Exponential to linear reduction: Over Z_2 coefficients, boundary relations depend only on the parity (odd/even) of the number of V-paths between critical cells, reducing the problem from exponential-time path enumeration to a linear-time parity computation.

• Smaller boundary matrix: The boundary matrix of the Morse complex is typically much smaller than that of the cubical complex, offering significant computational savings while yielding homology groups isomorphic to those of the original complex. Despite the additional auxiliary critical cells introduced by stratification-aware matching, the resulting boundary matrix remains considerably smaller than that of the original complex.

• Deep learning training cost: Integration into a topology-preserving deep network demonstrates that persistent homology computation is substantially reduced in cost during network training and is no longer the dominant cost.

• Collar width justification: By the cited prior work, determining the lower star of a voxel requires its 8-connected neighbors in 2D or 26-connected neighbors in 3D, so resolving a border cell's containing lower star reaches one voxel across the chunk border and one voxel beyond it, making a collar width of two the minimum required.

Methodology in Plain English

The approach combines streaming data decomposition with discrete Morse theory and GPU parallelism.

Partitioning and streaming. The input image is cut into chunks: 2D images along the y-axis, and 3D images along both the y- and z-axes. A 3D image of size (W, H, D) is cut c_y times along the y-axis and c_z times along the z-axis, yielding chunks of size (w, h, d) = (W, H/c_y, D/c_z). This layout respects C++ row-major storage and maximizes contiguous memory access. Each chunk is padded with a two-voxel-wide collar drawn from adjacent chunks: 2D chunks of dimensions (w, h) become (w, h+4), and 3D chunks of dimensions (w, h, d) become (w, h+4, d+4). Collar values come from adjacent chunks when available, or are set to +infinity if no adjacent chunk exists; no padding is applied along the x-axis. Chunks are loaded with the C++ Thread Pool Library (CTPL) using parallel CPU threads and stored in a single device buffer.

Per-block parallelism. Each chunk is covered by a grid of non-overlapping GPU blocks, excluding the collars. Each GPU block organizes threads in a 2D/3D structure mirroring the spatial layout of pixels/voxels, so each voxel is assigned exactly one GPU thread. Because the lower stars of all voxels form a disjoint partition of the cubical complex, threads execute independently with no dependencies across blocks. Filtration values are not explicitly stored for the cubical complex; they are determined on the fly from introducing voxels during Morse matching. After all blocks in a chunk finish, the device buffer is reused to load subsequent chunks.

Stratification-aware matching. The complex of a block of voxels is split into D+1 strata indexed by dimension. The D-dimensional stratum is the interior stratum; the lower-dimensional strata are border strata. Matchings are computed independently within each block's interior stratum and within its border strata of each dimension, and only among cells in the same lower star and same stratum. Cells in 0-dimensional border strata are always critical because no cells are available for pairing within the same stratum. The union of these matchings forms the global matching. Collar voxels are not given their own blocks or threads; the interior threads nearest the chunk border process them after finishing their own work.

Boundary derivation. With the matching graph built per block, boundary relations of the Morse complex are derived by counting V-paths between critical cells. Instead of enumerating paths (which can be exponential), the method uses parity over Z_2. This is done on the GPU with two algorithms: first a parallel topological sort to give the linear order required for boundary derivation, then a parallel path parity algorithm that traverses that order tracking cumulative path parities between critical cells. The matching graph is stored in a cube-map format. For a 3D chunk of size (w, h, d), the cube-map is an array of size (2w+1, 2h+1, 2d+1), where the dimension of an entry with index i = (i1, i2, i3) is the sum over k = 1 to 3 of (i_k mod 2). Each entry occupies 12 bits, partitioned into six consecutive 2-bit segments representing front, back, left, right, top, and bottom directions, each encoding incoming edge, outgoing edge, or disconnected.

CPU stage. After all blocks across all chunks have been processed and boundary relations finalized in host memory, critical cells are sorted by filtration value, then the boundary matrix is constructed and reduced to compute persistent homology.

Hardware constraints noted. On the NVIDIA RTX 4080 used in the experiments, per-block limits are up to 1,024 threads, 99 KB shared memory, and 64 K 32-bit registers.

Why This Matters

Impact on research. The work addresses a practical obstacle to topological deep learning: persistent homology computation is a bottleneck during network training, standing in sharp contrast to GPU-optimized deep learning workflows that keep tensors and computation resident on the device. By keeping the entire precomputation -- from cubical complex construction through Morse boundary derivation -- on the GPU under a streaming scheme and leaving only the boundary matrix reduction on the CPU, TopoGPU aims to catalyze broader adoption of persistent homology within deep learning pipelines. It is described as the first to do so in this configuration, distinguishing it from GUDHI, DIPHA, and CRipser (CPU-based and non-streaming), Cubicle (streams on the CPU but offers no GPU parallelism and no floating-point support), and CuMPerLay (a CUDA implementation restricted to 2D images via union-find and Alexander duality, which does not extend to 3D volumes).

Real-world applications:

• Medical imaging, such as CT scans, where images are acquired directly and analyzed topologically.
• Physical simulations, where scalar fields from numerical computations form the image-like input.
• Intermediate outputs of convolutional neural networks, notably likelihood maps, which need topological analysis during training.
• Topology-preserving deep learning, where topology-based losses and regularizers constrain networks to respect target topology during training.

Industry relevance. The paper states that because modern GPUs are increasingly powerful and GPU-optimized deep learning workflows dominate practice, a GPU-resident persistent homology pipeline reduces PCIe transfer and synchronization overhead that would otherwise arise from offloading tensors to the host. 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

Future Directions

• Scaling and memory bounds. TopoGPU can handle large datasets via streaming that keeps only one chunk in device memory at a time, and it can process arbitrarily large inputs provided the resulting boundary matrix fits in host memory. Extending beyond that host-memory constraint is a natural open problem.
• Moving boundary matrix reduction to the GPU. The pipeline currently leaves the boundary matrix reduction on the CPU; the paper does not report GPU reduction.
• Beyond cubical complexes. The paper explicitly scopes TopoGPU to cubical complexes induced by 2D and 3D images and states it does not target point clouds, Vietoris-Rips complexes, or general simplicial complexes.
• Chunking and collar tuning. The paper notes the trade-off that stratification-aware matching introduces additional critical cells compared to matching on the entire cubical complex, while the resulting boundary matrix remains considerably smaller than the original. How chunk and block layout parameters affect this trade-off is left open in the provided content.

Target Audience

Researchers and practitioners working at the intersection of topological data analysis and deep learning, particularly those who need persistent homology of 2D/3D image data at scale. It will be most useful to GPU/CUDA developers and parallel algorithm designers interested in discrete Morse theory, to computer vision and medical imaging researchers applying cubical persistence, and to developers integrating topology-preserving losses into neural network training. Readers without a background in algebraic topology, discrete Morse theory, or parallel programming will find the material demanding, since the contributions are algorithmic and the correctness arguments are stated as formal lemmas and corollaries over Z_2 coefficients.

Authors’ abstract

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.