Sparsification of the Generalized Persistence Diagrams for Scalability Through Gradient Descent

Motivated by the computational challenges of the GPD and the flexibility of this invariant upon selecting subsets of intervals as its domain, we propose a method for automatically sparsifying GPDs computed from ${ℝ}^2$-persistence modules based on gradient descent. Namely, we consider the following scenario.

Suppose we aim at classifying instances of a given dataset based on their topological features, and we have already computed a set of corresponding persistence modules ${\{M_i:{ℝ}^2\to {\mathrm{vect}}_k\}}_{1\le i\le t}$ from the dataset, where each persistence module corresponds to an individual data point. We consider the set ${\{{\mathrm{dgm}}_{M_i}^{ℐ}\}}_{1\le i\le t}$ of GPDs of $M_i$ over a large set $ℐ$ of intervals in ${ℝ}^2$ (cf. Definition 1). We refer to these as the full GPDs.

Let $n\gg 1$ be the cardinality of $ℐ$, and let $m$ be a sparsification parameter $m\in {\mathbb{N}}^{\ast }=\{1,2,\mathrm{\dots }\}$, which is typically significantly smaller than $n$. Let $\left(\genfrac{}{}{0pt}{}{\mathrm{Int}({ℝ}^2)}{m}\right)$ denote the set of $m$-subsets of $\mathrm{Int}({ℝ}^2)$, i.e., $\{𝒥\subset \mathrm{Int}({ℝ}^2):\left|𝒥\right|=m\}$. In order to identify an $m$-subset of intervals in ${ℝ}^2$ over which the new GPDs are computed (and subsequently used to classify the persistence modules ${\{M_i\}}_{1\le i\le t}$), we proceed as follows. Firstly, we identify a loss function defined on $\left(\genfrac{}{}{0pt}{}{\mathrm{Int}({ℝ}^2)}{m}\right)$:

where $d_{\mathrm{\ell }}$ is an appropriate dissimilarity function. Secondly, we search for a minimizer of the loss function. The goal of this search is to identify a subset ${𝒥}^{\ast }$ of $m$ intervals such that the GPDs of the ${\{M_i\}}_i$ over this sparse subset ${\{{\mathrm{dgm}}_{M_i}^{{𝒥}^{\ast }}\}}_{1\le i\le t}$, the sparse GPDs, best approximate their corresponding full counterparts overall. One natural way of searching for a (local) minimizer is through gradient descent, starting from a randomly chosen $m$-subset ${𝒥}_{\text{init}}$ of $ℐ$. To achieve this, the following requirements either must be met or highly desirable:

1. (I)

   (Distance) A suitable distance or dissimilarity function $d_{\mathrm{\ell }}$, utilized in constructing the loss function above, ideally satisfying certain stability guarantees w.r.t. the interleaving distance between persistence modules [12, 25],

2. (II)

   (Vectorization) A certain representation of the loss function ${ℒ}_{d_{\mathrm{\ell }},m,{\{M_i\}}_{1\le i\le t}}$ as a map defined on some subset $𝒟$ of Euclidean space,

3. (III)

   (Convexity) Convexity of the subset $𝒟$,

4. (IV)

   (Loss regularity) Lipschitz stability and differentiability of ${ℒ}_{d_{\mathrm{\ell }},m,{\{M_i\}}_{1\le i\le t}}$, and

5. (V)

   (Feasibility) Computational feasibility of ${ℒ}_{d_{\mathrm{\ell }},m,{\{M_i\}}_{1\le i\le t}}$ for practical implementation.

Our contributions can be listed according to Items (I)-(V) mentioned above.

