Persistent homology of unweighted complex networks via discrete Morse theory

TL;DR

This paper introduces a discrete Morse theory-based method for analyzing the topological features of unweighted networks using persistent homology, improving computational efficiency and ability to distinguish network structures.

Contribution

The paper presents a novel approach leveraging discrete Morse theory to efficiently compute persistent homology of unweighted networks, capturing higher-order topological features.

Findings

Achieves near-minimal critical simplices in various networks

Significantly improves computational efficiency of persistent homology

Can distinguish different network models based on topology

Abstract

Topological data analysis can reveal higher-order structure beyond pairwise connections between vertices in complex networks. We present a new method based on discrete Morse theory to study topological properties of unweighted and undirected networks using persistent homology. Leveraging on the features of discrete Morse theory, our method not only captures the topology of the clique complex of such graphs via the concept of critical simplices, but also achieves close to the theoretical minimum number of critical simplices in several analyzed model and real networks. This leads to a reduced filtration scheme based on the subsequence of the corresponding critical weights, thereby leading to a significant increase in computational efficiency. We have employed our filtration scheme to explore the persistent homology of several model and real-world networks. In particular, we show that our method can detect differences in the higher-order structure of networks, and the corresponding persistence diagrams can be used to distinguish between different model networks. In summary, our method based on discrete Morse theory further increases the applicability of persistent homology to investigate the global topology of complex networks.

Tables

Table 1. Table 2: The table lists the value of optimality indicator μ 𝜇 \mu for various model and real networks analyzed here. Equation 14 gives the definition of μ 𝜇 \mu which is an indicator for the proximity of our algorithm to the optimal case. The value of μ 𝜇 \mu ranges from 0 0 to 1 1 1 with μ = 1 𝜇 1 \mu=1 indicating that our algorithm achieves exactly the theoretical minimum number of critical simplices while μ = 0 𝜇 0 \mu=0 indicating that all simplices are critical. For model networks, the value reported in this table is the average of μ 𝜇 \mu across 10 samples of each model network for a chosen set of parameter values along with the corresponding standard deviations.

Network	Optimalitiy indicator $\mu$
ER model with $n=1000,p=0.004$	0.924 $\pm$ 0.004
ER model with $n=1000,p=0.006$	0.947 $\pm$ 0.004
ER model with $n=1000,p=0.008$	0.959 $\pm$ 0.006
WS model with $n=1000,k=4,p=0.5$	0.890 $\pm$ 0.003
WS model with $n=1000,k=6,p=0.5$	0.917 $\pm$ 0.007
WS model with $n=1000,k=8,p=0.5$	0.906 $\pm$ 0.003
BA model with $n=1000,m=2$	0.989 $\pm$ 0.003
BA model with $n=1000,m=3$	0.985 $\pm$ 0.004
BA model with $n=1000,m=4$	0.964 $\pm$ 0.006
Spherical model with $n=1000,\gamma =\infty ,T=0,k=4$	0.925 $\pm$ 0.007
Spherical model with $n=1000,\gamma =\infty ,T=0,k=6$	0.901 $\pm$ 0.003
Spherical model with $n=1000,\gamma =\infty ,T=0,k=8$	0.887 $\pm$ 0.003
Hyperbolic model with $n=1000,\gamma =2,T=0,k=4$	0.939 $\pm$ 0.013
Hyperbolic model with $n=1000,\gamma =2,T=0,k=6$	0.927 $\pm$ 0.007
Hyperbolic model with $n=1000,\gamma =2,T=0,k=8$	0.921 $\pm$ 0.005
US Power Grid	0.893937
Email communication	0.871847
Route views	0.952140
Yeast protein interactions	0.942157
Hamsterster friendship	0.793236
Euro road	0.840678
Human protein interaction	0.957924

Equations

$$\varnothing \subseteq K^0\subseteq K^1\subseteq ...\subseteq K^q=K$$

$$H_p^{i,j}=Z_p^i/(B_p^{i+j}\cap Z_p^i).$$

$${\beta }_p^{i,j}=\text{dim}(H_p^{i,j}).$$

$$U_{\alpha }^f=\{{\beta }^{p+1}|{\alpha }^p<{\beta }^{p+1}\text{ and }f(\beta )\le f(\alpha )\}$$

$$V_{\alpha }^f=\{{\gamma }^{p-1}|{\gamma }^{p-1}<{\alpha }^p\text{ and }f(\alpha )\le f(\gamma )\}$$

$$|U_{\alpha }^f|\le 1\text{ and }|V_{\alpha }^f|\le 1.$$

$$|U_{\alpha }^f|=0\text{ and }|V_{\alpha }^f|=0.$$

$$K(r)=\bigcup _{f(\beta )\le r}\bigcup _{\alpha \le \beta }\alpha$$

$$\varnothing \subseteq K(w_0)\subseteq K(w_1)\subseteq \cdots \subseteq K(w_{n-1})\subseteq K(w_n)=K.$$

$$\varnothing \subseteq K(w_{c_0})\subseteq K(w_{c_1})\subseteq \cdots \subseteq K(w_{c_{m-1}})\subseteq K(w_{c_m})\subseteq K.$$

$$W_{\infty }(X,Y)=\underset{\eta :X\to Y}{\inf }{\text{sup}}_{x\in X}||x-\eta (x)|{|}_{\infty }.$$

$$W_{q}(X,Y)=\bigg{[}\inf_{\eta:X\rightarrow Y}\sum_{x\in X}||x-\eta(x)||_{\infty}^{q}\bigg{]}^{\frac{1}{q}}.$$

$$g(v)=de{g}_{max}-\text{degree}(v)+ϵ$$

$$\mu =\frac{\sum _{p=0}^d{n}_p-\sum _{p=0}^d{m}_p}{\sum _{p=0}^d{n}_p-\sum _{p=0}^d{\beta }_p}$$

$${\partial }_p(\alpha )=\sum _{i=0}^p(-1{)}^i[v_0,v_1,\cdots ,\widehat{v_i},\cdots ,v_p],$$

$${\partial }_p(c)=\sum _i{n}_i{\partial }_p({\alpha }_i)$$

$$Z_p=\text{ker}({\partial }_p)=\{c\in C_p|{\partial }_p(c)=0\}.$$

$$B_p=\text{img}({\partial }_{p+1})=\{c\in C_p|\exists b\in C_{p+1},{\partial }_{p+1}(b)=c\}.$$

$$H_p=Z_p/B_p,$$

$$\chi (K)={\beta }_0-{\beta }_1+{\beta }_2-\cdots +(-1{)}^d{\beta }_d.$$

$$U_{\alpha }=\{{\beta }^{p+1}|{\alpha }^p<{\beta }^{p+1}\text{ and }f(\beta )\le f(\alpha )\}$$

$$V_{\alpha }=\{{\gamma }^{p-1}|{\gamma }^{p-1}<{\alpha }^p\text{ and }f(\alpha )\le f(\gamma )\}$$

