Construction of simplicial complexes with prescribed degree-size sequences

Capturing higher-order dependencies is essential for the structural interpretation of the organization and behavior of complex systems [1, 2, 3]. Simplicial complex modeling, among other methods in applied topology [4, 5, 6], provides a combinatorial description of the topology of the system, featuring algebraic redundancies that may be compressed out using equivalence relations. Similar to networks, simplicial complexes are composed of nodes that represent system observables, and high-dimensional analogs of edges, called simplices, that represent polyadic relationships among the nodes. Simplicial modeling made several recent discoveries in complex systems, including the emergence of discontinuous transitions [7] in epidemic spreading [8, 9] and synchronization [10, 11], the multiscale hierarchy in adaptive voter models [12], the role of simplex size fluctuations in temporal networks [13], localization of dynamics [14] and percolation transitions [15, 16, 17]. In a related role, simplicial complexes have been used to summarize static features, addressing questions about the local geometry of data, such as in distinguishing the voting patterns in densely populated cities [18] and in describing the shape of scientific collaborations [19].

Dynamical processes on networks depend crucially on the network structure [20, 21]. However, their generalizations to simplicial structures and nonpairwise interactions are relatively less understood. Notably, the basic question of which degree-size sequences can be realized by a complex is still open. In this Letter we make progress in this direction by extensive numerical experiments.

Let $V$ be a finite set of nodes. A simplicial complex on $V$ is a collection $K$ of nonempty subsets of $V$, called simplices, subject to two requirements: First, for each node $v\in V$, the singleton $\{v\}\in K$; second, for all simplices $\tau\in K$, all its proper subsets $\sigma\subset\tau$ must also be in $K$. A facet is a maximal simplex, i.e., one that is not a subset of any other simplex. Note that a simplicial complex can be fully specified by listing only its facets, and we follow that convention here. We define the degree $d_{i}$ of a node $v_{i}$ as the number of facets incident on $v_{i}$ and the size (cardinality) $s_{j}$ of a facet $\sigma_{j}$ as the number of nodes it contains. For any simplicial complex $K$, there is a corresponding degree-size sequence $\textbf{{t}}_{K}=(\textbf{{d}},\textbf{{s}})$, where $\textbf{{d}}=\{d_{1},d_{2},\dots,d_{n}\}$ and $\textbf{{s}}=\{s_{1},s_{2},\dots,s_{m}\}$ (see Fig. 1). However, the reverse statement is not always true. There are sequences t that cannot be realized by any simplicial complex. Hence, inspired by Ref. [22], we pose the simplicial complex realization problem: Given integer sequences $\textbf{{t}}=(\textbf{{d}},\textbf{{s}})$, does there exist a simplicial complex $K$ with that degree-size sequence? When the answer is affirmative, we call the sequence simplicial.

The degree-size sequence reflects the local coupling patterns of the system. Models that constrain these features can often be used as null models that detect salient structures. In particular, the simplicial configuration model (SCM) [22] specifies a distribution of simplicial complexes with fixed degree-size sequences, and can be sampled via a Markov chain Monte Carlo (MCMC) algorithm. The SCM extends the configuration model for graphs [23, 24] and similar efforts in simplicial complexes of equal-size facets [25]. Critically, the MCMC requires an initial simplicial configuration to work, restricting its use to empirical data which can act as the initialization. The algorithm we propose yields an initialization for arbitrary simplicial sequences, not necessarily taken from an observed data set.

Related work.—A key difficulty in simplicial complex realization is that no facet is allowed to be completely included in any other facet; we call this the “no-inclusion constraint.” In contrast, hypergraphs have no such constraint. For simple hypergraphs, Ref. [26] gives a rejection sampling algorithm that samples realizations with given degrees and hyperedge sizes. For their non-simple counterparts, Ref. [27] considers an MCMC approach for generating such hypergraphs and Ref. [28] ensures the existence of a starting realization which, however, needs not be simplicial. If we relax the notion of facets and consider instead the degree sequence of nodes partaking in given motifs, Refs. [29, 30] use generating functions to study such networks.