• $\mathrm{■}$

  A natural choice for $d_{\mathrm{\ell }}$ in Item (I) would be to use the erosion distance $d_{\mathrm{E}}$ [13, 33], a standard metric between GPDs,¹¹1 Note that $d_{\mathrm{E}}$ is also referred to as a metric between generalized rank invariants (GRIs), e.g., in [13, 21]. Since the GPD and the GRI determine one another (cf. Remark 2 (i) and (iii)), $d_{\mathrm{E}}$ can thus also be viewed as a metric between GPDs, as in [33, 36]. that is known to be stable under perturbations of input persistence modules w.r.t. the interleaving distance. However, the use of the erosion distance requires GPDs to be defined on the same set of intervals (Definition 3) that is closed under thickening, which implies the domain must be infinite. All these make it difficult to directly utilize $d_{\mathrm{E}}$. Hence, we introduce the sparse erosion distance $\widehat{d_{𝐄}}$ between GPDs relative to (possibly) different interval sets, as an adaptation of $d_{𝐄}$ (Definition 4, Proposition 6 (ii), and Corollary 7).²²2Another possibility would be to use bottleneck and Wasserstein distances, however we show in Remark 17 in the full version [11] that they fail to ensure continuity of the loss function.

• $\mathrm{■}$

  Regarding Item (V), we restrict our focus to intervals in ${ℝ}^2$ with a small number of minimal and maximal points (Remark 10). Then, if the sparse erosion distance $\widehat{d_{\mathrm{E}}}$ is computed between GPDs of the same persistence module $M$, we prove that $\widehat{d_{𝐄}}$ depends solely on the domains $ℐ\mathbf{,}𝒥$, and not on the input persistence module, i.e. $\widehat{d_{𝐄}}\mathbf{(}\mathbf{(}{\mathrm{𝐝𝐠𝐦}}_M\mathbf{,}ℐ\mathbf{)}\mathbf{,}\mathbf{(}{\mathrm{𝐝𝐠𝐦}}_M\mathbf{,}𝒥\mathbf{)}\mathbf{)}\mathbf{=}\widehat{d}\mathbf{(}ℐ\mathbf{,}𝒥\mathbf{)}$ (Proposition 6 (iii)). This fact significantly enhances the tractability of gradient descent as it allows to avoid recomputing the GPDs at every iteration. Moreover, we derive a closed-form formula for the computation of $\widehat{d}\mathbf{(}ℐ\mathbf{,}𝒥\mathbf{)}$ (Theorem 9 (i)).

• $\mathrm{■}$

  For achieving Items (II), (III), and (IV), we only consider intervals with at most two minimal points and one maximal point. This ensures the existence of natural embeddings of these intervals into Euclidean space ${ℝ}^6$, which can be obtained by stacking up the coordinates of the interval middle points $(x,y)$ together with the lengths $a,b,c,d\ge 0$ needed to define the interval lower boundaries (Figure 1). The main advantages of this vectorization method (w.r.t. the other natural ones) are that (1) it allows for a simple formulation of the loss function (Theorem 9 (ii)), and (2) its variables are all independent from each other, which makes the implementation of gradient descent easier (Remarks 10 and 11).³³3Alternatively, we can consider intervals with one minimal points and at most two maximal points. Either choice ensures that our restricted GPDs are not a weaker invariant than the rank invariant or the signed barcode [7]; see [13, Section 4]. Using this vectorization method, we also prove that $𝒟$ not only forms a convex subset of Euclidean space ${ℝ}^{\mathrm{𝟔}m}$ (Proposition 13), but also ensures the Lipschitz stability and almost everywhere differentiability of the loss function ${ℒ}_{d_{\mathbf{\ell }}\mathbf{,}m\mathbf{,}{\mathbf{\{}{M}_i\mathbf{\}}}_{\mathrm{𝟏}\mathbf{\le }i\mathbf{\le }t}}$ (Theorem 14 and Proposition 16). In fact, we prove that the loss function can actually be realized as

  (2)

  where ${𝕧}_{𝒥}$ is our $6m$-dimensional embedding of $𝒥$. Since $t$ (the size of the dataset) is a constant, the (local) minimizers of ${ℒ}_{\widehat{d_{\mathrm{E}}},m}$ coincide with those of the map $t^{-1}\cdot {ℒ}_{\widehat{d_{\mathrm{E}}},m}:{𝕧}_{𝒥}\mapsto \widehat{d}(ℐ,𝒥)$. Hence, searching for a (local) minimizer of ${ℒ}_{\widehat{d_{\mathrm{E}}},m}$ is essentially searching for the (locally) best $m$ intervals that represent the domain $ℐ$ of the full GPDs.

