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

II. BACKGROUND

• Persistent homology tracks the appearance and disappearance of features across a filtration and summarizes them as barcodes or persistence diagrams.A filtration is a nested sequence of subcomplexes; in image data, the paper focuses on cubical complexes induced by 2D and 3D images.Source: §I · p. 1, §II · pp. 2–3

I. INTRODUCTION

• Persistent homology remains a training bottleneck, and transferring intermediate tensors between GPU and CPU adds PCIe and synchronization overhead.The paper motivates a GPU pipeline that minimizes host-device transfers while keeping most computation on the GPU.Source: §I · p. 1
• TopoGPU combines stratification-aware Morse matching that preserves persistent homology under streaming with parallel algorithms for deriving Morse boundaries.The work also integrates TopoGPU into a topology-preserving deep-learning framework, where persistent-homology computation is no longer the dominant training cost.Source: §I · p. 2
• 53.24× average and 198.01× maximum end-to-end speedups over Cubical Ripser were achieved across evaluated benchmarks.TopoGPU outperformed Cubical Ripser on every evaluated benchmark.Source: §I · p. 2
• Parallel topological sorting achieves a 46.13× average and up to 54.81× speedup over classic Kahn’s algorithm.The sorting and V-path parity algorithms together enable efficient derivation of Morse-complex boundary relations.Source: §I · p. 2

II. BACKGROUND

• Discrete Morse matching pairs cells to form a matching graph, leaving critical cells whose V-path parity determines Morse-complex boundaries over Z2.The resulting Morse boundary matrix is typically much smaller than the cubical-complex boundary matrix while preserving isomorphic homology groups.Source: §II · p. 3

III. RELATED WORK

• Earlier image-persistence work includes approximation, exact boundary-matrix generation, and Morse-complex methods, while boundary-matrix construction remains a bottleneck.Morse-based approaches can reduce complex size; the paper describes streaming and hierarchical GPU partitioning as its approach to scaling computation.Source: §III · p. 3
• Perseus introduced streaming precomputation by simplifying stored level sets independently before merging them, unlike TopoGPU’s spatial volume decomposition.Cubicle instead partitions large image data into slices for streaming CPU processing.Source: §III · p. 3
• TopoGPU’s computational model runs precomputation from cubical-complex construction through Morse-boundary derivation on the GPU under streaming, leaving boundary-matrix reduction to the CPU.The paper contrasts this with CPU-based non-streaming packages and CPU-streaming Cubicle, which lacks GPU parallelism and floating-point support.Source: §III · p. 4
• Prior GPU work targets Vietoris–Rips computations, cubical complexes with Euler curves, and Morse-Smale path counting; TopoGPU targets V-path counting for image persistence.The paper notes that CuMPerLay computes 0D and 1D persistence for 2D images.Source: §III · p. 4

IV. OUR ALGORITHMS

• TopoGPU streams image chunks across grids of GPU blocks, with each block computing local matchings and Morse boundaries independently and in parallel.After all chunks are processed, boundary relations are finalized in host memory and the CPU sorts critical cells, constructs the boundary matrix, and reduces it.Source: §IV · pp. 4–5
• Stratification-aware matching handles block borders by adding auxiliary critical cells that let per-block local results be stitched without border-crossing alternating paths.These extra cells increase the number of critical cells compared with matching the entire complex, but the resulting boundary matrix remains considerably smaller than the original complex’s.Source: §IV · p. 5
• Morse-boundary computation counts V-paths between critical cells, using parity over Z2 to avoid explicitly enumerating potentially exponentially many paths.A parallel topological order is computed first; a parallel path-parity traversal then derives the boundary relations.Source: §IV · p. 5

A. Preprocessing

• Chunks are padded and loaded into a device buffer, then GPU threads process voxels independently to build lower stars and matching graphs, reusing the buffer for subsequent chunks.Each thread handles one voxel; lower stars partition the cubical complex, so filtration values are derived on the fly from introducing voxels.Source: A. Preprocessing · p. 5, 1) Cubical complex construction · p. 5