Testing whether a degree-size sequence is simplicial is a generalization of the graph realization problem, in which the main result is the Erdős–Gallai theorem [31] that exactly characterizes graphical sequences with a set of easy-to-test inequalities. Equivalently, a particular graph realization can be constructed by a recursive application of the Havel–Hakimi theorem [32, *hakimiRealizabilitySetIntegers1962]. The reformulation of these theorems expedites the direct sampling of networks and facilitates understanding of their properties (e.g., all scale-free networks are sparse [34]). Moreover, networks enjoy tractable expressions for many ensemble-averaged quantities of interest (e.g., degree correlation, which tends to be disassortative in heterogeneous networks [35, 36]), precisely because they are free from the no-inclusion constraint. For simpliciality testing, however, none of these methods applies. Recently, Deza, Levin, Meesum, and Onn [37] proved that deciding whether a given sequence is the degree sequence of a simple $s$-uniform hypergraph is $\mathsf{NP}$-complete when $s\geq 3$, through a reduction from the 3-partition problem [38]. Interestingly, simple $s$-uniform hypergraphs are equivalent to equidimensional simplicial complexes, because when all hyperedges are the same size, they automatically satisfy the no-inclusion constraint. This implies that deciding simpliciality is hard and there may not exist an efficient algorithm to enumerate and sample these instances. Yet, as our exploration reveals, not all sequence combinations are equally hard. For example, for any s, taking the trivial degree sequence $\textbf{{d}}=\{1,\dots,1\}$ immediately yields a simplicial realization.

In this Letter, we develop a deterministic, backtracking-based search algorithm that always correctly solves simpliciality, and present evidence that on most instances it runs in polynomial time. We then explore the topology of the constructed complex more generally via its Betti numbers, as well as using it as a seed for the SCM ensemble. With the randomized realizations, we reveal the regimes in which their expected topology changes as a function of the degree and size sequences.

Algorithm.— Our algorithm proceeds as follows. It is given as input a node degree sequence d and a facet size sequence s, where both sequences are nonincreasing. Let $n=\lvert\textbf{{d}}\rvert$ and $m=\lvert\textbf{{s}}\rvert$, where $|\cdot|$ stands for the cardinality. Simpliciality fails trivially if there are fewer $1$’s in d than in s, as it is doomed to violate the no-inclusion constraint, or if $d_{1}>m$ or $s_{1}>n$, or $\sum{\textbf{{d}}}\neq\sum{\textbf{{s}}}$, as it would be impossible to form an incidence matrix. As a preprocessing step, we pair the $1$’s in d with those in s to make a partial output.

Next, we pick $s_{1}$ nodes to make a candidate facet $\hat{\sigma}^{1}$. This selection is not done stochastically, but favors higher-degree nodes—the candidate facet $\hat{\sigma}$ with the largest sum of input degrees $w(\hat{\sigma})=\sum_{i\in\hat{\sigma}}{d_{i}}$ is preferred. We ensure that $\hat{\sigma}$ is not included in any existing facet, or we proceed with the next one in heuristic order. We will validate $\hat{\sigma}^{1}$ with a series of rules. If it fails any of them, we take the next candidate, until we accept and advance to pick $s_{2}$ nodes for $\hat{\sigma}^{2}$. For $\hat{\sigma}^{\ell}$ ($\ell\geq 2$), the facets with larger $w(\hat{\sigma}^{\ell})$ still attain precedence.

With this overall structure in mind, we now express the algorithm recursively. At each branching stage $\ell$, the input is a 3-tuple $(\textbf{{d}}^{\ell},\textbf{{s}}^{\ell},\textbf{{B}}^{\ell})$, where $\textbf{{d}}^{\ell}$ is the residual degree sequence and $\textbf{{s}}^{\ell}$ the residual size sequence, denoting the marginal sums of the incidence matrix that still need to be fulfilled. In addition, we have a collection of “blocking” facets $\textbf{{B}}^{\ell}\coloneqq\{\sigma^{k}:1\leq k<\ell\}$, where each accepted facet $\sigma^{k}$ plays a role in the no-inclusion constraint. After accepting $\hat{\sigma}^{\ell}$, the output is again a 3-tuple $(\textbf{{d}}^{\ell+1},\textbf{{s}}^{\ell+1},\textbf{{B}}^{\ell+1})$. The algorithm returns a simplicial realization if $\sum{\textbf{{s}}^{\ell}}=0$, or a negative result when the entire search tree has been traversed.