• $\mathrm{■}$

  Finally, regarding Item (V) again, and in order to showcase the efficiency of our proposed sparsification method, we provide numerical experiments on topology-based time series classification with supervised machine learning (Section 5 and Figure 2), in which we show that the sparse GPDs ${\{{\mathrm{dgm}}_{M_i}^{{𝒥}^{\ast }}\}}_{1\le i\le t}$ can be computed in much faster time than the full GPDs, while maintaining similar or better classification performances from random forests models. Our code is fully available at sparse GPDs.

While there exist many existing works that utilize multi-parameter persistent homology to enhance the performance of machine learning models [10, 14, 27, 29, 32, 34, 35, 37], our work takes a different approach as it focuses on methods to mitigate the computational overhead associated with multi-parameter persistence descriptors.

The works most closely related to ours are [29] and [37], which use the RI or GRI for multi-parameter persistence modules in a machine learning context. Firstly, [29] is based on an equivalent representation of the GPDs (computed from rectangle intervals⁴⁴4This types of GPDs are often called signed barcodes.) as signed measures [7], which allows to compare GPDs with optimal transport distances, as well as to deploy known vectorization techniques intended for general measures. Secondly, the approach proposed in [37] involves vectorizing the GRIs of 2-parameter persistence modules by evaluating them on intervals with specific shapes called worms. Our goal is different: we rather aim at vectorizing and sparsifying the domains of the GPDs in order to achieve good performance scores. In fact, the sparsification process that we propose can actually be used complementarily to both [29] and [37] by first sparsifying the set of rectangles or worms before applying their vectorization methods. Note that differentiability properties of both of these approaches w.r.t. the multi-parameter filtrations were recently established [32, 34]; in contrast, our work deals with the differentiability of a GPD-based loss function w.r.t. the interval domains (while keeping the multi-parameter filtrations fixed).

Section 2 reviews basic properties of the GPD and GRI. Section 3 introduces the sparse erosion distance, clarify its relation to the erosion distance given in [13], and presents its closed-form formula, which is specialized and useful in our setting. Section 4 establishes the Lipschitz stability and differentiability of our loss function. Section 5 presents our numerical experiments. Finally, Section 6 discusses future research directions.

Let $S$ be the collection of all $(1,1)$- and all $(2,1)$-intervals of ${ℝ}^2$. Consider the map from $S$ to $ℝ\times ℝ\times {ℝ}_{>0}\times {ℝ}_{\ge 0}\times {ℝ}_{\ge 0}\times {ℝ}_{>0}$ depicted in Figure 1. Clearly, this map is a surjection. It follows that the image of $S$ via the map $I\mapsto (x,y,a,b,c,d)$ depicted in Figure 1 is a convex subset of ${ℝ}^6$. Now, observe that the image of $\Im$ via the map $V$ is equal to the $n$-fold Cartesian product of the image of $S$, which is a subset of ${ℝ}^{6n}$. The fact that any Cartesian product of convex sets is convex implies our claim. $◀$

By the triangle inequality, we have $|\widehat{d}(𝒦,{𝒥}_1)-\widehat{d}(𝒦,{𝒥}_2)|\le \widehat{d}({𝒥}_1,{𝒥}_2)$ and thus it suffices to show that $\widehat{d}({𝒥}_1,{𝒥}_2)\le 2\cdot {\min }_{\pi }{\Vert {𝕧}_{\pi ({𝒥}_1)}-{𝕧}_{{𝒥}_2}\Vert }_{\mathrm{\infty }}$. Let ${\pi }_0$ be a permutation on ${𝒥}_1$ that attains the minimum. Let $ℐ:={\pi }_0({𝒥}_1)$ and $𝒥:={𝒥}_2$. By Lemma 12, ${ϵ}_{rr}$ (resp. ${ϵ}_{ss}$) $\le 2ϵ$ for all $r,s\in \{1,\mathrm{\dots },n\}$ and hence ${\min }_s{ϵ}_{rs}$ (resp. ${\min }_r{ϵ}_{rs}$) $\le 2ϵ$ for each $r$ (resp. $s$). Thus, we have ${\max }_r({\min }_s{ϵ}_{rs})\le 2ϵ\text{ and }{\max }_s({\min }_r{ϵ}_{rs})\le 2ϵ.$ By Theorem 9 (i), we are done. $◀$