$$f(\beta )=\{\begin{array}{ll}(f({\gamma }_0)+f({\gamma }_1))/2& \text{ if }\text{Flag}[{\gamma }_0]=0\text{ and }f({\gamma }_0)>f({\gamma }_1)\text{(Case A)}\\ f({\gamma }_0)+ϵ& \text{ otherwise}\text{(Case B)}\end{array}$$

$$\#\{{\beta }^{p+1}|{\alpha }^p<{\beta }^{p+1}\text{ and }f(\beta )<f(\alpha )\}\le 1.$$

$$\#\{{\beta }^{p+1}|{\alpha }^p<{\beta }^{p+1}\text{ and }f(\beta )<f(\alpha )\}=\#\{{\beta }^{p+1}|{\alpha }^p<{\beta }^{p+1}\text{ and }f(\beta )\le f(\alpha )\}$$

$$f(\alpha )=\{\begin{array}{ll}(f({\gamma }_0)+f({\gamma }_1))/2& \text{ if }\text{Flag}[{\gamma }_0]=0\text{ and }f({\gamma }_0)>f({\gamma }_1)\text{(Case A)}\\ f({\gamma }_0)+ϵ& \text{ otherwise}\text{(Case B)}\end{array}$$

$$\text{C0}:f(\alpha )<f(\beta )\iff g(\alpha )<g(\beta )\text{ for a p-simplex }\alpha \text{ and a (p+1)-simplex }\beta \text{ such that }{\alpha }^p<{\beta }^{p+1}$$

$$\text{C1}:f(\alpha )\le f(\beta )\iff g(\alpha )\le g(\beta )\text{ for any pair of simplices }\alpha ,\beta \in K.$$

$$\text{C2:}\text{ in addition to C1, we have }f(c_i)=g(c_i)\text{ for any }{c}_i\in C$$

$$W_{\infty }(X,Y)=\underset{\eta :X\to Y}{\inf }{\text{sup}}_{x\in X}||x-\eta (x)|{|}_{\infty }.$$

$${\mu }_k(w_i,w_j):=\underset{ϵ\to 0^{+}}{\lim }({\beta }_{w_i+ϵ}^{w_j-ϵ}-{\beta }_{w_i+ϵ}^{w_j+ϵ}+{\beta }_{w_i-ϵ}^{w_j+ϵ}-{\beta }_{w_i-ϵ}^{w_j-ϵ})$$

Full text

Introduction

In recent years, the field of topological data analysis (TDA) has rapidly grown to provide a set of powerful tools to analyze various important features of data Carlsson2009 . In this context, persistent homology has played a key role in bringing TDA to the fore of modern data analysis. It not only gives a way to visualize data efficiently, but also to extract relevant information from both structured and unstructured datasets. This crucial aspect has been used effectively in various applications from astrophysics (e.g., determination of inter-galactic filament structures Pranav2016 ) to imaging analysis (e.g., feature detection in 3D gray-scale images Gunther2011 ) to biology (e.g., detection of breast cancer type with high survival rates Nicolau2011 ). Informally, the essence of the theory is its power to extract the shape of data, as well as infer higher-order correlations between various parts of the data at hand which are missed by other classical techniques Carlsson2009 . The basic mathematical theory used in this subject is that of algebraic topology, and in particular the study of homology, developed by the French mathematician Henri Poincaré at the turn of the 20th century. The origins of persistent homology lie in the ideas of Morse theory Morse1934 , which gives a powerful tool to detect the topological features of a given space through the computation of homology using real-valued functions on the space. We refer the reader to the survey article Edelsbrunner2008 for further details.

On the other hand, the discretized version of Morse theory developed by Robin Forman Forman1995 ; Forman1998 ; Forman2002 , gives a way to characterize the homology group of a simplicial complex in terms of a real-valued function with certain properties, known as a discrete Morse function. Examples of such simplicial complexes associated with discrete spaces are the Vietoris-Rips complex corresponding to a discrete metric space, or the clique complex of a graph. Forman Forman1998 ; Forman2002 showed that given such a function, the so-called critical simplices completely determine the Euler characteristic of the space, which is a fundamental topological invariant.

The study of complex networks in the last few decades has also significantly raised our ability to understand various kinds of interactions arising in both natural and artificial realms Watts1998 ; Barabasi1999 ; Albert2002 ; Newman2010 . Understanding how different parts of networks behave and influence each other is therefore an important problem Watts1998 ; Barabasi1999 ; Albert2002 ; Newman2010 . However, for large networks, detecting higher-order structures remains a difficult task Bianconi2015 . Moreover, recent studies Kartun-Giles2019 ; Iacopini2019 ; Ritchie2017 indicate that these higher-order correlations are not captured by usual network measures such as clustering coefficients. While a graph representation captures binary relationships among vertices of a network, simplicial complexes also reflect higher-order relationships in a complex network De2007 ; Horak2009 ; Petri2013 ; Petri2014 ; Wu2015 ; Sizemore2016 ; Courtney2017 ; Courtney2018 ; Kartun-Giles2019 ; Iacopini2019 . In this context, persistent homology has been employed to explore the topological properties of complex networks De2007 ; Horak2009 ; Lee2012 ; Petri2013 ; Petri2014 ; Sizemore2016 . In this work, we present a systematic method to study the persistent homology of unweighted and undirected graphs or networks.

Previous work has investigated the persistent homology of weighted and undirected networks by creating a filtration of the clique complexes corresponding to threshold graphs obtained via decreasing sequence of edge weights Petri2013 ; Sizemore2016 . However, the lack of edge weights in unweighted networks does not permit a filtration based on threshold graphs Petri2013 ; Sizemore2016 . Thus for unweighted networks, Horak et al. Horak2009 propose a filtration scheme based on the dimension of the simplices in the clique complex corresponding to the unweighted network. Horak et al. Horak2009 do not assign weights to vertices, edges or higher-dimensional simplices in the clique complex corresponding to an unweighted graph. An unexplored filtration scheme involves transforming an unweighted network into a weighted network by assigning edge weights based on some network property, such as edge betweenness centrality Freeman1977 ; Girvan2002 or discrete edge curvature Sreejith2016 ; Samal2018 , and then employing the filtration scheme based on threshold graphs Petri2013 ; Sizemore2016 . As an alternative, we here use discrete Morse theory Forman1995 ; Forman1998 ; Forman2002 to create a filtration scheme for unweighted networks by assigning weights to vertices, edges, triangles and higher-dimensional simplices in the clique complex of the graph. In our method, the weight of a simplex is chosen such that it reflects the degree of the vertices which constitute the simplex while simultaneously satisfying the conditions for the weighing function to be a discrete Morse function. Moreover, as explained in the Results section, an equally important intuition behind the choice of these weights is based on the goal of reducing the number of so-called critical simplices.

In the context of TDA, classical Morse theory, which involves smooth functions defined on topological spaces that admit a smooth structure has been used to compute persistent homology, e.g. in statistical topology Bubenik2009 , astrophysics Pranav2016 . Since the clique complex of a weighted or unweighted graph does not permit a smooth structure in general, applying classical Morse theory is not possible. However, discrete Morse theory Forman1995 ; Forman1998 ; Forman2002 provides an efficient way of computing persistent homology. A discrete Morse function not only captures higher-order topological information of the underlying space, a “preprocessing” with respect to a suitable discrete Morse function leads to significant simplification of the topological structure. This makes computation of persistent homology groups or homology groups of filtered complexes much more efficient, see e.g. Mischaikow2013 . This is especially useful for large datasets where computationally efficient methods are key to compute their persistent homology.

We have applied discrete Morse theory to compute persistent homology of unweighted simple graphs. This is done by using the values given by the discrete Morse function to pass from an unweighted graph to a weighted simplicial complex (Figure 1). This transformation automatically produces a filtration that is needed for the computation of persistent homology, through the so-called level subcomplexes associated with weights of critical simplices (See Theory section and Figure 2). Moreover, this filtration is consistent with the topology of the underlying space and reveals finer topological features than the dimensional filtration scheme used in Horak et al.Horak2009 The combination of these techniques have been used Gunther2011 ; Delgado2014 with applications for image processing. However, to the best of our knowledge, this method has not been used for studying persistent homology in unweighted complex networks to date. Discrete Morse theory gives a theoretical lower bound on the number of critical simplices which can be attained by an optimal choice of the function on a simplicial complex. Interestingly, our method achieves close to the theoretical minimum number of critical simplices in several model and real networks analyzed here (See Results section). Furthermore, our algorithm for computing the discrete Morse function is easy to implement for complex networks.

Our results underline the potence of persistent homology to detect inherent topological features of networks which are not directly captured by homology alone. For instance, the $p$-Betti numbers of the clique complexes corresponding to small-world Watts1998 and scale-free Barabasi1999 networks with similar size and average degree, respectively, are of comparable magnitude and thus, homology reveals no deep insight into the differences between the topological features of these two model networks. On the other hand, our observations on the persistent homology of these two networks indicate a clear demarkation with respect to the evolution of topological features in the clique complexes corresponding to these model networks during the filtration process. This dissimilarity in the evolution of topological characteristics that resonates across dimensions and the average degree of the underlying network, indicates an inherent disparity in the persistent homology of small-world and scale-free graphs. In addition to unravelling higher-order relationships in networks, this ability to capture inherent topological differences between two dissimilar networks thus motivates the application of our methods to study the persistent homology of real-world networks.

The remainder of the paper is organized as follows. We begin with a Theory section which gives a brief overview of concepts in persistent homology and discrete Morse theory. We then proceed to describe the model networks and real-world networks that have been studied in this work in the Network datasets section. In the subsequent section on Results and Discussion, we present our algorithm to construct a discrete Morse function on a simplical complex associated with a network. In the same section, we present our results for model networks and real-world networks. The final section on Conclusions gives a summary and outlook of our findings. In Supplementary Information (SI) Appendix, we give a brief review of the mathematical theory of homology groups of a simplicial complex. In SI Appendix, we also provide a rigorous proof of concept for our algorithm to construct a discrete Morse function. We then follow up with two algorithms both of which illustrate key procedures that are essential to construct the filtration of a simplical complex associated with the investigated networks. In SI Appendix, time complexity for our algorithms as well as some theoretical results on stability of persistent homology and persistence diagrams are also given.

Theory

Graphs and Simplicial Complexes

Consider a finite simple graph $G(\mathcal{V},\mathcal{E})$ having vertex set $\mathcal{V}=\{v_{0},v_{1},\cdots,v_{n}\}$ and the edge set $\mathcal{E}$. Note that a simple graph does not contain self-loops or multi-edges Bollobas1998 . Such a simple graph $G$ can be viewed as a clique complex $K$ Zomorodian2005 . A clique simplicial complex $K$ is a collection of simplices where a $p$-dimensional simplex (or $p$-simplex) in $K$ is a set of $p+1$ vertices that form a complete subgraph. In other words, vertices correspond to [math]-simplices, edges to $1$-simplices, and triangles to $2$-simplices in the clique complex of a graph. Note that the dimension $p$ of simplices contained in $K$ is restricted to the range 0 to $(|\mathcal{V}|-1)$ in the graph $G$. The dimension $d$ of the clique complex $K$ is given by the maximum dimension of its constituent simplices. A face $\gamma$ of a $p$-simplex $\alpha$ is a subset of $\alpha$ with cardinality less than $p+1$. Note that by definition, a face $\gamma$ of a $p$-simplex $\alpha$ is a $l$-simplex where $0\leq l<p$ and this relationship is denoted as $\gamma^{l}<\alpha^{p}$. Formally, the clique complex $K$ corresponding to the simple graph $G$ satisfies the following condition which defines an abstract simplicial complex, namely, $K$ is a collection of non-empty finite sets or simplices such that if $\alpha$ is an element (simplex) of $K$ then so is every non-empty subset of $\alpha$. For additional details, the interested reader is referred to standard text in algebraic topology Munkres2018 .

Figure 1 displays an example of the correspondence between a simple graph and its clique complex. The ordering of the vertex set $\{v_{0},v_{1},v_{2},\cdots,v_{p}\}$ of a $p$-simplex $\alpha$ determines its orientation. Moreover, two orderings of the vertex set of $\alpha$ are considered to be equivalent if and only if they differ by an even permutation. If the dimension of a $p$-simplex is greater than 1, then all possible orderings of its vertex set fall under two equivalence classes, with each class being assigned an orientation Munkres2018 . An exception is the [math]-simplex with one vertex which has exactly one equivalence class and orientation. An oriented $p$-simplex $\alpha$ specifies the orientation of its $p+1$ vertices and is represented by $[v_{0},v_{1},v_{2},\cdots,v_{p}]$ Munkres2018 . In figure 1, the oriented $2$-simplices $[v_{2},v_{3},v_{4}]$ and $[v_{2},v_{4},v_{3}]$ have opposite orientations, i.e., $[v_{2},v_{3},v_{4}]=-[v_{2},v_{4},v_{3}]$.

Persistent homology of a simplicial complex

In SI Appendix, we give a brief review of the mathematical theory of homology groups of a simplicial complex. In particular, we define $p$-chain group, $p$-boundary operator, $p$-boundary, $p$-cycle, $p$-hole, $p$-homology group and $p$-Betti number.

A subset $K^{i}$ of a simplicial complex $K$ is called a subcomplex of $K$ if $K^{i}$ by itself is an abstract simplicial complex. Then a filtration of a simplicial complex $K$ is defined as a nested sequence of subcomplexes $K^{i}$ of $K$ where:

[TABLE]

Note that each subcomplex $K^{i}$ has an associated index $i$ in the filtration. Moreover, each subcomplex $K^{i}$ in the filtration has corresponding $p$-chain complexes $C_{p}^{i}$, $p$-boundary operators $\partial_{p}^{i}$, $p$-boundaries $B_{p}^{i}$ and $p$-cycles $Z_{p}^{i}$.

The $j$-persistent $p$-homology group of $K^{i}$ denoted as $H_{p}^{i,j}$ is defined as:

[TABLE]

In the above equation, $B_{p}^{i+j}$ is the subgroup of $C_{p}^{i+j}$ which constitutes the $p$-boundaries of the subcomplex $K^{i+j}$. The $j$-persistent $p$-Betti number of $K^{i}$ denoted as $\beta_{p}^{i,j}$ is defined as:

[TABLE]

An intuitive explanation of the above definitions of the $j$-persistent $p$-homology group and the corresponding Betti number is as follows. A $p$-hole of the subcomplex $K^{i}$ can potentially become the boundary of a $(p+1)$-chain of a later subcomplex $K^{i+j}$ with $j>0$, and thus, no longer constitute a $p$-hole of $K^{i+j}$. The $j$-persistent $p$-Betti number of $K^{i}$ represents the number of $p$-holes at the filtration index $i$ that persist at the filtration index $j+i$. Therefore, each $p$-hole that appears across the filtration has a unique index that corresponds to its birth and death, and the persistence of such a $p$-hole can thus be characterized by its corresponding birth and death indices. Studying persistent homology allows us to quantify the longevity of such $p$-holes during filtration, and thus, measures the importance of these topological features which appear and disappear across the filtration.

Discrete Morse Theory

Recalling from the preceeding section, to study the persistent homology corresponding to the clique complex $K$ of a simple graph $G$, the primary requirement is a filtration of $K$. We here present a systematic method to study the persistent homology of unweighted and undirected networks by utilizing a refined filtration of the clique complex $K$ based on discrete Morse theory Forman1995 ; Forman1998 ; Forman2002 . It is important to note that the order in which the simplices of the clique complex $K$ are added during the filtration affects the evolution of topological features which are observed by studying the persistent homology. Our proposed scheme which is based on discrete Morse theory developed by Forman Forman1995 ; Forman1998 ; Forman2002 tackles this by assigning weights to [math]-simplices (vertices), $1$-simplices (edges), $2$-simplices (triangles) and higher-dimensional simplices appearing in the clique complex corresponding to an unweighted and undirected network. Assigning weights to higher-dimensional simplices captures important higher-order correlations in addition to edges or $1$-simplices. Moreover, we leverage the following important features of the framework of discrete Morse theory in our new scheme. Firstly, the framework enables assignment of weights to $p$-simplices which are concordant with weights of $(p-1)$-simplices. Secondly, it captures the topology of a simplicial complex via the concept of critical simplices described below. Most importantly, the framework provides a natural way to create a filtration scheme to study persistent homology based upon the weights of the aforementioned critical simplices as will be described below.

We next provide the fundamental definitions in discrete Morse theory Forman1998 ; Forman2002 . We remark that a $p$-dimensional simplex $\alpha$ in a simplicial complex $K$ is denoted by $\alpha^{p}\in K$. Also, if a $p$-simplex $\alpha^{p}$ in $K$ is a face of a $(p+1)$-simplex $\beta^{p+1}$ in $K$ then this is represented as $\alpha^{p}<\beta^{p+1}$ in the sequel. Given a function $f:K\rightarrow\mathbb{R}$, for each simplex $\alpha^{p}\in K$, two sets $U_{\alpha}^{f}$ and $V_{\alpha}^{f}$ are defined as follows:

[TABLE]

Simply stated, the set $U_{\alpha}^{f}$ contains any $(p+1)$-simplex $\beta^{p+1}$ of which $\alpha^{p}$ is a face and the function value on $\beta$ is less than or equal to the function value on $\alpha$. The set $V_{\alpha}^{f}$ contains any $(p-1)$-simplex $\gamma^{p-1}$ which is a face of $\alpha^{p}$ and the function value on $\alpha$ is less than or equal to the function value on $\gamma$. A function $f:K\rightarrow\mathbb{R}$ is a discrete Morse function Forman1998 ; Forman2002 if and only if for each simplex $\alpha^{p}\in K$:

[TABLE]

Given a discrete Morse function $f$ on the simplicial complex $K$, a simplex $\alpha^{p}\in K$ is critical Forman1998 ; Forman2002 if and only if:

[TABLE]

Simply stated, a $p$-simplex $\alpha^{p}\in K$ is critical if the following conditions are simultaneously satisfied. The first condition being that if $\beta^{p+1}$ is any $(p+1)$-simplex in $K$ of which $\alpha^{p}$ is a face, then $f(\alpha)<f(\beta)$. The second condition being that if $\gamma^{p-1}$ is any $(p-1)$-simplex in $K$ which is a face of $\alpha^{p}$ then $f(\alpha)>f(\gamma)$. The concept of critical simplices in discrete Morse theory is in spirit a discrete analogue to the concept of critical points in classical Morse theory wherein the critical points corresponding to a smooth real valued function on $X\subseteq\mathbb{R}^{d}$ are the points where the gradient of the function vanishes.

We remark that once a discrete Morse function $f$ on a simplicial complex $K$ is fixed, the sets $U_{\alpha}^{f}$ and $V_{\alpha}^{f}$ are denoted by $U_{\alpha}$ and $V_{\alpha}$, respectively, to simplify the notation. A simple example for a discrete Morse function on a simplicial complex $K$ is the dimension function used in Horak et al Horak2009 . The value of the dimension function on a given simplex $\alpha$ is the dimension of the simplex $\alpha$. By a direct consequence of the definitions presented above in this section, for every simplex $\alpha\in K$, the sets $U_{\alpha}$ and $V_{\alpha}$ corresponding to the dimension function are empty. Thus, the dimension function is indeed a discrete Morse function and every simplex in $K$ is critical. In the results section, we present our new scheme and algorithm 1 to assign a discrete Morse function $f$ to a clique complex $K$ of an unweighted graph $G$.

We next describe the filtration of the clique simplicial complex $K$ based on the discrete Morse function $f$. Given a discrete Morse function $f$ on a simplicial complex $K$ and a real number $r$, a level subcomplex $K(r)$ is defined Forman1998 ; Forman2002 as follows:

[TABLE]

Simply stated, $K(r)$ contains all simplices $\beta$ in $K$ with the value of the discrete Morse function or assigned weight $f(\beta)\leq r$ along with any face $\alpha$ of $\beta$. Note that a face $\alpha$ of $\beta$ is included in $K(r)$ even if the discrete Morse function or assigned weight to a face $\alpha$ is greater than $r$.

Let $\{f(\sigma)\}_{\sigma\in K}$ denote the entire set of values assigned to simplices in $K$ using the discrete morse function $f$. Then, let $\{w_{k}\}_{k=0,\cdots,n}$ denote the finite increasing sequence of the unique values in the set $\{f(\sigma)\}_{\sigma\in K}$ associated with the finite simplicial complex considered here. We now have a sequence of inclusions of level subcomplexes corresponding to this increasing sequence $\{w_{k}\}$ as follows:

[TABLE]

This nested sequence gives a filtration of the simplicial complex $K$ which enables the study of persistent homology in the context of unweighted networks.

According to Lemma 2.6 by Forman Forman2002 , if there are no critical simplices $\alpha$ with $f(\alpha)\in(a,b],$ then $K(b)$ is homotopy equivalent to $K(a)$.

The implications of this Lemma are as follows. Let $\{f(\sigma_{c})\}$ denote the set of values assigned to critical simplices $\sigma_{c}$ in $K$ by the discrete Morse function $f$ and $\{w_{c_{k}}\}_{k=0,\cdots,m}$ denote the increasing sequence of the unique values in $\{f(\sigma_{c})\}$. We refer to the function values $\{f(\sigma_{c})\}$ assigned to the critical simplices $\sigma_{c}$ in $K$ as critical weights. Note that the set $\{f(\sigma_{c})\}$ defined for critical simplices is a subset of $\{f(\sigma)\}$ defined for all simplices in $K$ and the increasing sequence $\{w_{c_{k}}\}_{k=0,\cdots,m}$ is a subsequence of $\{w_{k}\}_{k=0,\cdots,n}$ with $m\leq n$. The above definition implies that there are no critical simplices $\alpha$ with $f(\alpha)\in(w_{c_{i}},w_{c_{i+1}})$. As homology is invariant under homotopy equivalence, Forman’s Lemma 2.6 gives us that for any $x$ and $y$ belonging to the real number interval $(w_{c_{i}},w_{c_{i+1}})$, the homology groups of $K(x)$ and $K(y)$ are isomorphic. Thus, in order to observe the changes in homology as the filtration proceeds, it suffices to study the persistent homology of a filtration which corresponds to the subsequence $\{w_{c_{k}}\}_{k=0,\cdots,m}$ of $\{w_{k}\}_{k=0,\cdots,n}$, where $m\leq n$, and this results in a potential decrease in the required number of filtration steps. The new filtration sequence can be represented as:

[TABLE]

Note that each simplex $\alpha$ in the clique complex $K$ is first introduced as part of certain level subcomplex $K(w_{c_{i}})$ in the above nested filtration sequence. Therefore, each simplex $\alpha$ in $K$ can be associated with a unique weight $w_{c_{i}}$ referred to as the filtration weight of $\alpha$. In SI Appendix, we present algorithm 2 and algorithm 3 which depict the procedure to compute the filtration weights of simplices in the clique complex $K$ of a graph $G$. In SI Appendix, we also give a sufficient condition for two discrete Morse functions to induce the same filtration and thus the persistent homology groups.

Using an example network in figure 2, we also show that the persistent homology observed using the filtration based on the entire sequence of weights satisfying discrete Morse function is equivalent to that observed using the filtration based on the subsequence of critical weights.

Let $m_{p}$ represent the number of critical $p$-simplices in a simplicial complex $K$ and let $\beta_{p}$ denote the $p$-Betti number of $K$. Then Theorem 2.11 by Forman Forman2002 can be stated as follows.

• (i)

For each $p=0,1,2,\cdots,d$ (where $d$ is the dimension of $K$), $m_{p}\geq\beta_{p}$.

• (ii)

$m_{0}-m_{1}+m_{2}-\cdots+(-1)^{d}m_{d}=\beta_{0}-\beta_{1}+\beta_{2}-\cdots+(-1)^{d}\beta_{d}$.

In other words, the above theorem gives a lower bound of the number of $p$-critical simplices $m_{p}$ for each dimension $p$ as the $p$-Betti number $\beta_{p}$ of $K$. In results section, we present our algorithm 1 to assign weights satisfying discrete Morse function to simplices in the clique complex $K$ of a graph $G$. Our choice of the function in algorithm 1 to assign weights to simplices in the clique complex $K$ tries to minimize the number of critical simplices (which has a lower bound given by Forman’s Theorem 2.11 Forman2002 ), and thus, reduces the number of filtration steps required to compute the persistent homology without loss of information. In the results section, we will show that our algorithm achieves near-optimal number of critical weights in clique complexes corresponding to many model and real networks analyzed here.

Comparing Persistence diagrams

Given a discrete Morse function $f$ and its associated filtration $\{K(w_{c_{k}})\}$ of the clique complex $K$ of a graph $G$ (equation 10), each $p$-hole has a critical weight $w_{c{{}_{birth}}}$ which corresponds to its birth index and $w_{c{{}_{death}}}$ which corresponds to its death index, with $w_{c{{}_{birth}}}<w_{c{{}_{death}}}$.

Persistence diagram D$(f)$ for a $d$-dimensional simplicial complex $K$ is the collection of points in $\mathbb{R}^{2}$ whose first and second coordinates, $x$ and $y$, respectively, correspond to the birth weight and death weight of a $p$-hole where $0\leq p\leq d$ Cohen-Steiner2007 . Since two different holes can have the same birth and death weights, each point in the persistence diagram has a corresponding multiplicity, we refer the reader to the SI Appendix for more details. Thus, the persistence diagram is a multiset of points in $\mathbb{R}^{2}$. The persistence of a $p$-hole which has birth and death weights, $w_{c{{}_{birth}}}$ and $w_{c{{}_{death}}}$, respectively, is defined as $w_{c{{}_{birth}}}-w_{c{{}_{death}}}$. Thus, the persistence diagram for a clique complex $K$ corresponding to a graph $G$ is a compact representation of the persistent homology of a network.

Given two persistence diagrams $X$ and $Y$ (which may correspond to two different networks), the $\infty$-Wasserstein distance between $X$ and $Y$, also known as the bottleneck distance Cohen-Steiner2007 , is defined as follows:

[TABLE]

Similarly, given two persistence diagrams $X$ and $Y$, the $q$-Wasserstein distance Kerber2017 between $X$ and $Y$ is defined as follows:

[TABLE]

In the above equations, $\eta$ ranges over all bijective maps from $X$ to $Y$, and given $(a,b)\in\mathbb{R}^{2}$, $||(a,b)||_{\infty}=\text{max}\{|a|,|b|\}$ is the $L_{\infty}$ norm. In this work, we use Dionysus 2 package (http://www.mrzv.org/software/dionysus2/) to compute the Wasserstein distance between two persistence diagrams corresponding to two different model networks (See results section). Note that it is not generally true that two persistence diagrams $X$ and $Y$ have the same number of off-diagonal points, i.e., features with non-zero persistence, and we refer the readers to Kerber et al. Kerber2017 for details on circumventing this issue and further information regarding how the computation of the Wasserstein distance is reduced to a bipartite graph matching problem in Dionysus 2 package. We remark that the bottleneck distance between two persistence diagrams which are subsets of the unit square is in the range 0 to 1.

Stability of persistence diagrams with respect to small changes in the discrete Morse function is a key property of persistent homology. The first such stability theorem was given by Cohen-Steiner et al. Cohen-Steiner2007 , and later refined by Chazal et al. Chazal2008 . In the SI Appendix, using results from Chazal et al.Chazal2008 , we give a stability result for persistence diagrams of discrete Morse functions with respect to the bottleneck distance.

Network datasets

Model networks. We have investigated the following models of unweighted and undirected networks, namely, the Erdös-Rényi (ER) Erdos1961 , the Watts-Strogatz (WS) Watts1998 , the Barabási-Albert (BA) Barabasi1999 and the Hyperbolic Graph Generator (HGG) Krioukov2010 . The ER model Erdos1961 is characterized by the property that the probability $p$ of the existence of each possible edge between any two vertices among the $n$ vertices in the graph $G$ is constant. The existence of edges in the ER model are independent of each other, and thus, the model produces random graphs $G(n,p)$ with average vertex degree $p(n-1)$. The WS model Watts1998 produces small-world graphs as follows. The WS model starts with an initial regular graph with $n$ vertices where each vertex is connected to its $k$ nearest neighbours. Next, the endpoint of each edge in the initial regular graph of the WS model is randomly chosen for rewiring based on a fixed rewiring probability $p$ and is rewired to another vertex in the graph which is chosen with uniform probability. The BA model Barabasi1999 produces scale-free graphs which are characterized by a degree distribution that follows a power law decay. The BA model utilizes a preferential attachment scheme to produce scale-free graphs. The BA model generates an initial graph of $m_{0}$ vertices, and then, at each successive iteration a new vertex is added with edges to $m$ already existing vertices which are chosen with probability proportional to their degree at that particular iteration. The iterations in the BA model cease when the graph has attained the requisite number $n$ of vertices. The HGG model Krioukov2010 ; Aldecoa2015 produces a random graph of $n$ vertices by initially fixing $n$ vertices to $n$ points on a hyperbolic disk. In the HGG model, the probability of existence of an edge between two vertices is proportional to the hyperbolic distance between the two points on the hyperbolic disk that correspond to these two vertices. By tuning the input parameter $\gamma$, the HGG model can produce either a hyperbolic or a spherical random graph Krioukov2010 ; Aldecoa2015 . Specifically, the HGG model produces hyperbolic random graphs for $\gamma\in[2,\infty)$ whereas spherical random graphs for $\gamma=\infty$.

Real networks. We have also studied seven real-world networks which are represented as unweighted and undirected graphs. We have considered two biological networks, namely, the Yeast protein interaction network Jeong2001 with 1870 vertices and 2277 edges, and the Human protein interaction network Rual2005 with 3133 vertices and 6726 edges. In both biological networks, each vertex represents a protein and an edge represents an interaction between the two proteins. We have considered two infrastructure networks, namely, the US Power Grid network Leskovec2007 and the Euro road network Subelj2011 . In the US Power Grid network, the 4941 vertices represent the generators, transformers and substations in the Western states of USA and the 6594 edges represent power links between them. The 1174 vertices of the Euro road network correspond to cities in Europe and the 1417 edges correspond to roads linking the cities. We have also studied the Email network Guimera2003 of the University of Rovira i Virgili with 1133 vertices representing users and 5451 edges, each representing the existence of at least one Email communication between the two users corresponding to the vertices anchoring the edge. We have also studied the Route views network Leskovec2007 which has 6474 autonomous systems as vertices and 13895 edges representing communication between the systems that are represented as vertices. We have considered a social network, the Hamsterster friendship network Kunegis2013 , containing 1858 vertices which represent the users and 12534 edges which represent friendships between the users. Note that we omit self-loops while constructing the clique complex $K$ corresponding to the undirected graph $G$ of a real-world network.

Results and Discussion

Algorithm to construct discrete Morse function on a simplicial

complex

From an unweighted and undirected graph $G(\mathcal{V},\mathcal{E})$ with vertex set $\mathcal{V}$ and edge set $\mathcal{E}$, it is straightforward to construct a clique simplicial complex $K$ with dimension $d$ (See Theory section). Figure 1 shows the construction of a clique complex starting from an example network. Given a simplicial complex $K$, its dimension $d$ and a non-negative real-valued function $g$ on the 0-simplices of $K$, the algorithm 1 assigns weights to any simplex in $K$, producing a discrete Morse function $f$ defined in equation 6. In the pseudocode of the algorithm 1, lines 2-6 initialize a variable Flag$[\alpha]$ for every simplex $\alpha$ in clique complex $K$ with the value 0. We remark that the variable Flag$[\alpha]$ associated with a simplex $\alpha$ in $K$ serves as a counter for the size of the set $U_{\alpha}$ defined in equation 4. Lines 7-9 assign weights to every [math]-simplex in $K$ based on the input non-negative function $g$. Lines 10-24 assign weights to $1$- or higher-dimensional simplices in $K$ in a manner which is consistent with the definition in equation 6 of a discrete Morse function. In summary, algorithm 1 outputs a discrete Morse function $f$ on $K$, and in SI Appendix, we present a rigorous proof for the following theorem which states the same.

Theorem. *Algorithm 1 produces a discrete Morse function $f$ on any simplicial complex $K$ of finite dimension $d$.

Given a simplicial complex $K$, its dimension $d$ and a discrete Morse function $f$ on $K$, the algorithm 2 in SI Appendix determines the weights of critical simplices in $K$. Given an unweighted and undirected graph $G$, we restrict the construction of clique complex $K$ by including simplices up to a maximum dimension $d$. Then, the algorithm 3 in SI Appendix creates the filtration of clique complex $K$ based on weights of critical simplices as described in the Theory section. In SI Table S1, we describe the role of key variables which appear in algorithms 1, 2 and 3. In SI Appendix, we also give a time complexity analysis for these algorithms.

Rationale for the choice of function on vertices

In order to construct a discrete Morse function $f$ on clique complex $K$ corresponding to a graph $G$ using our algorithm 1, a real-valued function $g$ has to be fixed on the [math]-simplices of $K$ (See lines 7-9 in algorithm 1). Let $deg_{max}$ denote the maximum degree of a vertex in the graph $G(\mathcal{V},\mathcal{E})$. Our choice for the function value on the vertices or [math]-simplices, $g:\mathcal{V}\rightarrow\mathbb{R}$, is as follows:

[TABLE]

where degree$(v)$ is the degree of the vertex $v\in G$ and $\epsilon$ corresponding to each vertex is a random number (noise) generated using the uniform distribution on the interval $(0,0.5)$.

In the Theory section, we had highlighted the Theorem 2.11 by Forman Forman2002 which gives a lower bound on the number of critical $p$-simplices, $m_{p}$, in a simplicial complex $K$ as the $p$-Betti number $\beta_{p}$. The choice of the real-valued function $g$ in algorithm 1 plays a key role in determining if $m_{p}$ is close to the theoretical minimum $\beta_{p}$ stated above. In the Theory section, we have shown that the number of critical simplices determines the effective number of filtration weights to study the persistent homology of a clique complex (See equation 10). This motivated our choice for the real-valued function (equation 13) which determines the weights of [math]-simplices, and the rationale for this choice is as follows.

Ignoring the noise term $\epsilon$ in equation 13, the reader can discern our intuition for choosing the function $g(v)=deg_{max}-\text{degree}(v)$ for any vertex $v$ in $G$ with the following example. Consider the simple example of the clique complex $K$ corresponding to a graph $G$ in figure 1. Here, we would like to obtain a discrete Morse function $f$ on $K$ such that the number of critical simplices is close to the theoretical minimum. This requirement applies to simplices of any dimension in $K$, and in the context of this example, we would like the number of critical $1$-simplices (edges) to be as close as possible to the $1$-Betti number $\beta_{1}$ of $K$. Note that $\beta_{1}=1$ for the example clique complex $K$ in figure 1.

Let us now examine lines 11-23 in algorithm 1. Consider any edge $e_{vw}=[v,w]$ such that $g(v)>g(w)$. While assigning the function value to the edge $e_{vw}$ in algorithm 1, the edge $e_{vw}$ and the vertex $v$ are guaranteed to be not critical provided that the if condition in the line 16 is satisfied. This is a consequence of the definition of a critical simplex (See equation 7). Thus, we would like to force this if condition to be True for as many edges as possible. Moreover, once the function value of the $1$-simplex $e_{vw}$ is set, we set the variable Flag$[v]$ to 1 in line 18, and this subsequently forces the if condition in the line 16 to fail for all other edges $e_{vz}=[v,z]$ in the graph that contain $v$ and have function value $g(v)>g(z)$.

Let us now examine the edge $e_{78}=[v_{7},v_{8}]$ in figure 1 which is anchored by vertices $v_{7}$ and $v_{8}$ with degree 5 and 1, respectively. As the degree of a vertex gives the number of edges that contain the vertex, $v_{7}$ is part of 4 other edges apart from $e_{78}$ while $v_{8}$ is part of only the edge $e_{78}$. Suppose $e_{78}$ is the first edge chosen for the function assignment in line 11 of algorithm 1 and both Flag[$v_{7}$] and Flag[$v_{8}$] for the anchoring vertices are [math]. We would then prefer that the if condition in the line 16 is satisfied for $e_{78}$, and as $v_{7}$ is part of 4 other edges apart from $e_{78}$ while $v_{8}$ is part of only $e_{78}$, ideally Flag[$v_{8}$] is set to 1 instead of Flag[$v_{7}$], in other words, we need the function value $g(v_{8})>g(v_{7})$. We emphasize that this choice of $v_{8}$ over $v_{7}$ prevents the forced failure (described in the previous paragraph) of the if condition for the 4 other edges apart from $e_{78}$ that contain $v_{7}$.

The above example suggests a need for a function $g$ on the vertices that has an inverse relationship with the degree of the vertices. Hence, our choice $g(v)=deg_{max}-\text{degree}(v)$ provides a simple and effective solution for the above requirement. As a consequence of this choice for the function $g$, the weights assigned to the simplices by algorithm 1 reflect the degree of its constituent vertices. The degree of a vertex can be thought of as a measure of its importance in the network. Hence, the intuition behind the assignment of weights to simplices by our method is to have an inverse or opposite relationship between ‘weight’ and ‘importance’ of a simplex while simultaneously satisfying the definition of a discrete Morse function. This inverse relationship instead of a proportional relationship between ‘importance’ and ‘weight’ plays a key role in our filtration scheme which is based on the sequence of level subcomplexes corresponding to the increasing sequence of critical weights (See Theory section and algorithm 3 in SI Appendix). Thus, within the constraint of ensuring that the definition for a simplicial complex is satisfied at each stage of the filtration, our scheme prioritizes the addition of simplices with higher ‘importance’ at earlier stages of the filtration due to their lower weights.

We now provide a rationale for the addition of a random noise $\epsilon$ in equation 13. As reasoned above, we would like to force the if condition in line 16 of algorithm 1 to be True for as many edges as possible. Consider an edge $e_{vw}=[v,w]$ such that $degree(v)=degree(w)$. The absence of a random noise $\epsilon$ in equation 13 forces $g(v)=g(w)$. Thus, irrespective of the state of Flag$[v]$ and Flag$[w]$, the if condition fails. This implies that the set $V_{e_{vw}}$ (See equation 5) corresponding to the edge $e_{vw}$ is forced to be empty since the only other possibility is $|V_{e_{vw}}|=2$ which cannot be the case because algorithm 1 produces a discrete Morse function (See Theorem and equation 6). Thus, provided $e_{vw}$ is not a face of any higher dimensional simplex, $e_{vw}$ would be a critical simplex irrespective of the state of Flag$[v]$ and Flag$[w]$. Hence, we would like $g(v)\neq g(w)$ while also retaining the inverse relationship of the function with the degree. Thus generating a small random noise $\epsilon$ in the range $(0,0.5)$ for each vertex as in equation 13 provides a simple resolution. We remark that the above argument can be generalized to higher-dimensional simplices, and thus, provides the intuition for the addition of noise $\epsilon$ in line 21 of algorithm 1.

We remind the readers that our initial motivation was not to develop a scheme to construct the optimal discrete Morse function on a clique complex corresponding to a graph. Rather, our main goal is to develop a systematic filtration scheme to study persistent homology in unweighted and undirected networks. In fact, constructing an optimal discrete Morse function in the general case has been shown to be MAX-SNP Hard Lewiner2003 . The primary utility of our scheme is to create a filtration by assigning weights to simplices in the clique complex $K$ of a graph $G$. However, we next report our empirical results from an exploration of model and real-world networks which underscore the following. Although our scheme is not optimal in the sense of minimizing the number of critical simplices, in practice, it achieves near-optimal results in several model and real-world networks (Table 1). Hence, our scheme based on discrete Morse theory reduces the number of filtration steps and increases the applicability of persistent homology to study complex networks.

Application to model and real networks

Given a model or real network $G$, we limit our study of persistent homology to the $3$-dimensional clique simplicial complex $K$ corresponding to $G$. In other words, during the construction of the clique complex, we only include $p$-simplices which have dimension $0\leq p\leq 3$ (See Theory section). On this $3$-dimensional clique complex $K$, we create the corresponding filtration based on the assigned weights to simplices using discrete Morse theory. In SI Appendix, we present algorithm 3 which outlines the procedure to compute the filtration weights of simplices. Thereafter, we make use of GUDHI Maria2014 , a C++ based library for Topological Data Analysis (http://gudhi.gforge.inria.fr/) to study the persistent homology of this filtration of $K$.

Each hole of any dimension in $K$ has a corresponding birth filtration weight and a death filtration weight (See Theory section). We normalize the birth and death filtration weights of all holes in $K$ by dividing with $w_{N}=1+\text{max}\{f(\alpha)\>|\>\alpha\in K\}$. In other words, $w_{N}$ is $1$ plus the maximum value among weights assigned to the simplices in $K$. We also make the convention that a normalized death filtration weight of $1$ for a hole in $K$ represents that the particular hole never dies.

A $H_{p}$-barcode diagram corresponding to a filtration of the clique complex $K$ is a graphical representation containing horizontal line segments, each of which represents a $p$-hole in $K$, plotted against the $x$-axis ranging from [math] to $1$ which corresponds to the normalized filtration weights of simplices in $K$ Ghrist2008 . A horizontal line in the $H_{p}$-barcode diagram of $K$ is referred to as a barcode. Thus, a barcode in the $H_{p}$-barcode diagram of $K$ which begins at a $x$-axis value of $w_{1}$ and ends at a $x$-axis value of $w_{2}$ represents a $p$-hole in $K$ whose birth and death weights are $w_{1}$ and $w_{2}$, respectively. In figures 3, 4 and 5, we display the barcode diagrams for model and real-world networks analyzed here.

Model Networks

In this work, we have investigated the persistent homology of unweighted and undirected graphs corresponding to five model networks, namely, ER, WS, BA, hyperbolic random graphs and spherical random graphs. We have considered model networks with 1000 vertices and expected average degree 4, 6 and 8. In main text, we report results for model networks with expected average degree 4, and in SI, those with expected average degree 6 and 8.

The $H_{0}$ barcode diagram of BA networks indicate a low number of [math]-holes in BA networks across the entire filtration (Figure 3 and SI Figures S1-S2). A standard result in algebraic topology Munkres2018 gives that the [math]-Betti number of a simplicial complex $K$ is equal to the number of connected components in $K$. In other words, the above observation indicates that the scale-free BA network has a strong tendency to maintain a low number of connected components during filtration (Figure 3 and SI Figures S1-S2). In contrast, both the random ER network and small-world WS network have a relatively high number of connected components at initial phase of the filtration and then progress towards a more connected network at later stages of the filtration (Figure 3 and SI Figures S1-S2). This indicates that the simplices in the clique complex which are key to the connectivity of the model networks are introduced very early into the filtration for the scale-free BA network while this is not the case for the random ER network or small-world WS network. The $H_{1}$ barcode diagram of BA networks also indicate late introduction of $1$-holes during filtration in contrast to both ER and WS networks where $1$-holes appear across a wider range of the filtration (Figure 3 and SI Figures S1-S2). Moreover, in ER and WS networks, it is interesting to observe that the $1$-holes start to appear at roughly the same stage of filtration which corresponds to a sharp reduction in the number of connected components (Figure 3 and SI Figures S1-S2).

In contrast to ER, WS and BA networks, the spherical and hyperbolic networks are characterized by a relatively high [math]-Betti number $\beta_{0}$ and a low $1$-Betti number $\beta_{1}$ (Table 1 and SI Table S2). Simply stated, this observation on the magnitude of $\beta_{0}$ indicates that both hyperbolic and spherical networks have a higher number of connected components in comparison to the ER, WS and BA networks of similar size, i.e., number of vertices, and average vertex degree (Table 1 and SI Table S2). Although, both spherical and hyperbolic networks exhibit a higher number of connected components, they differ from each other with respect to the evolution of these connected components during filtration (Figure 3 and SI Figures S1-S2). The hyperbolic model maintains a relatively low number of connected components until very late in the filtration wherein there is a sharp increase in the number of connected components (Figure 3 and SI Figures S1-S2). This is in contrast with the behavior of the $H_{0}$ barcode diagram of the spherical model which exhibits a more distributed evolution of connected components during filtration (Figure 3 and SI Figures S1-S2). In addition, the low $\beta_{1}$ for spherical and hyperbolic networks conveys the lack of $1$-holes in these networks. A possible reason for this observation is the incidence of higher number of $2$-simplices in the clique complex $K$ of the spherical and hyperbolic networks in comparison to the ER, WS and BA networks (Table 1 and SI Table S2). Note that the formation of a $2$-simplex can potentially fill in a $1$-hole, and thus, result in a low value for $\beta_{1}$ (See Theory section). Such a behaviour is also seen in the $H_{2}$ barcode diagrams of spherical and hyperbolic networks wherein the $2$-holes have very short persistence since the addition of $3$-simplices successively fill in the $2$-holes (Figure 5 and SI Figure S3). The $H_{3}$ barcode diagrams of spherical and hyperbolic networks (See SI Figures S3-S4) also indicate a clear difference in the evolution of their corresponding topological features during the filtration. The $H_{2}$ and $H_{3}$ barcode diagrams of ER, WS and BA networks do not provide any insight into network structure primarily due to a lack of higher-order correlations in these model networks that is essential for the formation of $2$-holes and $3$-holes.

A visual inspection of the barcode diagrams for the five model networks (Figures 3, 5 and SI Figures S1-S4) suggests that the different models can be distinguished based on their persistent homology. In Theory section, we had introduced the bottleneck distance which can be employed to quantify the differences between persistence diagrams from filtration of clique complexes corresponding to different model networks. Recall that the persistence diagram of a $d$-dimensional simplicial complex $K$ is a compact representation of the persistent homology of $K$ which encompasses topological information across all $d$ dimensions. Figure 6 and SI Table S3 give the bottleneck distance between different model networks with the same number of vertices and similar average vertex degree. For each of the five model networks, 10 random samples are generated by fixing the number of vertices $n$ and other parameters of the model. In figure 6 and SI Table S3, we report the distance between two different models as the average of the distance between each of the possible pairs of the 10 sample networks corresponding to the two models along with the standard error. We find a relatively higher distance between a random instance of a BA network and a random instance of a ER network with the same number of vertices and similar average degree (Figure 6 and SI Table S3). Similarly, we observe a relatively higher distance between a random instance of a BA network and a random instance of a WS network with similar size and average degree. In contrast, a relatively lower average distance is observed between a random instance of a ER network and a random instance of a WS network with similar size and average degree (Figure 6 and SI Table S3). These observations indicate a similarity between networks generated by ER and WS models in terms of their corresponding persistence diagrams and also show that the BA model exhibits topological properties that are different from the other two model networks with the same number of vertices and similar average degree. Moreover, the average distance between a random instance of a spherical network and a random instance of a hyperbolic network with similar size and average degree is very high (Figure 6 and SI Table S3). The last observation is a reflection of the differences in the persistent homology of the clique complexes corresponding to spherical and hyperbolic networks. Finally, the nature of differences observed between persistence diagrams of different model networks using bottleneck distance as shown in figure 6 remain consistent if the 1-Wasserstein or 2-Wasserstein distance is employed in place of the bottleneck distance (data not shown).

Table 1 and SI Table S2 list the empirical data on the number of $p$-simplices $n_{p}$, the number of critical $p$-simplices $m_{p}$, that our algorithm achieves and the $p$-Betti number $\beta_{p}$ of the clique complexes corresponding to the five model networks. In Table 2 we report a value $\mu$ which indicates the optimality of our algorithm with respect to reducing the number of critical simplicies for each model network that has been analyzed. The definition for $\mu$ is as follows.

[TABLE]

Here $d$ is the dimension of the corresponding clique complex. The value $\mu$ corresponding to a particular discrete Morse function on the clique complex of a given network is an indicator of its optimality with respect to minimizing the number of critical simplices. The value of $\mu$ ranges from [math] to $1$. $\mu=1$ indicates that the discrete Morse function achieves exactly the minimum number of critical simplices thereby corresponding to the most optimal situation, while $\mu=0$ indicates that all simplices are critical thereby corresponding to the least optimal case. Note that the value of $\mu$ for a particular network increases linearly with a decrease in the number of critical simplices. In Table 2, for model networks, the value reported is the average of $\mu$ across 10 samples of each model network for a chosen set of parameter values along with the corresponding standard deviations. Based on the data presented in Table 2, we report that our algorithm achieves near-optimal results in terms of reducing the number of critical simplices for each of the five model networks analyzed. Moreover for these model networks, data presented in SI Table S2 underlines the close proximity of $m_{p}$ to the theoretical lower bound $\beta_{p}$ across all dimensions $p=0,1,2,3$ (See Theory section). We remark that the worst case scenario corresponds to all simplices in the clique complex being critical which is never attained by our algorithm in the analyzed networks.

Comparison of our method with dimension based filtration in model networks

In SI Figures S5 and S6, we show the $H_{0}$, $H_{1}$, $H_{2}$ and $H_{3}$ barcode diagrams for model networks with expected average degree 4 obtained using the dimensional filtration scheme used in Horak et al. Horak2009 We restrict our investigation to the three-dimensional clique complex while computing the barcode diagrams for model networks using the dimensional filtration scheme of Horak et al. Horak2009 In SI Figures S5 and S6, we normalize the filtration index to be in the range 0 to 1, and $p$-holes with normalized filtration index 1 indicate that they never die. Moreover, we also report in SI Table S4 and SI Figure S7, the bottleneck distances between persistence diagrams of model networks obtained by the dimensional filtration scheme used in Horak et al. Horak2009 , and we find that the resultant barcode diagrams and the bottleneck distances between the persistence diagrams are inconclusive to distinguish between the five model networks which have similar size and average vertex degree.

The discrete Morse function produced by our algorithm is empirically near-optimal for the model networks analyzed here, in the sense of minimizing the number of critical simplices (See Table 2). On the other hand, it should be noted that the dimension function used in Horak et al. Horak2009 is the most non-optimal in this regard. Moreover, in terms of persistent homology, our methods based on the discrete Morse function constructed using algorithms 1, 2 and 3 have the distinct feature of being able to distinguish between various model networks with an efficient filtration.

Although the dimension function used in Horak et al. Horak2009 is algorithmically efficient in terms of having a low number of filtration steps, unlike our method it does not conclusively distinguish between different model networks based on the barcode diagrams and also the bottleneck distances between the corresponding persistence diagrams. Thus our discrete Morse function is a good candidate for applications in both, computational aspects of discrete Morse theory as well as in persistent homology of unweighted networks, and achieves a tradeoff between efficiency and applicability.

Real networks

In this work, we have investigated the persistent homology of unweighted and undirected graphs corresponding to seven real-world networks and the barcode diagrams for five of these real networks is shown in figures 4, 5 and SI figure S4. Based on the $H_{0}$ barcode diagrams, the behavior of most real networks considered here can be broadly classified into two categories. Real networks such as the Email communication, the Hamsterster friendship and the Route views exhibit a relatively low number of connected components across the entire range of filtration (Figure 4). On the other hand, the two biological networks, namely, the Yeast protein interaction and the Human protein interaction, exhibit a sharp increase in the number of connected components at the later stages of filtration. The $H_{0}$ barcode diagrams for the US Power Grid and Euro road do not conform with either of the above characterizations. The $H_{0}$ barcode diagram of the US Power Grid network reveals that though there exists only a single connected component at the end of the filtration, there are a considerable number of non-persisting connected components that appear and subsequently disappear during filtration (Figure 4). The $H_{0}$ barcode diagram of the Euro road network shows a more distributed increase in the number of connected components (data not shown).

In the context of the $H_{1}$ barcode diagrams, the real networks considered here exhibit similar properties with $1$-holes appearing late in the filtration (Figure 4). The $H_{2}$ barcode diagrams reveal a lack of $2$-holes with long persistence in both biological networks, as well as the Route views network, the Euro road network and the US Power Grid network (Figure 5). In contrast, from the $H_{2}$ and $H_{3}$ barcode diagrams (Figure 5, SI Figure S4) we find that the social network, Hamsterster friendship, and the Email communication network exhibit a relatively high number of $2$-holes and $3$-holes with longer persistence.

Table 2 lists the value of the optimality indicator $\mu$ (See Equation 14) for each of the seven real networks analyzed here. This data indicates near-optimal performance of our algorithm with respect to minimizing the number of critical simplices for each of these seven real networks. Table 1 also lists the empirical data on the number of critical $p$-simplices, $m_{p}$, that our algorithm achieves and the $p$-Betti number $\beta_{p}$ of the clique complexes across each dimension $p$, corresponding to the seven real networks analyzed here.

Conclusions

To conclude, we have proposed a systematic scheme based on discrete Morse theory to study the persistent homology of unweighted and undirected networks. Our methods leverage the concept of critical simplices to permit a reduced filtration scheme while simultaneously admitting a finer inspection of the changes in topology across the filtration of a clique complex corresponding to an unweighted network. Moreover, our proposed algorithm to construct a discrete Morse function on the clique complex of a simple graph achieves close to optimal number of critical simplices for several model and real networks that have been studied here. Furthermore, based on visual representations of persistent homology such as the barcode diagrams as well as quantitative information in the form of distance between persistence diagrams, our methods successfully distinguish various model networks that exhibit inherently different properties. This motivates the application of our methods to real-world networks. We report the results obtained for seven real-world networks that are well studied in the network science community and observe certain patterns in the evolution of their topological features across the filtration. For instance both biological networks, namely the Yeast protein interaction network and the Human protein interaction network exhibit similar characteristics with respect to the $H_{0}$, $H_{1}$, $H_{2}$ and $H_{3}$ barcode diagrams. Similarly, both the Email network and the Hamsterster friendship network, exhibit shared features with respect to $H_{0}$, $H_{1}$, $H_{2}$ and $H_{3}$ barcode diagrams that vary from the characteristics of the two biological networks considered here. Our observations hint at the ability and possible applications of our methods to detect and classify real-world networks that are inherently different.

Future directions and ongoing work include examining the significance of critical simplices in the context of real-world networks. In other words, we aim to determine whether a critical edge in the context of discrete Morse theory holds any key significance when it is viewed as a link between two real entities in a real-world network. We also intend to explore the presence or absence of a correlation between the notion of critical simplices and network curvature. Since discrete Morse theory captures information about the Euler characteristic of the clique complex corresponding to a graph, the presence of such a correlation could potentially signify a close relationship between the discrete curvature of a graph and its topology, much like in the case of smooth, compact surfaces wherein the Gauss-Bonnet theorem relates the Gaussian curvature of a surface to its Euler characteristic.

Acknowledgments

We thank Amritanshu Prasad for fruitful discussions. A.S. would like to acknowledge support from the Max Planck Society, Germany, through the award of a Max Planck Partner Group in Mathematical Biology, and I.R. from the Science and Engineering Research Board (SERB) of the Department of Science and Technology (DST) India through the award of a MATRICS grant [MTR/2017/000835]. We thank the anonymous reviewers for their comments which have helped improve this manuscript.

Author contributions

H.K., I.R. and A.S. designed the study. H.K. performed the simulations. H.K., E.S., I.R. and A.S. analyzed results. H.K., I.R. and A.S. wrote the manuscript. All authors reviewed and approved the manuscript.

Homology of a simplicial complex

In this section, we provide an overview of the mathematical theory of homology groups of a simplicial complex.

We first describe a mathematical group which provides the machinery to represent paths in a simplicial complex. The $p$-chain group $C_{p}$ of a simplicial complex $K$ is the Abelian group generated by the oriented $p$-simplices in $K$ with coefficients in a field $\mathbb{F}$. Elements $c$ of the $p$-chain group $C_{p}$ are referred to as $p$-chains and have the form $c=\sum_{i}n_{i}\alpha_{i}$ where we use the same notation to represent both the oriented $p$-simplex and its corresponding generator in $C_{p}$ and $n_{i}$ are scalars from the field $\mathbb{F}$. Note that the identity element [math] in $C_{p}$ is the unique $p$-chain for which all the coefficients $n_{i}$ are zero in $\mathbb{F}$. Note also that if two $p$-simplices, $\alpha$ and $\beta$, have the same vertex set but opposite orientations, then the generators of $C_{p}$ corresponding to the two simplices are inverse of each other (i.e., $\alpha=-\beta$). In figure 1 of the main text, the cycle of length 4 formed by the edges or 1-simplices, $[v_{1},v_{2}]$, $[v_{2},v_{5}]$, $[v_{5},v_{7}]$ and $[v_{1},v_{7}]$, can be represented as an element of the $1$-chain group which is $c=[v_{1},v_{2}]+[v_{2},v_{5}]+[v_{5},v_{7}]-[v_{1},v_{7}]$.

The boundary operator $\partial_{p}$ on the generator corresponding to an oriented $p$-simplex $\alpha=[v_{0},v_{1},\cdots,v_{p}]$ is defined as follows Munkres2018 :

[TABLE]

where $\hat{v_{i}}$ refers to the absence of $v_{i}$ in the $(p-1)$-simplex. By linear extension, the boundary operator $\partial_{p}$ on any element $c=\sum_{i}n_{i}\alpha_{i}$ of $p$-chain group $C_{p}$ gives:

[TABLE]

Note that the boundary operator in the above equation maps a $p$-chain to a $(p-1)$-chain. In figure 1 of the main text, the boundary operator applied to the $1$-chain $c=[v_{0},v_{1}]$ representing an edge gives $v_{1}-v_{0}$ while the boundary operator applied to the $1$-chain $c=[v_{1},v_{2}]+[v_{2},v_{5}]+[v_{5},v_{7}]-[v_{1},v_{7}]$ representing a cycle of length 4 gives 0.

This motivates the definition of $p$-cycles $Z_{p}$ and $p$-boundaries $B_{p}$. The $p$-cycles $Z_{p}$ are the elements of the $p$-chain group $C_{p}$ which are mapped to 0 by the boundary operator $\partial_{p}$, and thus:

[TABLE]

The $p$-boundaries $B_{p}$ is defined as follows:

[TABLE]

Thus, $p$-boundaries $B_{p}$ are the elements of the $p$-chain group $C_{p}$ which also happen to be the boundary of an element in the $(p+1)$-chain group $C_{p+1}$. Note that both $Z_{p}$ and $B_{p}$ are subgroups of the $p$-chain group $C_{p}$. It can be shown that the composition $\partial_{p}\circ\partial_{p+1}=0$. This implies that $B_{p}$ is a subgroup of $Z_{p}$ Munkres2018 . Simply stated, the $p$-boundary operator $\partial_{p}$ when applied on a $p$-boundary gives 0. Hence, it follows that every $p$-boundary is a $p$-cycle but not necessarily vice versa.

The $p$-homology group is defined as Munkres2018 :

[TABLE]

where $Z_{p}/B_{p}$ is the quotient group Dummit2003 of $Z_{p}$ over $B_{p}$. We informally refer to the elements of $p$-homology group $H_{p}$ as $p$-holes. To provide an intuition for the definition of homology groups, a natural way to describe a $p$-hole would be to characterize it as a $p$-cycle which is not a $p$-boundary. In figure 1 of the main text, the $1$-cycle $c=[v_{1},v_{2}]+[v_{2},v_{5}]+[v_{5},v_{7}]-[v_{1},v_{7}]$ is a $1$-hole as it is not a $1$-boundary of any $2$-chain while the $1$-cycle $c=[v_{2},v_{3}]+[v_{3},v_{4}]+[v_{4},v_{2}]$ is not a $1$-hole as it is the $1$-boundary of the $2$-simplex $[v_{2},v_{3},v_{4}]$. Thus, the concept of quotient groups provide the mathematical machinery to characterize such $p$-holes in simplicial complexes.

We remark that though it is customary to call $H_{p}$ a group, since we use coefficients from a field $\mathbb{F}$, it satisfies additional properties which make it a vector space over $\mathbb{F}$. The $p$-Betti number is defined as the dimension of the homology group $H_{p}$ viewed as a vector space over field $\mathbb{F}$ Munkres2018 . Informally, the $p$-Betti number $\beta_{p}$ represents the number of $p$-holes of the simplicial complex. We remark that the Euler characteristic of the clique complex $K$ with dimension $d$ corresponding to a graph $G$ is given by the alternating sum of Betti numbers Munkres2018 , namely,

[TABLE]

In this work, we use the finite field $\mathbb{F}=\mathbb{Z}/2$, i.e., the field with two elements.

Proof of concept for Algorithm 1

In this section, we present a rigorous proof for the following theorem.

Theorem 1. *Algorithm 1 produces a discrete Morse function $f$ on any simplicial complex $K$ of finite dimension $d$.

Proof. Let $K$ denote a simplicial complex of dimension $d$. Recall that for the function $f$ on $K$ which is constructed by algorithm 1, for each simplex $\alpha^{p}\in K$, the two sets $U_{\alpha}$ and $V_{\alpha}$ are defined as follows:

[TABLE]

To prove that $f$ is a discrete Morse function we need to show that for each simplex $\alpha\in K$, both $|V_{\alpha}|\leq 1$ and $|U_{\alpha}|\leq 1$.

Firstly, we show that for each simplex $\alpha\in K$, $|V_{\alpha}|\leq 1$. Consider a [math]-simplex $\alpha\in K$. Since, the dimension of a simplex cannot be less than [math], for each [math]-simplex $\alpha\in K$ the corresponding set $V_{\alpha}$ is empty. In other words, for each [math]-simplex $\alpha\in K$, we have that $|V_{\alpha}|=0$. Also, for each $p$-simplex $\alpha^{p}\in K$ such that $1\leq p\leq d$, Lemma 6 below shows that $|V_{\alpha}|\leq 1$. Thus, for every simplex $\alpha\in K$, we have shown that $|V_{\alpha}|\leq 1$.

Secondly, we show that for each simplex $\alpha\in K$, $|U_{\alpha}|\leq 1$. For each $p$-simplex $\alpha^{p}\in K$ such that $0\leq p\leq(d-1)$, we prove in Lemma 5 below that $|U_{\alpha}|\leq 1$. Now consider a $d$-simplex $\alpha^{d}\in K$. Since, by assumption $K$ is a $d$-dimensional simplicial complex, there are no $(d+1)$-simplices in $K$, and thus, the set $U_{\alpha}$ for each $d$-simplex $\alpha^{d}$ in $K$ is empty. In other words, for each $d$-simplex $\alpha^{d}\in K$, we have that $|U_{\alpha}|=0$. Thus, for every simplex $\alpha\in K$ we have shown that $|U_{\alpha}|\leq 1$.

In summary, we have shown that for each simplex $\alpha\in K$, $|V_{\alpha}|\leq 1$ and $|U_{\alpha}|\leq 1$. Thus, $f$ satisfies the definition of a discrete Morse function on the simplicial complex $K$.

∎

We next prove the Lemmas used in the proof of the theorem above. This is done in the following sequence of Lemmas 2 to 6. We assume that $K$ is a $d$-dimensional simplicial complex and $f$ is the output function on $K$ obtained from algorithm 1. We remark that a $p$-dimensional simplex $\alpha$ of a simplicial complex $K$ is denoted by $\alpha^{p}\in K$. Also, if $p$-simplex $\alpha$ is a face of a $(p+1)$-simplex $\beta$ then this is represented as $\alpha^{p}<\beta^{p+1}$ in the sequel.

Lemma 2. *For each $p$ where $0\leq p\leq(d-1)$, if $\alpha^{p}\in K$ and $\beta^{p+1}\in K$ such that $\alpha^{p}<\beta^{p+1}$, then, $f(\alpha)\neq f(\beta)$.

Proof. Let $\gamma_{0},\gamma_{1},\cdots,\gamma_{p+1}$ denote the $p$-dimensional faces of $\beta^{p+1}$ such that $f(\gamma_{0})\geq f(\gamma_{1})\geq\cdots\geq f(\gamma_{p+1})$. Note that $\alpha^{p}$ is one such $\gamma_{i}$ since by assumption $\alpha^{p}$ is a $p$-dimensional face of $\beta^{p+1}$. Based on lines 11-23 in algorithm 1, we have that:

[TABLE]

where $\epsilon>0$.

Case $A$ implies $f(\gamma_{0})>f(\beta)>f(\gamma_{1})\geq f(\gamma_{2})\cdots\geq f(\gamma_{p})\geq f(\gamma_{p+1})$.

Case $B$ implies $f(\beta)>f(\gamma_{0})\geq f(\gamma_{1})\geq f(\gamma_{2})\cdots\geq f(\gamma_{p})\geq f(\gamma_{p+1})$.

Thus, for both cases we have that $f(\beta)\neq f(\gamma_{i})$ for each $i\in\{0,1,2,\cdots,(p+1)\}$. Since $\alpha^{p}$ is one such $\gamma_{i}$ for some $i\in\{0,1,2,\cdots,(p+1)\}$, we have $f(\alpha)\neq f(\beta)$.

∎

Lemma 3. For each $p$ where $0\leq p\leq(d-1)$, if $\alpha^{p}\in K$ and $\beta^{p+1}\in K$ such that $\alpha^{p}<\beta^{p+1}$, then $f(\alpha)>f(\beta)$ if and only if Flag$[\alpha]$ *changes value from [math] to $1$ while assigning function value for $\beta$.

Proof. Given $\alpha^{p}\in K$ and $\beta^{p+1}\in K$  such that  $\alpha^{p}<\beta^{p+1}$, we first assume $f(\alpha)>f(\beta)$. Let $\gamma_{0},\gamma_{1},\cdots,\gamma_{p+1}$ denote the $p$-dimensional faces of $\beta^{p+1}$ such that $f(\gamma_{0})\geq f(\gamma_{1})\geq\cdots\geq f(\gamma_{p+1})$. Then, the value $f(\beta)$ is given by equation A.7.

Case $A$ in equation A.7 implies $f(\gamma_{0})>f(\beta)>f(\gamma_{1})\geq f(\gamma_{2})\cdots\geq f(\gamma_{p})\geq f(\gamma_{p+1})$.

Case $B$ in equation A.7 implies $f(\beta)>f(\gamma_{0})\geq f(\gamma_{1})\geq f(\gamma_{2})\cdots\geq f(\gamma_{p})\geq f(\gamma_{p+1})$.

Since by assumption, $\alpha^{p}$ is a face of $\beta^{p+1}$ and $f(\alpha)>f(\beta)$, Case $A$ is applicable, and we have $\gamma_{0}$ equals $\alpha^{p}$. Thus, based on line 18 in algorithm 1, while assigning the function value on $\beta^{p+1}$, Flag$[\alpha]$ changes value from [math] to $1$.

Now, given $\alpha^{p}\in K$ and $\beta^{p+1}\in K$  such that  $\alpha^{p}<\beta^{p+1}$, we assume that Flag$[\alpha]$ changes value from [math] to $1$ while assigning function value for $\beta^{p+1}$. Let $\gamma_{0},\gamma_{1},\cdots,\gamma_{p+1}$ denote the $p$-dimensional faces of $\beta^{p+1}$ such that $f(\gamma_{0})\geq f(\gamma_{1})\cdots\geq f(\gamma_{p+1})$. Based on lines 11-23 in algorithm 1, Flag$[\alpha]$ changes value from [math] to $1$ while assigning function value on $\beta^{p+1}$ implies that $\gamma_{0}$ equals $\alpha^{p}$. Thus, we have that $f(\alpha)>f(\beta)=(f(\alpha)+f(\gamma_{1}))/2$.

∎

Lemma 4. Let $\alpha$ be a simplex of $K$. Then, the number of times Flag$[\alpha]$ *changes value from [math] to $1$ is $\leq 1$.

Proof. Flag$[\alpha]$ is initially set to [math] in algorithm 1. In algorithm 1, if Flag$[\alpha]$ transitions to 1, its value never changes. In other words, there is no procedure in our algorithm 1 which changes the Flag of a simplex from 1 to 0. Thus, either Flag$[\alpha]$ remains 0 throughout or changes value from 0 to 1 exactly once in algorithm 1.

∎

Lemma 5. *For each $p$ where $0\leq p\leq(d-1)$, if $\alpha^{p}\in K$, then $|U_{\alpha}|\leq 1$.

Proof. From Lemma 3, we have that for each $\alpha^{p}\in K$, the number of $(p+1)$-simplices $\beta^{p+1}$ such that $\alpha^{p}<\beta^{p+1}$ and $f(\beta)<f(\alpha)$ is equal to the number of times Flag$[\alpha]$ changes value from [math] to $1$. Applying Lemma 4, we get that for each $\alpha^{p}\in K$,

[TABLE]

Furthermore, Lemma 2 tells us that,

[TABLE]

Thus, we have that for each $\alpha^{p}\in K,\ |U_{\alpha}|$ = #$\{\beta^{p+1}\>|\>\alpha^{p}<\beta^{p+1}\text{ and }f(\beta)\leq f(\alpha)\}\;\leq\;1$.

∎

Lemma 6. *For each $p$ where $1\leq p\leq d$, if $\alpha^{p}\in K$, then $|V_{\alpha}|\leq 1$.

Proof. Let $\gamma_{0},\gamma_{1},\cdots,\gamma_{p}$ denote the $(p-1)$-dimensional faces of $\alpha^{p}\in K$ such that $f(\gamma_{0})\geq f(\gamma_{1})\geq\cdots\geq f(\gamma_{p})$. Based on lines 11-23 in algorithm 1, we have that:

[TABLE]

where $\epsilon>0$.

Case $A$ implies $f(\gamma_{0})>f(\alpha)>f(\gamma_{1})\geq f(\gamma_{2})....\geq f(\gamma_{p-1})\geq f(\gamma_{p})$, and thus, $|V_{\alpha}|=1$.

Case $B$ implies $f(\alpha)>f(\gamma_{0})\geq f(\gamma_{1})\geq f(\gamma_{2})....\geq f(\gamma_{p-1})\geq f(\gamma_{p})$, and thus, $|V_{\alpha}|=0$.

Hence, for each $\alpha^{p}\in K$ with $1\leq p\leq d$, we have that $|V_{\alpha}|\leq 1$.

∎

Filtration of the clique complex based on weights of critical simplices

Given a simplicial complex $K$, its dimension $d$ and a discrete Morse function $f$ on $K$, the algorithm 2 determines the weights of critical simplices in $K$. In the pseudocode of the algorithm 2, lines 2-6 initialize a variable IsCritical$[\alpha]$ associated to every simplex $\alpha$ in clique complex $K$ to be True. Lines 7-17 determine the critical simplices in $K$ by checking for the condition in equation 7 of the main text which defines a critical simplex. Lines 18-28 determine the weights of critical simplices or critical weights in $K$. Finally, the algorithm 2 outputs an array $w_{c}[]$ which contains an increasing sequence of critical weights in $K$. Subsequently, this increasing sequence of critical weights will be used for the filtration of the clique complex $K$.

Given an unweighted and undirected graph $G$, we restrict the construction of clique complex $K$ by including simplices up to a maximum dimension $d$. Then, the algorithm 3 creates the filtration of clique complex $K$ based on weights of critical simplices as described in the Theory section of the main text. In the pseudocode of the algorithm 3, lines 2-6 assigns a non-negative function $g$ to [math]-simplices in clique complex $K$. Line 7 calls the algorithm 1 for the assignment of weights satisfying discrete Morse function to every simplex in $K$. Line 8 calls the algorithm 2 to obtain an increasing sequence of unique weights corresponding to critical simplices in $K$. Lines 9-11 initialize a variable IsAdded$[\alpha]$ associated to every simplex $\alpha$ in $K$ which tracks whether the simplex $\alpha$ has been added to the filtration or not. Lines 12-31 compute the filtration weight of each simplex $\alpha$ in $K$ as described in the Theory section of the main text.

In SI Table S1, we describe the role of key variables which appear in algorithms 2 and 3.

Computational aspects of our algorithms

Time complexity of Algorithm 1

We briefly discuss here some computational aspects of the algorithm 1 in the main text. Given a simplicial complex $K$, its dimension $d$ and a non-negative real-valued function $g$ on the 0-simplices of $K$, algorithm 1 assigns weights to any simplex in $K$, producing a discrete Morse function $f$ on $K$.

We remark that the algorithms 1, 2 and 3 were implemented using the C++ programming language. The simplices of a $d$-dimensional simplicial complex were stored using the set container which is a part of the C++ Standard Template Library (STL). The association of the discrete Morse function values to the corresponding simplices was implemented using the map container of the C++ STL. Implementing the storage of the Flag variable associated with each simplex was also done using the map container of C++ STL. We remark that accessing a particular value of a map container has a complexity which is logarithmic with respect to the size of the container, namely, the number of simplices $n$.

Assuming that the dimension $d$ of $K$ is constant, we here present a brief discussion about the complexity of algorithm 1 with respect to the number of simplices $n$ in $K$. We remark that since the dimension $d$ is constant and the complexity of accessing values using map containers is logarithmic with respect to $n$, lines 2-9 in algorithm 1 have a complexity of $O(n\log{}n)$. The complexity of the operation corresponding to line 12 in algorithm 1, namely, finding the $(p-1)$-dimensional faces of a $p$-simplex is independent of $n$. Such an operation depends on $p$ which has an upper bound of $d$ that has been assumed to be a constant.

The operation in line 13 of algorithm 1 corresponds to arranging the simplices computed in the preceding line based on their function values. As noted earlier, accessing values using map containers is logarithmic with respect to $n$. Thus, since $p\leq d$ and $d$ is fixed, the operation represented by line 13 of algorithm 1 has a complexity of $O(\log{}n)$. Lines 14-22 of algorithm 1 except the random number generation in line 20 of algorithm 1 require accessing values of map containers. Random number generation based on a uniform distribution using C++ can be implemented with amortized constant complexity. Thus, the block of pseudocode represented by lines 14-22 in algorithm 1 have complexity of $O(\log{}n)$. Since $d$ is constant and for any $p$ the number of $p$-simplices is bounded by the total number of simplices $n$, the pseudocode represented in lines 10-24 of algorithm 1 has a complexity of $O(n\log{}n)$. Hence the complexity of algorithm 1 with respect to the number of simplices $n$ can be described as $O(n\log{}n)$.

Time complexity of algorithms 1, 2 and 3 with respect to the number of critical simplices

Firstly, it is apparent that the complexity of algorithm 1 is independent of the number of critical simplices since the notion of critical simplices is relevant only after the discrete Morse function has been constructed.

Algorithm 2 accepts as inputs a simplicial complex $K$, its dimension $d$ and a discrete Morse function $f$ on $K$ and outputs the sorted sequence of weights which correspond to the critical simplices in $K$. Let $m$ denote the total number of critical simplices in $K$ corresponding to $f$. We assume henceforth that the simplicial complex $K$ and its dimension $d$ are fixed. This implies that the number of simplices $n$ and the dimension $d$ are treated as constants. In such a situation, the complexity of all lines in the pseudocode of algorithm 2 except for line 28 is constant and independent of $m$. Line 28 of algorithm 2 represents sorting the array which contains the function values of all critical simplices, i.e., the critical weights, and removing any duplicate entries. The number of critical weights is less than or equal to the number of critical simplices $m$, and thus, the operation represented by line 28 of algorithm 2 can be implemented with a complexity of $O(m^{2})$. Hence, assuming that the simplicial complex $K$ and its dimension $d$ are fixed, algorithm 2 can be implemented with a complexity of $O(m^{2})$ in respect to the number of critical simplices $m$.

Algorithm 3 accepts a graph as input and determines the filtration weights of the simplices in the $d$-dimensional clique complex $K$ of $G$. Similar to the analysis of algorithm 2, we assume that the graph $G$ and the dimension $d$ of $K$ is fixed which implies that the number of simplices in the $d$-dimensional clique complex $K$ of $G$ is also constant. In such a situation, lines 1-7 of algorithm 3 have constant complexity which is independent of the number of critical simplices $m$. As seen in the preceding paragraph, provided $K$ and $d$ are fixed, the complexity of algorithm 2 and hence the complexity of line 8 of algorithm 3 is $O(m^{2})$. Lines 9-11, lines 13-24 and lines 26-31 in the pseudocode of algorithm 3 have constant complexity which is independent of the number of critical simplices $m$. The number of times the operations of lines 13-24 in algorithm 3 are repeated is equal to the number of critical weights which is less than or equal to the number of critical simplices $m$. Thus, assuming that the graph $G$ and the dimension $d$ of the corresponding clique complex $K$ are fixed, algorithm 3 can be implemented with a complexity of $O(m^{2})$ in respect to the number of critical simplices $m$.

Theoretical results on stability of persistent homology and persistence diagrams of discrete

Morse functions

Stability of persistent homology

In Algorithm 1, some choices are made in assigning values to every simplex in order to get a discrete Morse function. A natural question arises which is the following. Under what conditions do two discrete Morse functions give the same persistent homology groups. We remark that in Forman’s theory Forman2002 , the actual values of the discrete Morse functions are less important than the gradient vector field that the discrete Morse function induces. For the definition of a gradient vector field, we refer the reader to Definition 3.3 in Forman’s Forman2002 article. In his expository paper Forman2002 , Forman asserts on page 15 that: “In fact, this gradient vector field contains all of the information that we will need to know about the function for most applications.”

Note that two Morse functions $f$ and $g$ induce the same gradient vector field if and only if they satisfy the following condition (see Forman Forman2005 or Theorem 3.1 in Ayala et al.Ayala2009 ):

[TABLE]

We show that for two discrete Morse functions to have the same filtration given by their level subcomplexes on critical weights, the actual values of the discrete Morse function are not as important as the order relation induced by the function on the simplices. The discrete Morse functions with the same filtration clearly give the same persistent homology groups. In particular, we show below that a sufficient condition for this to happen is that both functions have the same order relation with respect to the values assigned to the simplices.

Lemma 7. Let $f$ and $g$ be two discrete Morse functions on a finite simplicial complex $K$. Assume that the following condition is satisfied:

[TABLE]

Then $f$ and $g$ have the same set of critical simplices $C:=\{c_{1},c_{2},\cdots,c_{k}\}$, and the level subcomplex of $f$ corresponding to the critical cell $c_{i}$ is the same as the level subcomplex of $g$ corresponding to $c_{i}$, for any $i=1,2,\cdots,k.$

Proof. The condition C1 guarantees that the functions $f$ and $g$ have the same gradient vector field, since C1 implies C0. Therefore they have the same set of critical simplices $C$. Denote the level subcomplex of $f$ (respectively, of $g$) corresponding to the critical cell $c_{i}$ by $K^{f}(c_{i})$ (respectively, by $K^{g}(c_{i})$). Now, if $\alpha\in K^{f}(c_{i})$, then there exists $\beta\in K$ such that $\alpha\leq\beta$ and $f(\beta)\leq f(c_{i})$. Note that $\alpha\leq\beta$ represents that either two simplices $\alpha$ and $\beta$ are same or $\alpha$ is a face of $\beta$. By assumption C1 this implies that $g(\beta)\leq g(c_{i})$, thus $\alpha\in K^{g}(c_{i})$. So we get $K^{f}(c_{i})\subseteq K^{g}(c_{i})$. Reversing the roles of $f$ and $g$ in the above argument, we also get $K^{g}(c_{i})\subseteq K^{f}(c_{i})$. Thus $K^{f}(c_{i})=K^{g}(c_{i})$ for any critical cell $c_{i}\in C$. ∎

A simple example of two discrete Morse functions $f$ and $g$ which satisfy C1 is the following. Starting with a discrete Morse function $f$ obtain $g$ by sliding the values of $f$ by a fixed real number $r$. A slightly more sophisticated example of two functions which satisfy C1 is the following. Given a discrete Morse function $f$ on a finite simplicial complex $K$, one can arrange the values of $f$ on the simplices of $K$ in non-decreasing order. Suppose that $\alpha$ and $\beta$ are simplices of $K$ which are adjacent in this ordering, i.e., we have $f(\alpha)\leq f(\beta)$ and no other simplex $\sigma$ exists such that $f(\sigma)$ lies in the interval $[f(\alpha),f(\beta)]$, and let us call such pairs of simplices $(\alpha,\beta)$ as $f$-successive. If for each $f$-successive pair of simplices $(\alpha,\beta)$, one adds a non-negative number $\delta(\alpha,\beta)$ to $f(\alpha)$ such that $\delta(\alpha,\beta)\leq f(\beta)-f(\alpha)$, we obtain a new discrete Morse function $g$ given by $g(\alpha):=f(\alpha)+\delta(\alpha,\beta)$, for which any pair $(\alpha,\beta)$ is $f$-successive if and only if it is $g$-successive. It is easy to deduce that $f$ and $g$ then satisfy the condition C1. Thus, $f$ and $g$ will have the same persistent homology groups.

Stability of persistence diagrams

We remark that two discrete Morse functions may have the same filtration (and thus the same persistent homology groups) but their persistence diagrams may differ if their critical weights (i.e., values on the critical simplices) are different. A sufficient condition that assures that two discrete Morse functions $f$ and $g$ have the same persistence diagrams is the following, keeping the notations of Lemma 7:

[TABLE]

Even if the condition C2 is not satisfied, there are stability theorems which give conditions under which the bottleneck distance between the respective persistence diagrams does not change much, i.e., the persistence diagrams are close in the bottleneck metric. The first such stability theorem was given by Cohen-Steiner et al. Cohen-Steiner2007 . This result Cohen-Steiner2007 is not directly applicable to discrete Morse functions, as they are not continuous. However, we show below that stability theorems for discrete Morse functions for finite regular CW-complexes can be inferred from stability results in Chazal et al. Chazal2008 (see also Bauer et al.Bauer2012 ). For simplicity we give the arguments for discrete Morse functions on a finite simplicial complex, although all the results below are also valid for finite regular CW-complexes.

Recall the $\infty$-Wasserstein distance or bottleneck distance between two multisets $X$ and $Y$ in $\mathbb{R}^{2}$:

[TABLE]

In the above equation, the supremum is taken over all bijections $\eta:X\rightarrow Y$ (with the convention that a point with multiplicity $k\in\mathbb{N}$ is considered as $k$ individual points) and for $(a,b)\in\mathbb{R}^{2}$, $||(a,b)||_{\infty}:=\max\{|a|,|b|\}$.

For a discrete Morse function $f$ with critical values $w_{i}:=f(c_{i}),i=1,2,\cdots,n$, the $k^{th}$ persistence diagram $D^{k}f$ (or simply $Df$ when the value of $k$ is fixed), is defined as follows. Consider the multiset of points $W^{k}_{f}:=\{(w_{i},w_{j}):w_{i}<w_{j},i,j=1,2,\cdots,n\}$ with each point $(w_{i},w_{j})$ endowed with the multiplicity $\mu_{k}(w_{i},w_{j})$ given by (see e.g. Di Fabio et al. DiFabio2015 ):

[TABLE]

where $\beta_{x}^{y}:=\operatorname{rank}(H_{k}(K^{f}(x))\rightarrow H_{k}(K^{f}(y)))$ for $x,y\in\mathbb{R}$ with $x<y$. Denote by $\Delta$ the diagonal in $\mathbb{R}^{2}$ considered as a multiset with infinite multiplicity given to each of its points.

Definition 8. The persistence diagram $D^{k}f$ is the subset of $W^{k}_{f}\cup\Delta$ consisting of points $(u,v)$ with $\mu_{k}(u,v)>0$.

We state below the stability results from Chazal et al. Chazal2008 that we shall need, and to do so, we recall the concepts of persistence modules and $\epsilon$-interleaving:

Definition 9. (Persistence modules (Definition 2.2. in Chazal et al. Chazal2008 )) Let $R$ be a commutative ring with unity. A collection ${\mathcal{F}}=\{F^{\alpha}\}_{\alpha\in\mathbb{R}}$ of $R$-modules $F^{\alpha}$ together with homomorphisms $\phi_{\alpha\beta}:F^{\alpha}\rightarrow F^{\beta}$ for $\alpha\leq\beta$, is called a persistence module if

[TABLE]

A persistence module is called $\delta$-tame if for any $\alpha<\alpha+\delta<\beta$ one has $\operatorname{rank}\phi_{\alpha\beta}<\infty$.

Definition 10. ($\epsilon$-interleaving of functions (Definition 4.1 in Chazal et al. Chazal2008 )) Let $\epsilon>0$, and let $f$ and $g$ be two real-valued functions on a finite simplicial complex $K$. For $x\in\mathbb{R}$, denote the sub-level sets $f^{-1}(-\infty,x]$ of $f$ (respectively, $g$) by $F^{x}$ (respectively, $G^{x}$). The functions $f$ and $g$ are said to be strongly $\epsilon$-interleaved if for any $a\in\mathbb{R}$ and $n\in\mathbb{Z}$, one has

[TABLE]

The main result for $\epsilon$-interleaved functions that we shall need is as follows:

Lemma 11. (Lemma 4.2(iii) in Chazal et al. Chazal2008 ) Given $\epsilon>0$, two functions $f$ and $g$ are strongly $\epsilon$-interleaved if and only if $||f-g||_{\infty}<\epsilon$.

Definition 12. ($\epsilon$-interleaving of persistence modules (Definition 4.3 in Chazal et al. Chazal2008 )) Let ${\mathcal{F}}=\{F^{\alpha}\}_{\alpha\in\mathbb{R}}$ and $\mathcal{G}=\{G^{\alpha}\}_{\alpha\in\mathbb{R}}$ be persistence modules over a commutative unital ring $R$. Let $\epsilon>0$, then ${\mathcal{F}}$ and $\mathcal{G}$ are said to be strongly $\epsilon$-interleaved if for all $a\in\mathbb{R}$ and $n\in\mathbb{Z}$, there exist homomorphisms $\Phi_{a+2n\epsilon}:F^{a+2n\epsilon}\rightarrow G^{a+(2n+1)\epsilon}$ and $\Psi_{a+(2n+1)\epsilon}:G^{a+(2n+1)\epsilon}\rightarrow F^{a+(2n+2)\epsilon}$ such that the following diagram commutes:

[TABLE]

If $f$ and $g$ are strongly $\epsilon$-interleaved, then the associated persistence modules given by their $k^{th}$ singular homology groups are also strongly $\epsilon$-interleaved for any $k\in\mathbb{N}\cup\{0\}$, the associated homomorphisms in Definition 10 are induced by the corresponding inclusion maps. Note that these persistence modules are [math]-tame since the functions are defined on a finite simplicial complex $K$.

One can define persistence diagrams for $\delta$-tame persistence modules (see Definitions 3.6 and 3.10 in Chazal et al. Chazal2008 ). Denote the persistence diagram of a $\delta$-tame persistence module ${\mathcal{F}}$ by $D_{\delta}{\mathcal{F}}$. The persistence diagram of persistence modules generalizes the classical persistence diagrams associated with real-valued functions, i.e., if $f$ is a tame real-valued function on a finite simplicial complex $K$, the associated [math]-tame persistence module given by the $k^{th}$ singular homology groups of its sublevel sets $F^{\alpha}$, denoted ${\mathcal{F}}_{k}:=\{H_{k}(F^{\alpha})\}_{\alpha\in\mathbb{R}}$, has the same persistence diagram as the $k^{th}$ persistence diagram of $f$, i.e. $D^{k}f=D_{0}{\mathcal{F}}_{k}$ (see Remark 1, page 12 in Chazal et al. Chazal2008 for a short proof).

The next lemma shows that if $f$ is a discrete Morse function and $K^{f}(\alpha),\alpha\in\mathbb{R}$ are the level subcomplexes of $f$, then the persistence diagram of the [math]-tame persistence module $\mathcal{K}_{k}=\{H_{k}(K^{f}(\alpha))\}_{\alpha\in\mathbb{R}}$ is the same as the $k^{th}$ persistence diagram $D^{k}f$ of $f$ corresponding to the filtration of $K$ given by the level subcomplexes at critical weights, i.e., $D^{k}f=D_{0}\mathcal{K}_{k}$. The proof is essentially the same as the arguments given in Remark 1, page 12 of Chazal et al. Chazal2008 , except that we use Forman’s Lemma 2.6 Forman2002 for comparing homology groups across critical weights. Note that due to Forman’s Lemma 2.6 Forman2002 , the set of $k$-homological critical values of $f$, i.e., points $c\in\mathbb{R}$ such that for any $\epsilon>0$ the homomorphism on the $k$-homology groups induced by inclusion map $K^{f}(c-\epsilon)\hookrightarrow K^{f}(c+\epsilon)$ is not an isomorphism, is a subset of its critical weights in general. However, it is easy to see that only those points $(w_{i},w_{j})$ appear in the persistence diagram of $f$ for which both $w_{i}$ and $w_{j}$ are $k$-homological critical values of $f$ (otherwise the multiplicity of such points is zero).

Lemma 13. We have $D^{k}f=D_{0}\mathcal{K}_{k}$.

Proof. Let $k$ be fixed, to simplify the notation we denote the $k$-persistence diagrams of $f$ and $\mathcal{K}_{k}$ by $Df$ and $D\mathcal{K}$, respectively. From the definition of the persistence diagrams of persistence modules (Definition 3.6 in Chazal et al. Chazal2008 ), it suffices to show that there exists $\epsilon>0$ $x\in(0,\epsilon)$, such that $W_{\infty}(Df,D\mathcal{K}_{\frac{\epsilon}{2^{n}},x})<\epsilon/2^{n}$ for any $n\geq 1$, where $D\mathcal{K}_{\frac{\epsilon}{2^{n}},x}$ is the multi-subset of $\mathbb{R}^{2}$ given by points of non-zero multiplicity in the union of the grid

[TABLE]

and the diagonal $\Delta$, where the multiplicity $\mu(m,m^{\prime})$ of a point $(x+\frac{m\epsilon}{2^{n}},x+\frac{m^{\prime}\epsilon}{2^{n}})\in G_{\epsilon/2^{n},x}$ is given by the following formula from Definition 3.1(ii) in Chazal et al. Chazal2008 :

[TABLE]

where we have used the notation $H(m):=H_{k}(K^{f}(x+\frac{m\epsilon}{2^{n}})),m\in\mathbb{Z}$ for simplicity.

Since the number of critical weights is finite, one can choose $\epsilon>0$, $x\in(0,\epsilon]$ such that none of the grids $G_{\epsilon/2^{n},x}$ have a point outside the diagonal that coincides with a point in $Df$. Then, for each $n\geq 1$ a point $(w_{i},w_{j})\in Df$ is contained in exactly one of the cells $C$ in the grid $G_{\epsilon/2^{n},x}$. It is easy to check using Forman’s Lemma 2.6 Forman2002 that the multiplicity $\mu(w_{i},w_{j})$ is then equal to the multiplicity of the upper-right corner of $C$, while all other off diagonal points in $D\mathcal{K}$ which are not upper-right corners of cells containing a point of $Df$, have multiplicity zero. Therefore the bijection $\Phi_{\epsilon,n}:Df\rightarrow D\mathcal{K}_{\epsilon/2^{n},x}$ which maps $(w_{i},w_{j})$ in $Df$ to the upper right corner of the unique cell in $G_{\epsilon/2^{n},x}$ containing $(w_{i},w_{j})$ and any point in the region $\{(x,y)\in\mathbb{R}^{2}:y\in[x,x+\epsilon/2^{n}]\}$ to the nearest point on the diagonal, moves points by a distance at most $\frac{\epsilon}{2^{n}}$. Thus, we get $W_{\infty}(Df,D\mathcal{K}_{\frac{\epsilon}{2^{n}},x})<\epsilon/2^{n}$ for any $n\geq 1$. ∎

We can now state the stability theorem for $\epsilon$-interleaved persistence modules:

Theorem 14. (Theorem 4.4 in Chazal et al. Chazal2008 ) Let ${\mathcal{F}}$ and $\mathcal{G}$ be two strongly $\epsilon$-interleaved persistence modules which are $\delta$-tame for some $\delta\geq 0$. Then,

[TABLE]

We are now ready to state a stability theorem for persistence diagrams for discrete Morse functions:

Theorem 15. Let $f$ and $g$ be discrete Morse functions on a finite simplicial complex $K$ which induce $k^{th}$ persistence diagrams $D^{k}f$ and $D^{k}g$, respectively. Suppose that there exists $\epsilon>0$, such that $||f-g||_{\infty}<\epsilon$. Then,

[TABLE]

Proof. Let us show that the level subcomplexes of $f$ and $g$ are strongly $\epsilon$-interleaved. Indeed, since $||f-g||_{\infty}<\epsilon$, $f$ and $g$ are strongly $\epsilon$-interleaved by Lemma 11. Let $c\in\mathbb{R}$ and let $\sigma\in K^{f}(c)$ and denote by $F^{c}$ the sub-level set $f^{-1}((-\infty,c])$ of $f$; $G^{c}$ is defined similarly for the discrete Morse function $g$. Then, by Lemma 3.2 of Forman Forman1998 , there exists a co-face $\tau$ of $\sigma$ such that $f(\tau)\leq c$, i.e. $\tau\in F^{c}$. By the strong $\epsilon$-interleaving property, one gets $F^{c}\subseteq G^{c+\epsilon}$, so we have $g(\tau)\leq c+\epsilon$, which implies that $\sigma\in K^{g}(c+\epsilon)$, so that we get $K^{f}(c)\subseteq K^{g}(c+\epsilon)$. Similarly one easily gets $K^{g}(c)\subseteq K^{f}(c+\epsilon)$ for any $c\in\mathbb{R}$. This implies that the persistence modules of $f$ and $g$ given by the $k^{th}$ homology groups of the level subcomplexes are strongly $\epsilon$-interleaved for any $k\in\mathbb{N}\cup\{0\}$. Hence by Theorem 14, one gets:

[TABLE]

This concludes the proof. ∎

As a concluding remark to this section, we observe that the bound of $3\epsilon$ can be improved to $\epsilon$, following the same arguments as given in the proof of Theorem 4.8 in Chazal et al. Chazal2008 , we omit the details here.

Illustration of the stability of persistence diagrams using a real network example

As an illustration of the stability property, we compute the bottleneck distance between persistence diagrams of two discrete Morse functions $f_{1}$ and $f_{2}$ on the Yeast protein interaction network Jeong2001 generated by making small variations in Algorithm 1. The construction of $f_{1}$ uses a uniform distribution on $[0.1,0.4)$ for generating the values of $\epsilon$ in Algorithm 1, while the construction of $f_{2}$ uses a normal distribution with mean $\mu=0.25$ and standard deviation $\sigma=0.05$ for generating the values of $\epsilon$ in Algorithm 1. The $L^{\infty}$-norm of the difference $||f_{1}-f_{2}||_{\infty}$ is $18.3782$, while the bottleneck distance between their total persistence diagrams (i.e., the persistence diagram of the persistence module $\mathcal{K}_{*}:=\{\oplus_{k=0}^{\dim K}H_{k}(K^{f}(\alpha))\}_{\alpha\in\mathbb{R}}$), is computed to be $W_{\infty}(Df_{1},Df_{2})=0.0061$. The barcode diagrams for the $0^{th}$ and $1^{st}$ homology groups in the Yeast protein interaction network Jeong2001 with the two discrete Morse functions, $f_{1}$ and $f_{2}$, are shown in Appendix figure A1.

**Supplementary Information (SI)

Figures**