To validate a candidate facet $\hat{\sigma}^{\ell}$, we assume that $\hat{\sigma}^{\ell}$ is accepted and virtually move to the next branching stage to obtain a number of intermediate data, which must obey the validation rules before we actually branch. Precisely, we compute the 3-tuple $(\textbf{{d}}^{\ell+1},\textbf{{s}}^{\ell+1},\textbf{{B}}^{\ell+1})$ at stage $\ell+1$ and the set of non-shielding nodes $Q^{\ell}\coloneqq V^{\ell}\setminus\hat{\sigma}^{\ell}$, where $V^{\ell}$ is the set of nodes at stage $\ell$. For clarity, we drop all superscripts in the following developments.

Rule 1 (for d and s) ¹¹1This rule is similar to, but stronger than, the Gale–Ryser characterization of bigraphic pairs of sequences [63, *ryser_1957], which would allow $s_{\text{max}}=\lvert\textbf{{d}}\rvert$ for non-terminal facets, thus violating the no-inclusion constraint.. We ask that $d_{\text{max}}\leq\lvert\textbf{{s}}\rvert$, and that $s_{\text{max}}\leq\lvert\textbf{{d}}\rvert$, where equality in the latter requirement is only allowed when $\lvert\textbf{{s}}\rvert=1$.

Rule 2 (for d, s, and $Q$). We require that $\lvert\textbf{{s}}\rvert\leavevmode\nobreak\ \leq\leavevmode\nobreak\ \sum_{i\in Q}{d_{i}}$, as every subsequent facet must contain at least one non-shielding node (i.e., element of $Q$) in order to not be included in $\hat{\sigma}$.

If equality holds, each subsequent facet must take exactly one non-shielding node. To secure its availability, we can thus further require that $s_{\text{max}}-1\leq\lvert\textbf{{d}}\rvert-\lvert Q\rvert$.

Rule 3 (for d and B). Let $V^{\prime}$ be the set of nodes with positive residual degree. We require that, for every blocking facet $\sigma\in\textbf{{B}}$, $V^{\prime}$ be not a subset of $\sigma$.

While Rule 1 is essentially the precondition, Rules 2–3 are meant to cut the combinatorial explosion as much as possible, but not prohibit any realizable sequence from being realized. If accepting a candidate facet leads to $d_{i}=\lvert\textbf{{s}}\rvert$ for any node $i$, these nodes are required to be used for the remaining facets—we invoke a subroutine to recursively consume those “forced” degrees. A Python implementation of the algorithm is freely available [40]. The pseudocode is given in the Supplemental Material [41].

Easy & hard instances.—To benchmark the algorithm, we define the running time as

$$\tau=\tau_{\text{b}}+\tau_{\text{r}}\ ,$$

where $\tau_{\text{b}}$ records the number of times the algorithm backtracks (due to candidate depletion) and $\tau_{\text{r}}$ records the number of times a proposed facet is rejected. These numbers are correlated: If we come up with better validation rules, we reduce rejections and prevent backtracking. We call a degree-size sequence easy (or polynomial) if $\tau=0$, meaning that no backtracking is necessary, the algorithm either finds a realization in near-linear time or rejects simpliciality immediately. Otherwise, we call it hard. Of course, this distinction is dependent on the choice of algorithm and “hard” instances with $\tau$ small will still be easy in practice, but for the purposes of understanding the problem, we find it useful to distinguish between those cases that are solved greedily by our algorithm and those that are not.

In Fig. 2(a), we show the realizability of all combinations of degree-size sequences of a fixed instance size $E=\sum{\textbf{{d}}}=13$. The self-similar pattern reflects the recursive nature of the algorithm. However, we are unable to conclude a recurrence relation in the spirit of Erdős–Gallai, due to the inherent complexity of solving particular instances. Indeed, there are two distinct populations of easy and hard instances, where a major fraction of inputs falls in the polynomial region. The easy majority can be understood as the iterative application of the heuristic policy. For 3-uniform, $\mathsf{NP}$-hard instances [Fig. 2(b)], most inputs are polynomial and, on the average, the non-simplicial instances are harder than simplicial ones, as tree exhaustion is needed. This highlights a useful property that false negatives are kept at minimum when applying a reasonable cutoff.