The three maps ${ℝ}^2\to ℝ$ given by $(x,y)\mapsto \delta (x,y)|y-x|$ (cf. Equation (9)), $(x,y)\mapsto \max \{x,y\}$, and $(x,y)\mapsto \min \{x,y\}$ are finitely segmented piecewise linear. Therefore, the functions given in Equations (13)–(16) are all finitely segmented piecewise linear on Euclidean spaces, and thus so are ${ϵ}_{rs}$ given in Theorem 9 (ii). Hence, trivially, these functions are differentiable almost everywhere. Now, Proposition 6 (iii) implies our claim. $◀$

The datasets we consider are taken from the UCR repository [15], and correspond to classification tasks that have already been studied with persistent homology before [29, Section 4.2]. More precisely, instances in these datasets take the form of labelled time series, that we pre-process with time-delay embedding in order to turn them into point clouds. Specifically, each labelled time series $T=\{f(t_1),\mathrm{\dots },f(t_n)\}$ of length $n$ is transformed into a point cloud $X_T\subset {ℝ}^3$ of cardinality $n-2$ with

Then, we compute both the Vietoris-Rips filtration and the sublevel set filtration induced by a Gaussian kernel density estimator (with bandwidth $\sigma =0.1\cdot \mathrm{diam}(X_T)$, where $\mathrm{diam}(X_T):={\max }_{x,y\in X_T}{\Vert x-y\Vert }_2$ is the diameter of the point cloud) using the PointCloud2FilteredComplex function of the multipers library [28] (a tutorial is available in [26]). Both filtrations are then normalized so that their ranges become equal to the unit interval $[0,1]$.

In order to compute GPDs out of these 2-parameter filtrations and modules, one needs a subset of intervals. As explained in the introduction, we then minimize

where the full domain $ℐ$ (resp. the sparse domain $𝒥$) is comprised of $n=1,600$ (resp. $m=400$) $(2,1)$-intervals obtained from a grid in ${ℝ}^6$ computed by evenly sampling $10$ (resp. $5$) values for $x$ and $y$ and $2$ values for $a,b,c,d$ within their corresponding filtration ranges.⁶⁶6We also tried random uniform sampling for the initial distributions of $x,y,a,b,c,d$, but this resulted in slower convergence and lower performances. Note that while the $(2,1)$-intervals in $𝒥$ are treated as parameters to optimize, the $(2,1)$-intervals in $ℐ$ are fixed throughout the optimization process. Moreover, the formula that we provided in Theorem 9 (ii) can be readily implemented in any library that uses auto-differentiation, such as pytorch. In particular, we run stochastic gradient descent with momentum $0.9$ on $ℒ$ for $750$ epochs using learning rate $\eta =0.001$ with exponential decay of factor $0.99$ to achieve convergence, and obtain our sparse subset ${𝒥}^{\ast }$. See Figure 5 for a visualization of the loss decrease. Note how the Lipschitz stability proved in Theorem 14 translates into a smooth decrease with small oscillations.