2) Stratification-Aware Morse Matching:

• The global matching pairs cells only within the same lower star and stratum, preventing border-crossing V-paths and yielding an acyclic matching.Auxiliary critical cells split paths that would cross block borders; they contribute only zero-persistence pairs.Source: 2) Stratification-Aware Morse Matching · p. 6
• The global discrete Morse complex encodes the same persistent homology as the original, up to zero-persistence pairs.The border blocking lemma ensures a path that reaches a border complex remains in it.Source: 2) Stratification-Aware Morse Matching · p. 6
• Two-voxel collars make adjacent chunks induce consistent matchings on border complexes, with width two required to resolve border cells’ lower stars.Collars extend into neighboring chunks along each cutting direction; their values come from adjacent chunks or are set to +∞ where none exists.Source: 2) Stratification-Aware Morse Matching · p. 7

C. Discrete Morse Boundary Computation

• Directly computing a critical cell’s boundary by enumerating V-paths can be infeasible because their number may grow exponentially with complex size.A boundary contains critical cells connected by an odd number of V-paths.Source: C. Discrete Morse Boundary Computation · p. 7
• The prior two-stage method sorts the matching-graph DAG and counts V-paths sequentially, creating challenges for GPU parallelization.The paper addresses these stages with parallel topological sorting and parallel V-path parity computation.Source: C. Discrete Morse Boundary Computation · p. 7

1) Parallel Topological Sorting:

• The proposed topological sorter assigns one GPU thread per vertex and iteratively identifies predecessors, traces toward sources, and assigns orders for bounded-degree matching DAGs.The method is restricted to Morse matching graphs from 2D/3D image data, with degree four in 2D and six in 3D.Source: 1) Parallel Topological Sorting · p. 7

3 Initialize hist as an empty hop count histogram

• Each vertex thread identifies a predecessor, traces until it reaches a source or an unidentifiable predecessor, and records the source and hop count.Successful traces remove outgoing edges and increment the corresponding histogram bin; failed traces resume in later iterations from the furthest reached vertex.Source: §3 · p. 8
• The algorithm assigns topological indices by ordering each iteration’s traceable vertices by hop count and appending them to the prior order.A prefix-sum over the hop-count histogram assigns index ranges to vertices at each distance while preserving chain sequences.Source: §3 · pp. 8–9
• The algorithm has a general upper bound of O(N) parallel time, while its iteration cost depends on predecessor-tracing depth.The bound assumes sufficient thread parallelism and counts synchronization and atomic updates in the standard parallel-time model.Source: §3 · p. 9
• Large trace depth can still yield logarithmic time when handled in one iteration, while shallow traces can take constant time if completed in two iterations.The single-source large-depth regime has D1 = Θ(N) and T = Θ(log N); the bipartite constant-depth regime has r = 2, Dt = Θ(1), and T = Θ(1).Source: §3 · p. 9

2) Parallel Path Parity Computation

• Each GPU thread tracks parity for batches of up to 32 critical cells, packing their path parities into bits updated as vertices are visited in topological order.One 32-bit word stores the parity for 32 critical-cell-to-vertex pairs.Source: 2) Parallel Path Parity Computation · p. 10
• A critical cell belongs to another critical cell’s boundary exactly when the number of V-paths between them is odd.The parallel algorithm replaces path counting with path-parity computation while traversing in topological order.Source: 2) Parallel Path Parity Computation · pp. 9–10
• O(|V|+|E|) time is achieved when the critical-cell count satisfies c ≤ 32·T, where T is the number of GPU threads.The passage states this condition holds in practice for a 1024-thread block and at most 12,825 cubical-complex cells.Source: 2) Parallel Path Parity Computation · p. 10
• Choosing T as a multiple of 32 gives each warp contiguous 32-word accesses for coalesced reads and writes.A warp comprises 32 GPU threads executing the same instruction simultaneously.Source: 2) Parallel Path Parity Computation · p. 10

