Distributed Discrete Morse Sandwich: Efficient Computation of Persistence Diagrams for Massive Scalar Data

TL;DR

The paper addresses the scalability challenge of computing persistence diagrams for massive 3D scalar fields by extending the fast Discrete Morse Sandwich (DMS) algorithm to distributed-memory systems, yielding the Distributed Discrete Morse Sandwich (DDMS). It introduces a hybrid MPI+OpenMP implementation with a self-correcting distributed pairing for $\mathcal{D}_0$ and $\mathcal{D}_2$, a distributed extension of PairCriticalSimplices for $\mathcal{D}_1$, computation tokens, and a dedicated communication thread to overlap communication and computation. The approach delivers strong and weak scaling up to 16 nodes (512 cores), achieving up to an 8× speedup over DIPHA and handling a 6 billion-vertex dataset in 174 seconds on 512 cores, with open-source availability for reproducibility. These results establish DDMS as a practical, scalable solution for topological data analysis on massive datasets and lay groundwork for future enhancements in memory efficiency and heterogeneous architectures.

Abstract

The persistence diagram, which describes the topological features of a dataset, is a key descriptor in Topological Data Analysis. The "Discrete Morse Sandwich" (DMS) method has been reported to be the most efficient algorithm for computing persistence diagrams of 3D scalar fields on a single node, using shared-memory parallelism. In this work, we extend DMS to distributed-memory parallelism for the efficient and scalable computation of persistence diagrams for massive datasets across multiple compute nodes. On the one hand, we can leverage the embarrassingly parallel procedure of the first and most time-consuming step of DMS (namely the discrete gradient computation). On the other hand, the efficient distributed computations of the subsequent DMS steps are much more challenging. To address this, we have extensively revised the DMS routines by contributing a new self-correcting distributed pairing algorithm, redesigning key data structures and introducing computation tokens to coordinate distributed computations. We have also introduced a dedicated communication thread to overlap communication and computation. Detailed performance analyses show the scalability of our hybrid MPI+thread approach for strong and weak scaling using up to 16 nodes of 32 cores (512 cores total). Our algorithm outperforms DIPHA, a reference method for the distributed computation of persistence diagrams, with an average speedup of x8 on 512 cores. We show the practical capabilities of our approach by computing the persistence diagram of a public 3D scalar field of 6 billion vertices in 174 seconds on 512 cores. Finally, we provide a usage example of our open-source implementation at https://github.com/eve-le-guillou/DDMS-example.

Figures (17)

• Figure 1: Filtration example for a PL scalar field representing the elevation on a terrain. As simplices are added by increasing scalar value, the topology of the domain changes: in (a), the terrain is composed of two connected components. In (b), only one connected component remains, but a handle has appeared. In (c), the smaller hill at the front has been filled and there is still a handle. In (d), there are two handles. In (e), only one handle is left. Finally, in (f), the whole domain has been swept (one connected component).
• Figure 2: Example of discrete gradient field for a given scalar field (left). The larger spheres represent critical simplices. The light yellow arrows are edge-triangle pairs. The light blue arrows are vertex-edge pairs. Two examples of v-paths are shown in transparent blue (right), going from a critical edge and ending in critical vertices. These form the unstable set of the critical edge.
• Figure 3: Examples of $1$-cycles on a sub-complex of a filtration. The PL scalar field represents the elevation on a terrain. Both the orange and the blue cycles are $1$-cycles: the computation of their boundaries results in summing all the vertices of each cycle twice (as all vertices are faces of exactly two edges of the cycle), yielding a boundary equal to $0$ for each cycle. However, the two cycles are not homologous: it is not possible to go from one to the other by adding a $p$-boundary. The blue cycle can be chosen as the representative of its class, representing the hole at the top of the hill it surrounds.
• Figure 4: Persistence diagram for a clean (left) and noisy (right) dataset. Critical simplices are denoted by spheres (light blue: minima, light yellow: maxima, others: saddles). The persistence of each pair is the height of the bar. $\mathcal{D}_{0}$ corresponds to minimum-saddle pairs and $\mathcal{D}_{2}$ corresponds to saddle-maximum pairs. Salient features (long pairs) can easily be distinguished from the noise (short pairs).
• Figure 5: Overview of the DMS algorithm for the computation of $\mathcal{D}_0$. First, the unstable set is computed from the vertices of the critical $1$-simplex $\sigma$ by tracing v-paths (white curves, left). The set is then collapsed into an extremum graph$\mathcal{G}_{{{0}}}$ (middle). Each critical simplex is represented by a node in the graph, with the arcs representing the v-paths. Finally, the graph is processed with a Union-Find structure to produce the persistence pair in $\mathcal{D}_0$ (right).