To investigate the dependence of the simpliciality and hardness of degree-size sequences on the instance size, we perform extensive numerical calculations, generating uniform ensembles of sequences of random integers, $(\textbf{{d}},\textbf{{s}})$, with a range up to $E=1010$. We test each sequence for simpliciality by applying the algorithm and compute for each $E$ the simplicial fraction $S/N$, where $S$ is the total number of simplicial sequences in the ensemble and $N$ is the ensemble size, chosen at $N=10^{6}$. We also compute the fractions $P/N$ and $S_{p}/S$, where $P$ is the total number of polynomial instances in the ensemble and $S_{p}$ is the number of instances that are both simplicial and polynomial. The results, plotted in Fig. 2(c), clearly demonstrate the persistent existence of easy and hard populations at much larger sizes. An open question is to what extent we can further separate these classes in polynomial time, perhaps through better validation rules.

Heuristic policy and topology—When data are encoded as a simplicial complex $K$, we can characterize their homotopical information by the Betti numbers $\beta_{k}(K)$ [44, 5]. They quantify the $k$-dimensional connectivity of objects by comparing their volume and boundary at each dimension—$\beta_{0}$ is the number of connected components, $\beta_{1}$ the number of homological cycles (or loops), while higher $\beta_{k}$ effectively counts the number of $k$-dimensional cavities. The topological cavities can be geometric in physical space, such as in granular packings [45], or abstract structures in experimental measurements [18, 19]. For instance, a cavity in diffusion MRI readings could indicate axonal dropout, a neurological disorder [46]. In the following, we focus on the first two Betti numbers because they are easier to interpret.

An important feature of our heuristic policy is that high-degree nodes tend to form larger facets, resulting in a core-periphery structure [47] with dangling and isolated facets on the fringe. Therefore, the heuristic tends to find a realization with a relatively large number of connected components and few loops ²²2The latter case is due to an algorithmic choice to favor the facets with small residual degrees, from among the candidates with the same priority $w$ [41].. We show in Fig. 3(a) this feature on an empirical degree-size sequence extracted from the human diseasome network [49, 22]. Indeed, our algorithm can discover a realization with $\beta_{0}=643$ and $\beta_{1}=5$ compatible with the degree-size sequence, whereas typical realizations sampled from the SCM have much lower $\beta_{0}$ and much higher $\beta_{1}$, see Fig. 3(b). This algorithmic trait is consistent across other datasets we examined [40]. More broadly, this priority policy shows a minimal example where adding degree-size correlations can introduce atypical Betti numbers. This sheds light on growth mechanisms that generate structures with specific homology, such as being totally connected [50] or containing many cycles [51].

SCM ensembles.—To test the method, we generate random degree-size sequences and test them for simpliciality. We generate from two Poisson distributions, modified so that all facets of size $1$ are guaranteed to be matched with a degree-$1$ node. For each simplicial sequence, we use the constructed complex to initialize an MCMC sampler for the SCM ensemble [22] and compute its mean Betti numbers $\bar{\pmb{\beta_{0}}}$ and $\bar{\pmb{\beta_{1}}}$. Then we take the average over the random partitions. Note that finding an initialization is inefficient with a rejection sampler. The result, shown in Fig. 4, is a systematic study of $\left<\bar{\pmb{\beta_{0}}}\right>$ and $\left<\bar{\pmb{\beta_{1}}}\right>$ with respect to $\bar{\textbf{{d}}}$ and $\bar{\textbf{{s}}}$, where $\bar{\textbf{{s}}}$ is the mean facet size and $\bar{\textbf{{d}}}$ is the mean node degree, excluding the nodes that are created to match the facets.

In Fig. 4(a), we observe that the average number of connected components $\left<\bar{\pmb{\beta_{0}}}\right>$, in the SCM ensemble, decreases with increasing $\bar{\textbf{{s}}}$, likely due to the reduction of facets with cardinality $1$ in the size sequence. However, the distribution of cycles $\left<\bar{\pmb{\beta_{1}}}\right>$ is considerably more complicated. In the low $\bar{\textbf{{s}}}$ regime, the simplicial complex acts as a dyadic network, where denser networks contain more loops. By contrast, in the high $\bar{\textbf{{s}}}$ regime, the system is abundant in large simplices. Once there is a fraction of higher-degree nodes, we have no other option but to bundle the large facets with those nodes, resulting in a tree-like, few loops complex. We supply a parallel analysis on $d$-regular degree distributions in Fig. 4(b), which tend to entail more loops than Poissonian ones with the same average degree, as they possess fewer high-degree nodes.