D. Boundary Matrix Reduction on CPU

• After all chunks are processed, TopoGPU finalizes the complete set of boundary relations in host RAM, avoiding GPU memory capacity limits.The CPU then performs boundary matrix reduction, as stated in the section context.Source: D. Boundary Matrix Reduction on CPU · p. 10
• Critical cells are ordered by filtration value, and their boundary relations are encoded as a binary matrix over Z2.The matrix size is measured by its number of ones; the passage identifies the twist algorithm as the reduction method.Source: D. Boundary Matrix Reduction on CPU · p. 10
• Construction of the boundary matrix and extraction of persistence pairs are parallelized with CTPL.The reduction implementation is identified as PHAT.Source: D. Boundary Matrix Reduction on CPU · p. 10

V. EXPERIMENTS

• Experiments ran on Windows 11 with an Intel Core i5-13600K, 64 GB RAM, and an NVIDIA RTX 4080 GPU with 16 GB memory.The implementation used CUDA 12.7 and C++17; reported tables give standard deviations over five runs.Source: §V · pp. 10–11
• The evaluation included 28 scientific-visualization volumes, with 22 representative volumes reported in the main paper, and simulated breast phantoms spanning four tissue scenarios.The VICTRE dataset contains five MRI volumes per scenario.Source: §V · p. 10
• Gaussian Random Fields were evaluated at 512^3 resolution across spectral-index values, alongside DIPHA, GUDHI, CRipser, and Cubicle as baselines.The passage notes that Cubicle does not support floating-point inputs, while DIPHA, GUDHI, and CRipser failed on some large or complex inputs.Source: §V · p. 11
• TopoGPU produced results identical to tcripser, the CRipser implementation using the T-construction.The evaluation covered end-to-end pipeline performance and stage-wise runtime and peak host-RAM usage.Source: §V · p. 11

A. Performance Comparisons in End-to-End Persistent Homology Computation

• 198.01× maximum and 53.24× average end-to-end speedups over CRipser accompany TopoGPU’s win over all baselines across every dataset.Relative to the best-performing baseline, TopoGPU reaches up to 71.95× speedup, averaging 11.53×.Source: A. Performance Comparisons in… · p. 11
• 321.46× larger cubical-complex boundary matrices and 8.99× larger partitioned Morse-complex matrices, on average, contrast with TopoGPU’s boundary matrices.DIPHA, GUDHI, and CRipser crashed on some large inputs because of non-streaming designs or oversized boundary matrices exceeding available RAM.Source: A. Performance Comparisons in… · p. 11
• 10.77 MVox/s average throughput, peaking at 32.90 MVox/s on the Aneurysm volume, was reported for TopoGPU.Throughput is measured in millions of voxels processed per second.Source: A. Performance Comparisons in… · p. 12

B. Runtime Breakdown

• 52.39% of runtime went to boundary matrix reduction on average, followed by path parity at 27.96%, I/O at 5.68%, and device-host transfers at 5.48%.Topological sorting and Morse matching were minor stages.Source: B. Runtime Breakdown · p. 12, Table II · p. 12
• Path parity often dominated small and medium volumes, while boundary matrix reduction dominated large and complex volumes.Examples include 78.76% path parity for Shockwave and 61.59% reduction for Pawpawsaurus.Source: B. Runtime Breakdown · p. 12
• Peak host RAM ranged from 0.19 GB on the smallest volume to 54.04 GB on Pawpawsaurus; TopoGPU ranked second-lowest among baselines.Only Cubicle used less RAM, by offloading boundary relations to external storage.Source: B. Runtime Breakdown · p. 12
• Log-runtime and reduction-share association was positive and statistically significant: Pearson r = 0.64 (p = 0.0076) and Spearman ρ = 0.62 (p = 0.010).The passage interprets the association as indicating that matrix reduction is the primary bottleneck in this pipeline.Source: B. Runtime Breakdown · p. 12
• Figure 6’s stacked bars show stage shares and its black curve shows total runtime; path parity dominates at low persistence-pair counts, while reduction dominates at high counts.The yellow remainder comprises initialization, array allocation, and resource release.Source: B. Runtime Breakdown · p. 12, Figure 6 · p. 13