In order to quantify the running time improvements as well as the information loss (if any) when switching from the full domain $ℐ$ to the optimized sparse one ${𝒥}^{\ast }$, we computed the full and sparse GPDs in homology dimensions $0$ and $1$ with the zigzag persistence diagram implementation in dionysus [31], and then we trained random forest classifiers on these GPDs to predict the time series labels.⁷⁷7Recall that optimizing the sparse domain can be performed without computing GPDs, so the optimization process described in the previous paragraph is done once and for all, before computing the GPDs. In order to achieve this, we first turned both the full and optimized GPDs into Euclidean vectors by first binning every GPD (seen as a point cloud in ${ℝ}^6$) with a fixed $6$-dimensional histogram, and then convolving this histogram with a $6$-dimensional Gaussian kernel in order to smooth its values, with a procedure similar to the one described in [29, Section 3.2.1]. Both the histogram bins and the kernel bandwidths were found with $3$-fold cross-validation on the training set (see the provided code for hyperparameter values). Then, random forest classifiers were trained on these smoothed histograms to predict labels; in Table 1 we report the accuracy scores of these classifiers for the initial (before optimization) sparse domain ${𝒥}_{\mathrm{init}}$, the optimized sparse domain ${𝒥}^{\ast }$, and the full domain $ℐ$. Moreover, in Table 2, we report the running time needed to compute all GPDs using either ${𝒥}_{\mathrm{init}}$, ${𝒥}^{\ast }$, or $ℐ$, as well as the improvement when switching from ${𝒥}^{\ast }$ to $ℐ$. Note that our goal is not to improve on the state-of-the-art for time series classification, but rather to assess whether optimizing the loss ${ℒ}_{\widehat{d_{\mathrm{E}}},m}$ based on our upper bound is indeed beneficial for improving topology-based models. See Figure 2 for a schematic overview of our full pipeline.

As one can see from Table 1, there is either a clear improvement or a comparable performance in accuracy scores after optimizing $𝒥$. Indeed, by minimizing ${ℒ}_{\widehat{d_{\mathrm{E}}},m}$, one forces the sparse domain $𝒥$ to be as close as possible to the full domain $ℐ$, and thus to retain as much topological information as possible. Hence, either the full GPDs are more efficient than the initial sparse GPDs, in which case the optimized GPDs perform much better than the initial sparse ones, or the full GPDs are less efficient than the initial sparse GPDs,⁸⁸8Recall that the accuracy score is only an indirect measure: While the full GPDs are always richer in topological information, we hypothesize that their scores might still be lower due to a large number of intervals on which the full GPDs vanish – a question that we leave for further study. in which case the optimized GPDs have comparable or slightly worse performances than the initial sparse ones thanks to the small sizes of their domains (except for the PC and IPD datasets, for which the optimized GPDs still perform interestingly better). In all cases, we emphasize that optimized GPDs provide the best solution: they maintain scores at levels that are either comparable or better than the best solution between the initial sparse and full GPDs, while avoiding to force users to choose interval domains (as their domains are obtained automatically with gradient descent). Note in particular that whenever optimized GPDs do not achieve the best scores, their performances remain of the same order as those of their competitors (see, e.g., DPC or SC), whereas when optimized GPDs are the best, the performances of the other two competitors are sometimes far behind (see, e.g., GPA).

As for running times, we observe a slight increase from the running times associated to the initial sparse GPDs (except for the IPD dataset), which we hypothesize is due to the use of optimized intervals that intersect more with each other, and thus induce more poset relations when computing Möbius inversions, and a strong improvement over the running times associated to the full GPDs, with ratios ranging between $5$ and $15$ times faster. Our optimized GPDs thus achieve the best of both worlds: they are strinkingly fast to compute while keeping high accuracy scores. Our code was run on a 2x Xeon SP Gold 5115 @ 2.40GHz, and is fully available at sparse GPDs.

In our experiments, we optimized the GPDs over intervals with at most two minimal points and exactly one maximal point. Allowing more complex intervals, as well as dataset-dependent terms in the loss function, could improve performance, but finding suitable embeddings of such intervals into Euclidean space remains a challenge (cf. Remark 11) and would increase the computational cost (cf. Remark 10). It would also be interesting to quantify the discriminating power of GPDs with measures that are more direct than the score of a machine learning classifier (or at least to further study the dependencies between GPD sparsification and classifier scores), and to investigate on heuristics (based on drops in the loss values) for deciding whether sparse GPDs are sufficiently good so that optimization can be stopped.

While we only focused on sparsifying GPDs in this work, it would be interesting to measure the extent to which our interval domain optimization adapts to other descriptors based on the GRI from the TDA literature; of particular interest are the GPDs/signed barcodes coming from rank exact decompositions, which have recently proved to be stable [8] (see also Section 8.2 in [7] which discusses the influence of the choice of the interval domains on the resulting invariants), as well as GRIL, which focuses on specific intervals called worms [37].