The decay of the average number of cycles $\left<\bar{\pmb{\beta_{1}}}\right>$ when the average degree $\bar{\textbf{{d}}}$ is increased is reminiscent of the law of large numbers for Betti numbers in random simplicial complexes (e.g., the Linial–Meshulam model [52] or the random clique complex [53]). We note that in these studies, the limiting behavior of Betti numbers is discussed in the context of increasing facet density, where the decay is driven by filling the $k$-dimensional cycles with $(k+1)$-simplices. Here, the system has a fixed number of simplices and the decay is driven by the no-inclusion constraints that prevent the realization of any such complexes.

We also observe that $\left<\bar{\pmb{\beta_{1}}}\right>$ is unimodal with respect to mean facet size $\bar{\textbf{{s}}}$ [Figs. 4(a) and 4(b)]. The unimodality comes from the competitive relationship between the gradually multifaceted local geometry and the depletion of available facets. When $\bar{\textbf{{s}}}$ is small, there are fewer large facets to avoid for smaller ones and, thus, increasing the facet sizes has a similar effect as densifying the degrees, which creates more loops [cf. Fig. 4(c), \Circled1 to \Circled2]. That said, the loops are destroyed as facets merge into larger ones [cf. Fig. 4(c), \Circled2 to \Circled3]. This suggests the existence of an optimal facet size distribution for loopy configurations, as the number of facets is limited.

Discussion.—In computational complexity, many graph problems are $\mathsf{NP}$-hard in general, but may be in $\mathsf{P}$ for certain classes of graphs [54]. Simplicial complex realization is yet another addition to the list. We present a depth-bounded branching algorithm whose complexity presents a rich structure. In particular, simpliciality can be solved in time $f(m)E^{c}$ for some constant $c\approx 1$, where $E^{c}$ is the time spent in validations and the prefactor counts the number of nodes in the execution tree. For easy instances, $f(m)=m$, and the algorithm is linear in instance size. Otherwise, we have $f(m)=b^{m+1}$ in general, where $b$ is the branching ratio that grows with instance size. We find that $b$ is highly heterogeneous as a function of the branching stage—the searcher stalls at mid-stages, not at the beginning or the end. It remains open to accurately parametrize the complexity of hard instances of the simpliciality problem, and to prove rigorously, if the algorithm runs in linear time on average. Finally, we note that the branching design has opened up an avenue to systematically improve the algorithm, for example, through stronger validation rules to reduce the branching ratio, or introducing variants by non-chronological backjumping or clause learning techniques, as critically used in modern Boolean satisfiability (SAT) solvers [55].

Aside from these computational complexity questions, the boundary between easy and hard instances deserves further attention. We find that the instances tend to be harder when $(\textbf{{d}}$, $\textbf{{s}})$ contain numerous uniform entries, whereas a Poisson-Poisson combination yields very little backtracking. The understanding of when and why their hardness differs is poor compared to what is known for constraint satisfaction problems [56] or the inference of stochastic block models [57]. This raises a number of open questions, for example, is there an algorithmic phase transition that separates easy and hard regions? Here, hard instances could mean either $\tau>0$ or really hard in some other sense, e.g., takes exponential time, as seen in SAT. Or, would there even be two different phase transitions? It is also known that graph isomorphism can be solved in linear time for random graphs [58, 59], by leveraging the fact that in random graphs the degree distribution is essentially never uniform, so that high-degree nodes help break symmetries. For simpliciality, it could be useful to investigate the dependency of algorithmic hardness on the degree sequence among equidimensional sequences. We believe that the solutions to these problems will require new insights in the statistical physics of computation [60, 61], notably, to identify the “order parameter” that characterizes the algorithmic barrier to hard instances [62].

In closing, we have developed a constructive algorithm to allow faster realization of simplicial complexes with arbitrary degree-size sequences. Our algorithm paves the way for a more comprehensive analysis of higher-order phenomena in terms of local structural attributes, revealing their roles in various dynamical systems.

Acknowledgements.—I am particularly grateful to Joshua A. Grochow for significant discussions. I am grateful to Alice Patania and Jean-Gabriel Young for helpful correspondence, as well as to Daniel B. Larremore for support. I acknowledge the BioFrontiers Computing Core at the University of Colorado Boulder for computing facilities. This work was funded in part by Grochow startup funds.

\foreach\x

in 1,…,0 See pages \x, of YEN_PRE_sm.pdf