C. Comparative Analysis of Topological Sorting

• TopoSort targets bounded-degree matching graphs and can order vertices whose dependencies are resolved through predecessor tracing, exposing parallelism beyond zero-in-degree frontiers.The subsection evaluates it as a standalone algorithm that outputs a topological ordering.Source: C. Comparative Analysis of Topological… · p. 13
• TopoSort beats CPUKahn by 46.13× on average and up to 54.81×, and GPUKahn by 1.61×–9.80×; 89.0% of vertices are sorted before reaching zero in-degree.Across the 16 benchmark graphs, sorting takes 11–17 iterations, averaging 13.25.Source: C. Comparative Analysis of Topological… · p. 13

D. Comparison with GPU-Based 2D Cubical Persistence

• TopoGPU is 3.3× faster than CuMPerLay on ISIC 2018 and 24.6× faster on the higher-resolution Glaucoma dataset.The comparison uses CuMPerLay in its native 2D setting on two datasets.Source: D. Comparison with GPU-Based 2D… · p. 13

E. Deep Learning Pipeline Integration

• With TopoGPU, the topology-preserving framework reaches loss 0.011 in 6.45 hours, versus CRipser’s 0.014 after 1.71 days; PH takes 19.09% rather than 87.31% of training time.The comparison uses Hu et al.’s framework and reports averages of 48.56 seconds per epoch with TopoGPU and 309.76 seconds with CRipser.Source: E. Deep Learning Pipeline Integration · p. 14

VI. LIMITATIONS AND FUTURE WORK

• TopoGPU’s four stages are cubical-complex construction, stratification-aware Morse matching, GPU Morse-boundary computation, and CPU boundary-matrix reduction.Parallel topological sorting and path-parity computation derive the Morse boundaries.Source: §VI · pp. 14–15
• Boundary-matrix reduction fails when a single volume requires more RAM than available; very small inputs can also cost more to process with TopoGPU than with a lightweight CPU method.The reported failure volumes require more than 201.9 GB of RAM for matrix reduction; overhead can exceed computation on small inputs.Source: §VI · p. 15

I. FULL END-TO-END RUNTIME AND SPEEDUP RESULTS

• Against CRipser, TopoGPU reaches up to 198.01× speedup, averaging 53.63×; versus the fastest successful baseline per dataset, it reaches up to 71.95×, averaging 10.19×.TopoGPU has the lowest end-to-end runtime among successful methods across all datasets.Source: §I · p. 18
• DIPHA, GUDHI, and CRipser fail on several large or complex inputs, while TopoGPU succeeds across the benchmark suite; the cubical-complex boundary matrix averages 321.46× the TopoGPU matrix.The TopoGPU boundary matrix averages 8.99× the size of the non-partitioned Morse-complex matrix.Source: §I · p. 18
• TopoGPU averages 9.97 MVox/s throughput across datasets and reaches 32.90 MVox/s on Aneurysm.Throughput is measured in millions of voxels processed per second.Source: §I · p. 18

II. BOUNDARY MATRIX SIZE COMPARISON WITH NON-PARTITIONED MORSE COMPLEX

• Partitioning makes TopoGPU’s boundary matrix 8.99× larger on average than the non-partitioned Morse matrix, but it remains 321.46× smaller than the cubical-complex matrix.The comparison uses datasets where the non-partitioned baseline is available.Source: §II · p. 18

III. FULL STAGE-WISE RUNTIME BREAKDOWN

• Boundary-matrix reduction accounts for 54.31% of runtime on average, followed by path-parity computation at 25.63% and host-device transfer at 7.02%.Small and medium volumes can instead be path-parity heavy, while large and complex volumes are typically reduction dominated.Source: §III · p. 19
• Peak RAM ranges from 0.19 GB to 54.67 GB across evaluated datasets; TopoGPU has the second-lowest peak RAM usage, behind Cubicle.RAM use scales with the resulting boundary-matrix size.Source: §III · p. 19

IV. RUNTIME SCALING WITH TOPOLOGICAL COMPLEXITY

• On log-log axes, runtime has a near-linear relationship with both critical-cell count and boundary-matrix size, making resulting Morse-complex size a stronger cost predictor than voxel count.The plotted results cover the first 16 datasets; two equal-size volumes can differ in topological complexity and resulting matrix size.Source: §IV · pp. 19–20, Figure S.1 · p. 20

V. EMPIRICAL ITERATION COUNTS OF PARALLEL TOPOLOGICAL SORTING

• Across 16 matching graphs ranging from 9.4 · 10^5 to 2.2 · 10^8 vertices, TopoSort takes 11–17 iterations, averaging 13.25.The observed behavior indicates that the evaluated image-derived matching graphs do not exhibit large iteration counts in practice.Source: §V · pp. 20–21

• Small inputs incur device setup and transfer costs, and one thread per voxel leaves much of the GPU idle; concurrent batches can saturate the device and amortize overhead.TopoGPU’s overheads include CUDA context creation, memory allocation and initialization, and host-device transfers.Source: §VI · pp. 21–22
• At 256 × 256 and larger, TopoGPU beats CRipser even one input at a time, and batching widens the advantage; on constant 1024 × 1024 inputs it is 22× faster.On constant inputs, the nearly trivial Morse complex leaves little work for sequential CPU boundary-matrix reduction.Source: §VI · p. 22
• Batching improves single-input processing on random inputs by 2–4×, including a change from 408.6 ms to 100.9 ms at 1024 × 1024.Batching helps only modestly on constant inputs because GPU overheads dominate their nearly trivial Morse-complex computation.Source: §VI · p. 22
• At 32 × 32, CRipser remains faster than TopoGPU in single-input and batched modes because fixed overheads outweigh parallelism gains.Batching narrows, but does not close, the gap.Source: §VI · p. 22

VII. DEEP LEARNING PIPELINE INTEGRATION

• The integration compares otherwise identical CRipser- and TopoGPU-based training frameworks across three 3D U-Net architectures and the PARSE and MSD Brain datasets.PARSE volumes are resized to 128×128×128; MSD Brain masks combine three tumor sub-regions into one binarized Whole Tumor label.Source: §VII · p. 23
• Runtime measurements average 480 epochs, each computing persistent homology for 96 likelihood maps and 96 segmentation masks of size 32×32×32.The basic configuration uses batch size 4, learning rate 10^-4, and Adam on an Ubuntu 18.04 multi-GPU workstation with NVIDIA GTX 1080 GPUs.Source: §VII · p. 23
• Across datasets and architectures, TopoGPU reduces average epoch time, including PARSE PH time from about 245–272 seconds to 8–9 seconds and MSD Brain PH time from about 135–155 seconds to 13–16 seconds.The main performance gain comes from accelerating persistent-homology computation; Table S.IX reports the Attention 3D U-Net comparison on PARSE and MSD Brain.Source: §VII · p. 23, PARSE · p. 24
• Residual 3D U-Net generally has the highest runtime, followed by Attention 3D U-Net and standard 3D U-Net, while TopoGPU’s relative advantage remains consistent across all three.The paper attributes the runtime ordering to added computations from residual blocks and attention gates, and says acceleration is primarily determined by the persistent-homology backend.Source: §VII · p. 23

PARSE

• For PARSE, Table S.IX compares CRipser and TopoGPU under different configurations using average epoch and persistent-homology time in seconds and peak GPU memory in MB.GPU1 is the network-training device and GPU2 is used by TopoGPU for persistent-homology computation.Source: PARSE · p. 24