Compressing network populations with modal networks reveals structural diversity

For our clustering method, we rely on the minimum description length (MDL) principle: the best model among a set of candidate models is the one that permits the greatest compression—or shortest description—of a dataset rissanen1978 . The MDL principle provides a principled criterion for statistical model selection and has consequently been employed in various applications ranging from regression to time series analysis to clustering hansen2001mdl . A large body of research uses the MDL principle for clustering data, including studies on MDL-based methods for mixture models that accommodate continuous georgieva2011cluster ; tabor2014cross and categorical data li2004entropy , as well as methods that are based on more general probabilistic generative models narasimhan2005q . The MDL approach has also been applied to complex network data, most notably for community detection algorithms to cluster nodes within a network Rosvall07 ; Peixoto14a ; kirkley2022spatial and for decomposing graphs into subgraphs koutra2014vog ; wegner2014subgraph ; bloem2020large ; young2021hypergraph ; bouritsas2021partition , but also for clustering entire partitions of networks Kirkley22Reps . Our methods are similar in spirit to the one presented in ref. Kirkley22Reps for identifying representative community divisions among a set of plausible network partitions. Both approaches involve transmitting first a set of representatives and then the dataset itself by describing how each partition or network differs from its corresponding representative. However, the methods differ substantially in their details since they address fundamentally different questions.

We consider an experiment in which the initial data are a population of networks consisting of $S$ undirected, unweighted networks $\mathcal{D}=\{\bm{D}^{(1)},...,\bm{D}^{(S)}\}$ on a set of $N$ labeled nodes. The networks record, for instance, the co-location patterns among kids in a class of $N$ students over $S$ class periods.

We aim to summarize these data with $K$ modal networks $\mathcal{A}=\{\bm{A}^{(1)},...,\bm{A}^{(K)}\}$ (also undirected and unweighted) on the same set of nodes, with associated clusters of networks $\mathcal{C}=\{C_{1},...,C_{K}\}$, where $C_{k}$ comprises networks that are similar to the mode $\bm{A}^{(k)}$. This summary would allow researchers to, for instance, perform all downstream network analyses on a small set of representative networks—the modes—instead of a large set of networks likely to include measurement errors and from which it is difficult to draw valid conclusions.

We assume for simplicity of presentation that all networks $\mathcal{D}$ and $\mathcal{A}$ have no self- or multi-edges, although we can account for them straightforwardly. While $K$ can be fixed if desired, we assume that it is unknown and must be determined from regularities in the data.

To select among all the possible modes and assignments of networks to clusters, we employ information theory and construct an objective function that quantifies how much information is needed to communicate the structure of the network population $\mathcal{D}$ to a receiver. Clustering networks in groups of mostly similar instances allows us to communicate the population $\mathcal{D}$ efficiently in three steps: first the modes, then the clusters, and finally the networks $\mathcal{D}$ themselves as a series of small modifications to the modes $\mathcal{A}$. The MDL principle tells us that any compression achieved in this way reveals modes and clusters that are genuine regularities of the population rather than noise rissanen1978 .

We first establish a baseline for the code length: the number of bits needed to communicate $\mathcal{D}$ without using any regularities. One way to do this is to first communicate the parameters of the population at a negligible information cost (size $S$, number of nodes $N$, and the total number of edges $E$ in all networks of $\mathcal{D}$) and then transmit the population $\mathcal{D}$ directly. There are ${N\choose 2}$ possible edge positions in each of the $S$ undirected networks in $\mathcal{D}$, or $S{N\choose 2}$ possible edge positions for the whole population. So these networks can be configured in ${S{N\choose 2}\choose E}$ ways. It thus takes approximately

$\displaystyle\mathcal{L}_{0}(\mathcal{D})=\log{S{N\choose 2}\choose E}$

(1)

bits to transmit these networks to a receiver. (We use the convention $\log\equiv\log_{2}$ for brevity.) Applying Stirling’s approximation $\log x!\approx x\log x-x/\ln(2)$, we obtain

$\displaystyle\mathcal{L}_{0}(\mathcal{D})\approx S{N\choose 2}H_{b}\left(\frac{E}{S{N\choose 2}}\right)$

(2)

written in terms of the binary Shannon entropy

$\displaystyle H_{b}(p)=-p\log p-(1-p)\log(1-p).$

(3)

In practice, we expect to need many fewer bits than $\mathcal{L}_{0}$ to communicate $\mathcal{D}$, because the population of networks will often have regularities. We propose a multi-part encoding that identifies such regularities by grouping similar networks in clusters $\mathcal{C}$ with modes $\mathcal{A}$, which proceeds as follows. First, we send a small number of modes $\mathcal{A}$ in their entirety, which ideally captures most of the heterogeneity in the population $\mathcal{D}$. This step is costly but will save us information later. We then send the network clusters $\mathcal{C}$ by transmitting the cluster label of each network $s\in\mathcal{D}$. Finally, we transmit the edges of networks in each cluster, using the already transmitted modes as a starting point to compress this part of the encoding significantly. The expected code length can be quantified using simple combinatorial expressions, and the configuration of modes $\mathcal{A}$ and clusters $\mathcal{C}$ that minimizes the total expected code length—the MDL configuration—provides a succinct summary of the data $\mathcal{D}$. Figure 1 summarizes the transmission process and the individual description length contributions.

The expected length of this multi-part encoding is the sum of the description length of each part of the code that has significant communication costs. The modes are the first objects that incur such costs. Following the same reasoning as before, we denote the number of edges in mode $k$ as $M_{k}$ and conclude that we can transmit the positions of the occupied edges in mode $\bm{A}^{(k)}$ using approximately

$\displaystyle\mathcal{L}(\bm{A}^{(k)})=\log{{N\choose 2}\choose M_{k}}\approx{N\choose 2}H_{b}\left(\frac{M_{k}}{{N\choose 2}}\right)$

(4)

bits, where the second expression results from a Stirling approximation as in Eq. (2). We can therefore transmit all the modes with a total code length of

$\displaystyle\mathcal{L}(\mathcal{A})=\sum_{k=1}^{K}\mathcal{L}(\bm{A}^{(k)})$

(5)

bits.

The next step is to transmit the cluster label $k$ of each network in $\mathcal{D}$. For this part of the code, we first send the number of networks $S_{k}$ in each cluster $k=1,...,K$ at a negligible cost and then specify a particular clustering compatible with these constraints. The multinomial coefficient ${S\choose S_{1}\ S_{2}\ \cdots\ S_{k}}$ gives the total number of possible combinations of these cluster labels. The information content of this step is thus

$\displaystyle\mathcal{L}(\mathcal{C})=\log{S\choose S_{1}\ S_{2}\ \cdots\ S_{k}}\approx S\ H\Big{(}\{S_{k}/S\}\Big{)},$

(6)

where we again use the Stirling approximation and where

$\displaystyle H\Big{(}\{q_{k}\}\Big{)}=-\sum_{k=1}^{K}q_{k}\log q_{k}$

(7)

is the Shannon entropy of a distribution $\{q_{k}\}$.

Finally, we transmit the network population $\mathcal{D}$ by sending the differences between the networks in each cluster and their associated mode. To calculate the length of this part of the code, we focus on a particular cluster $C_{k}$ and count the number of times we will have to remove an edge from the mode $\bm{A}^{(k)}$ when specifying the structure of networks in its cluster using $\bm{A}^{(k)}$ as a reference. We call these edges false negatives, and count them as

$\displaystyle n_{k}=\sum_{s\in C_{k}}\left|\bm{A}^{(k)}\setminus\bm{D}^{(s)}\right|,$

(8)

where we interpret $\bm{D}^{(s)}$ and $\bm{A}^{(k)}$ as sets of edges, so the summand is the number of edges in mode $k$ that are not in the network $s$. Similarly, we also require the number of edges that occur in the networks of cluster $k$ but not in the mode—the number of false positives,

$\displaystyle p_{k}=\sum_{s\in C_{k}}\left|\bm{D}^{(s)}\setminus\bm{A}^{(k)}\right|.$

(9)

Like the cluster sizes $S_{k}$ and edge counts per cluster $M_{k}$, the pairs $(n_{k},p_{k})$ can be communicated to the receiver at a comparatively negligible cost, and we ignore them in our calculations.

To estimate the information content of this part of the transmission, we count the number of configurations of false-negative and false-positive edges in $C_{k}$. Focusing first on the false negatives—the edges that must be deleted—we count that of the $S_{k}M_{k}$ edges in the $S_{k}$ copies of the mode of cluster $k$, $n_{k}$ will be false-negative edges that can be configured in ${S_{k}M_{k}\choose n_{k}}$ ways. Similarly, using the shorthand $M_{k}^{\ast}={N\choose 2}-M_{k}$ to denote the unoccupied pairs of nodes in the mode $k$, there are $S_{k}M_{k}^{\ast}$ locations in which we must place $p_{k}$ false-positive edges, for a total of ${S_{k}M_{k}^{\ast}\choose p_{k}}$ possible configurations of false-positive edges. The total information content required for transmitting the locations of the false-negative and false-positive edges of every network in cluster $k$ is thus

	$\displaystyle\ell_{k}$	$\displaystyle:=\mathcal{L}(\{\mathcal{D}^{(s)}|s\in C_{k}\}|\bm{A}^{(k)})$
		$\displaystyle=\log{S_{k}M_{k}\choose n_{k}}+\log{S_{k}M_{k}^{\ast}\choose p_{k}},$		(10)

which we approximate as

$$\ell_{k}\approx S_{k}M_{k}H_{b}\left(\frac{n_{k}}{S_{k}M_{k}}\right)+S_{k}M_{k}^{\ast}H_{b}\left(\frac{p_{k}}{S_{k}M_{k}^{\ast}}\right).$$

(11)

Summing over all clusters,

$\displaystyle\mathcal{L}(\mathcal{D}|\mathcal{A},\mathcal{C})=\sum_{k=1}^{K}\ell_{k},$

(12)

we obtain the total information content of the final step in the transmission process.

We obtain the total description length $\mathcal{L}(\mathcal{D})$ by adding the contributions of Eqs. (5), (6), and (12), as

$\displaystyle\mathcal{L}(\mathcal{D})=\mathcal{L}(\mathcal{A})+\mathcal{L}(\mathcal{C})+\mathcal{L}(\mathcal{D}|\mathcal{A},\mathcal{C}).$

(13)

This objective function allows for efficient optimization because we can express it as a sum of the cluster-level description lengths

$$\mathcal{L}_{k}(\bm{A}^{(k)},C_{k})=\mathcal{L}(\mathcal{A}^{(k)})+S\log\left(\frac{S}{S_{k}}\right)+\ell_{k},$$

(14)

giving

$\displaystyle\mathcal{L}(\mathcal{D})=\sum_{k=1}^{K}\mathcal{L}_{k}(\bm{A}^{(k)},C_{k}).$

(15)

Equations (4) and (10) provide explicit expressions for $\mathcal{L}(\mathcal{A}^{(k)})$ and $\ell_{k}$.

Equation (15) gives the total description length of the data $\mathcal{D}$ under our multi-part transmission scheme. By minimizing this objective function we identify the best configurations of modes $\mathcal{A}$ and clusters $\mathcal{C}$. A good configuration $\{\mathcal{A},\mathcal{C}\}$ will allow us to transmit a large portion of the information in $\mathcal{D}$ through the modes alone. If we use too many modes, the description length will increase as these are costly to communicate in full. And if we use too few, the description length will also increase because we will have to send lengthy messages describing how mismatched networks and modes differ. Hence, through the principle of parsimony, Eq. (15) favors descriptions with the number of clusters $K$ as small as possible but not smaller.

This framework can be modified to accommodate populations of bipartite or directed networks. For the bipartite case, we make the transformations ${N\choose 2}\to N_{1}N_{2}$ and $M_{k}^{\ast}\to N_{1}N_{2}-M_{k}$, where $N_{1}$ and $N_{2}$ are the numbers of nodes in each of the two groups. This modification reduces the number of available positions for potential edges. Similarly, for the directed case, we can make the transformations ${N\choose 2}\to N(N-1)$ and $M_{k}^{\ast}\to N(N-1)-M_{k}$, which increases the number of available edge positions.

Since Eq. (15) has large support, is not convex, and has many local optima, a stochastic optimization method is a natural choice for finding reasonable solutions rapidly. We exploit the objective function’s decoupling into a sum over clusters $k$ and implement an efficient merge-split Monte Carlo method for the search peixoto2020merge ; Kirkley22Reps . The method greedily optimizes $\mathcal{L}(\mathcal{D})$ using moves that involve merging and splitting clusters of networks $D^{(s)}\in\mathcal{D}$.

Our merge-split algorithm minimizes the description length in Eq. (15) by performing one of the following moves selected uniformly at random and accepting the move as long as it results in a reduction of the description length (15):

1. 1.

   Reassignment: Pick a network $s$ at random and move it from its current cluster $C_{k}$ to the cluster $C_{k^{\prime}}$ that results in the greatest decrease in the description length. Compute the modes $\bm{A}^{(k)}$ and $\bm{A}^{(k^{\prime})}$ that minimize the cluster-level description lengths $\mathcal{L}_{k}(\bm{A}^{(k)},C_{k})$ and $\mathcal{L}_{k^{\prime}}(\bm{A}^{(k^{\prime})},C_{k^{\prime}})$ using Eq. (14) and the procedure described below, conditioned on the networks in $C_{k}$ and $C_{k^{\prime}}$.

2. 2.

   Merge: Pick two clusters $C_{k^{\prime}}$ and $C_{k^{\prime\prime}}$ at random and merge them into a single cluster $C_{k}$. Compute the mode $\bm{A}^{(k)}$ that minimizes the cluster-level description length $\mathcal{L}_{k}(\bm{A}^{(k)},C_{k})$ using Eq. (14) and the procedure described below, conditioned on the networks in $C_{k}$. Finally, compute the change in the description length that results from this merge.

3. 3.

   Split: Pick a cluster $C_{k}$ at random and split it into two clusters $C_{k^{\prime}}$ and $C_{k^{\prime\prime}}$ using the following $2$-means algorithm. First assign every network in $C_{k}$ to the cluster $C_{k^{\prime}}$ or $C_{k^{\prime\prime}}$ at random. Refine the assignments by successively moving every network to the cluster $C_{k^{\prime}}$ or $C_{k^{\prime\prime}}$ that results in a greater decrease in the description length and compute the modes $\bm{A}^{(k^{\prime})}$ and $\bm{A}^{(k^{\prime\prime})}$ that minimize the cluster-level description lengths $\mathcal{L}_{k^{\prime}}(\bm{A}^{(k^{\prime})},C_{k^{\prime}})$ and $\mathcal{L}_{k^{\prime\prime}}(\bm{A}^{(k^{\prime\prime})},C_{k^{\prime\prime}})$, conditioned on the networks now in $C_{k^{\prime}}$ and $C_{k^{\prime}}$. After convergence of the $2$-means style algorithm, compute the change in the description length that results from this split of cluster $C_{k}$.

4. 4.

   Merge-split: Pick two clusters at random, merge them as in move 2, then perform move 3 on this merged cluster. These two moves in direct succession help reassign multiple networks simultaneously; their addition to the move set improves the algorithm’s performance.

Since these moves modify only one or two clusters, the change in the global description length $\mathcal{L}(\mathcal{D})$ can be recomputed quickly as updates to the cluster-level description lengths in Eq. (14). Every time a mode is needed for these calculations, we use the mode that minimizes the cluster-level description length $\mathcal{L}_{k}(\bm{A}^{(k)},C_{k})$ in Eq. (14). To find this optimal mode efficiently, we start with the “complete” mode

$\displaystyle\bm{A}_{\mathrm{comp}}^{(k)}=\bigcup_{s\in C_{k}}\bm{D}^{(s)},$

(16)

with an edge between nodes $i$ and $j$ if at least one network in the cluster contains the edge. We then greedily remove edges from $\bm{A}_{\mathrm{comp}}^{(k)}$ in increasing order of occurrence in the networks of $C_{k}$—starting first with edges only found in a single network and going up from there—and update the cluster-level description length as we go. After removing all edges from $\bm{A}_{\mathrm{comp}}^{(k)}$, the mode giving the lowest cluster-level description length is chosen as the mode for the cluster. This approach is locally optimal under a few assumptions about the sparsity of the networks and the composition of edges in the clusters (see Appendix A for details).

We run the algorithm by starting with $K_{0}$ initial clusters (this choice has a negligible effect on the results, see Appendix B) and stop when a specified number of consecutive moves all result in rejections, indicating that the algorithm has likely converged. The worst-case complexity of this algorithm is roughly $O(NS)$ (the worst case is a split move right at the start, see Appendix B). Appendix A details the entire algorithm, and Appendix B provides additional tests of the algorithm, such as its robustness for different choices of $K_{0}$.

To diagnose the quality of a solution, we compute the inverse compression ratio

$\displaystyle\eta(\mathcal{D})=\mathcal{L}_{\mathrm{MDL}}(\mathcal{D})/\mathcal{L}_{0}(\mathcal{D}),$

(17)

where $\mathcal{L}_{\mathrm{MDL}}(\mathcal{D})$ is the minimum value of $\mathcal{L}(\mathcal{D})$ over all configurations of $\mathcal{A},\mathcal{C}$, given by the algorithm after termination, and $\mathcal{L}_{0}$ is given in Eq. (2). Equation (17) tells us how much better we can compress the network population $\mathcal{D}$ by using our multi-step encoding than by using the naïve fixed-length code to transmit all networks individually. If $\eta(\mathcal{D})<1$, our model compresses the data $\mathcal{D}$, and if $\eta(\mathcal{D})>1$, it does not because we waste too much information in the initial transmission steps.

In Sec. II.2, we described a merge-split Monte Carlo algorithm to identify the clusters $\mathcal{C}$ and modes $\mathcal{A}$ that minimize the description length in Eq. (15). This algorithm samples the space of unconstrained partitions $\mathcal{C}$ of the network population $\mathcal{D}$. However, in many applications, particularly in longitudinal studies, we may only be interested in constructing contiguous clusters, where each cluster is now a set of networks where adjacent indexes $s\in\{1,...,S\}$ indicate contiguity of some form (temporal, spatial, or otherwise). Such constraints reduce the space of possible clusterings $\mathcal{C}$ drastically, and we can minimize the description length exactly (up to the greedy heuristic for the mode construction) using a dynamical program jackson2005algorithm ; bellman2013dynamic ; PAY2022 .

Before we introduce an optimization for this problem, we require a small modification to Eq. (14) for the cluster-level description length to accurately reflect the constrained space of ordered partitions $\mathcal{C}$ that we are considering. In our derivation of the description length, we assumed that the receiver knows the sizes $\{S_{k}\}$ of the clusters in $\mathcal{C}$. If we transmit these sizes in the order of the clusters they describe, the receiver will also know the exact clusters $\mathcal{C}$, since knowing the sizes $\{S_{k}\}$ is equivalent to knowing the cluster boundaries in this contiguous case. We can therefore ignore the term $S\log(S/S_{k})$ in Eq. (14) that tells us how much information is required to transmit the exact cluster configuration. This modification results in a new, shorter description length

$\displaystyle\mathcal{L}^{(\mathrm{cont})}_{k}(\bm{A}^{(k)},C_{k})=\mathcal{L}(\mathcal{A}^{(k)})+\ell_{k}$

(18)

and a new global objective

$\displaystyle\mathcal{L}_{\mathrm{cont}}(\mathcal{D})=\sum_{k}\mathcal{L}^{(\mathrm{cont})}_{k}(\bm{A}^{(k)},C_{k}).$

(19)

Since the objective in Eq. (19) is a sum of independent cluster-level terms, minimizing this description length for contiguous clusters admits a dynamic programming algorithm solution jackson2005algorithm ; bellman2013dynamic ; PAY2022 that can identify the true optima in polynomial time.

The algorithm is constructed by recursing on $\mathcal{L}^{(i)}_{\mathrm{MDL}}$, the minimum description length of the first $i$ networks in $\mathcal{D}$ according to Eq. (19). Since the objective function decomposes as a sum over clusters, for any $j\in[1,S]$, the MDL can be calculated as

$\displaystyle\mathcal{L}^{(j)}_{\mathrm{MDL}}=\underset{i\in[1,j]}{\operatorname{min}}\left\{\mathcal{L}^{(i-1)}_{\mathrm{MDL}}+\mathcal{L}^{(\mathrm{cont})}_{k}([i,j])\right\},$

(20)

where we set the base case to $\mathcal{L}^{(0)}_{\mathrm{MDL}}=0$ and define $\mathcal{L}^{(\mathrm{cont})}_{k}([i,j])$ as the description length of the cluster of networks with indices $\{i,...,j\}$, according to Eq. (18) with the mode computed with the greedy procedure described in Sec. II.2. Once we recurse to $\mathcal{L}^{(S)}_{\mathrm{MDL}}$, we have found the MDL of our complete dataset, and keeping tab of the minimizing $i$ in Eq. (20) for every $j$ allows us to reconstruct the clusters.

In practice, the recursion can be implemented from the bottom up, starting with $\mathcal{L}^{(1)}_{\mathrm{MDL}}$, then $\mathcal{L}^{(2)}_{\mathrm{MDL}}$, and so on. The computational bottleneck for calculating $\mathcal{L}^{(j)}_{\mathrm{MDL}}$ is finding the modes of a cluster $j$ times for each evaluation of Eq. (20) (once for each $i=1,...,j$), leading to an overall complexity $O(jN\log N)$ for this step. Summing over $j\in[1,S]$, the overall time complexity of the dynamic programming algorithm is $O(S^{2}N\log N)$, which we verify numerically in Appendix B.

To demonstrate that the methods presented in Secs. II.2 and II.3 can effectively identify modes and clusters in network populations, we test their ability to recover the underlying modes and clusters generated from the heterogeneous population model introduced in ref. YKN22Clustering . We examine the robustness of these methods under varying noise levels that influence the similarity of the generated networks with the cluster’s mode.

The generative model in YKN22Clustering supposes (using different notation) that we are given $K$ modes $\mathcal{A}$ as well as the cluster assignments $\mathcal{C}$ of the networks $\mathcal{D}$. Each network $s\in C_{k}$ is generated by first taking each edge $(i,j)\in\bm{A}^{(k)}$ independently and adding it to $\bm{D}^{(s)}$ with probability $\alpha_{k}$ (the true-positive rate). Then, each of $M_{k}^{\ast}$ possible edges absent from $\bm{A}^{(k)}$ is added to $\bm{D}^{(s)}$ with probability $\beta_{k}$ (the false-positive rate). After performing this procedure for all clusters, the end result is a heterogeneous population of networks $\mathcal{D}$ with $K$ underlying modes, with noise in the networks $C_{k}$ surrounding each mode $\bm{A}^{(k)}$ determined by the rates $\alpha_{k}$ and $\beta_{k}$. The higher the true-positive rate $\alpha_{k}$ and the lower the false-positive rate $\beta_{k}$, the closer the networks in cluster $C_{k}$ resemble their corresponding mode $\bm{A}^{(k)}$.

Employing Bayesian inference of the modes and cluster assignments as in ref. YKN22Clustering involves adding prior probability distributions over the modes $\mathcal{A}$ and cluster assignments $\mathcal{C}$ to the heterogeneous network model YKN22Clustering . With a specific choice of priors on the modes and cluster sizes, Eq. (15) is precisely the equation giving us the Maximum A Posteriori (MAP) estimators of $\mathcal{A}$ and $\mathcal{C}$ in this model. We defer the details of this correspondence to Appendix D.

For our experiments, we use two modes, mode 1 and mode 3 from Fig. 1, as the planted modes $\mathcal{A}_{true}$ we aim to recover. To provide a single intuitive parameter quantifying the noise level in the generative model, we choose the true- and false-positive rates to satisfy $p=\beta_{1}=\beta_{3}=1-\alpha_{1}=1-\alpha_{3}$ for each run. Viewing the networks as binary adjacency matrices, the parameter $p$ corresponds to the probability of flipping entries of the matrix from 0 to 1 and vice-versa when constructing a network from its assigned cluster. We denote the parameter $p$ as the “flip probability” to emphasize this interpretation (same formulation as in ref. YKN22Clustering ). A flip probability $p=0$ corresponds to clusters of networks identical to the cluster modes, and a flip probability of $p=0.5$ corresponds to completely random networks with no clustering in the population. We thus expect it to be easy to recover the planted modes $\mathcal{A}_{true}$ and clusters $\mathcal{C}_{true}$ for $p=0$, and the problem becomes more and more difficult as we approach $p=0.5$.

We run three separate recovery experiments to test both the unconstrained and contiguous description length objectives in Eq. (15) and Eq (19), respectively. For the unconstrained objective, in each run we generate a population of $S$ networks from the model described above, with each network generated from either mode 1 or mode 3 at random with equal probability. We then identify the modes $\mathcal{A}_{\mathrm{MDL}}$ and clusters $\mathcal{C}_{\mathrm{MDL}}$ that minimize the objective in Eq. (15) using the merge-split algorithm detailed in Sec. II.2 and Appendix A. For the recovery of contiguous clusters, in one experiment we generate $S/2$ consecutive networks from each mode so that the population consists of $K=2$ adjacent contiguous clusters. And in another experiment, we generate $S/4$ networks from mode 1, $S/4$ networks from mode 3, and repeat this so that there are $K=4$ adjacent contiguous clusters of the $S$ networks generated from the two distinct modes. For these two experiments, we run the dynamic programming algorithm detailed in Sec. II.3 to identify the modes $\mathcal{A}_{\mathrm{MDL}}$ and clusters $\mathcal{C}_{\mathrm{MDL}}$ that minimize the objective in Eq. (19). In all three experiments, we generate a population of $S=100$ networks, each constructed from its corresponding mode using the single flip probability $p$ to introduce true- and false-positive edges.

To quantify the mode recovery error, we use the network distance quantified by the average Hamming distance between the inferred modes $\mathcal{A}$ and the planted modes $\mathcal{A}_{true}$. As both of our algorithms automatically select the optimal number of clusters $K$, the number of modes we infer can differ from the true number ($K=2$ or $K=4$, depending on the experiment). In each experiment, we therefore choose the $K$ inferred modes in $\mathcal{A}$ with the largest corresponding clusters and compute the average Hamming distance between these and the true modes in $\mathcal{A}_{true}$. (Since there are $K!$ ways to choose the inferred mode labels, we choose the labeling that produces the smallest Hamming distance.) To measure the error between our inferred clusters $\mathcal{C}$ and the planted clusters $\mathcal{C}_{true}$ (the “partition distance”), we use one minus the normalized mutual information vinh2010information . We also compute the inverse compression ratio (Eq. (17)) to measure how well the network population can be compressed. We pick a range of values of $p$ to tune the noise level in the populations, and at each value of $p$ we average these three quantities over $200$ realizations of the model to smooth out noise due to randomness in the synthetic network populations. We choose $K_{0}=1$ for these experiments, but this choice has little to no effect on the results (see Appendix B).

Figure 2 shows the results of our first reconstruction experiment. The reconstruction performance gradually worsens as $p$ increases due to the increasing noise level in the sampled networks relative to their corresponding modes (Fig. 2a). In all experiments, the network distance reaches that expected for a completely random guess of the mode networks—a 50/50 coin flip to determine the existence of each edge, denoted by the dashed black line—when $p=0.5$. The results in Fig. 2a indicate that in both the unconstrained and contiguous cases, our algorithms are capable of recovering the modes underlying these synthetic network populations with high accuracy, even for substantial levels of noise (up to $p\approx 0.3$, corresponding to an average of $30\%$ of the edges/non-edges differing between each mode and networks in its cluster).

The partition distance shows similar gradual performance degradation, with substantial increases in the distance beginning at $p\approx 0.3$ for the contiguous experiments and $p\approx 0.15$ for the unconstrained experiment (Fig. 2b). The partition distance levels off at different values across the three experiments, with the unconstrained case exhibiting significantly worse performance than the contiguous cases. We expect this result since contiguity simplifies the reconstruction problem by reducing the space of possible clusterings. Because information-based measures account for the entire space of possible clusterings instead of the highly constrained set produced by contiguous partitions, they overestimate the similarity of partitions in this constrained set. This overestimation intensifies with more clusters kirkley2022spatial .

The inverse compression ratio (Eq. (17)) for these experiments gradually approaches $1$ (no compression relative to transmitting each network individually, denoted by the dashed black line) as the noise level $p$ increases (Fig. 2c). This result is consistent with the intuition that noisier data will be harder to compress, while data with strong internal regularity will be much easier to compress, as the homogeneities can be exploited for shorter encodings. When $p$ is small, we can achieve up to $10$ times compression over the naive baseline by using the inferred underlying modes and clusters to transmit these network populations.

The results in Fig. 2 indicate that our algorithms can recover the underlying modes and their clusters in synthetic network populations. However, these results also depend on how distinguishable the underlying modes are. For identical modes, $\bm{A}^{(1)}=\bm{A}^{(2)}$, it is impossible to recover the cluster labels of the individual network samples $\bm{D}^{(s)}$. To investigate the dependence of the recovery performance on the modes themselves, we repeat the experiment in Fig. 2, except this time we systematically vary the mode networks $\mathcal{A}$ for each trial to achieve various levels of distinguishability. In each trial, we set $\bm{A}^{(1)}$ equal to the corresponding mode in Fig. 1 as before, but then generate the edges in $\bm{A}^{(2)}$ from $\bm{A}^{(1)}$ using the flip probability $\gamma$, which we call the “mode separation”. For mode separations $\gamma\approx 0$, it is challenging to recover the correct cluster labels of the individual sample networks because $\bm{A}^{(2)}$ will closely resemble $\bm{A}^{(1)}$. On the other hand, for mode separations $\gamma\approx 0.5$, the modes will typically be easily distinguished since $\bm{A}^{(2)}$ will have many edges/non-edges that have flipped relative to $\bm{A}^{(1)}$.

Figure 3 shows the results of this second experiment. The panels show the partition distance between the true and inferred cluster labels for a range of mode separations $\gamma$. In all experiments, the recovery becomes worse for lower values of the separation $\gamma$, but the algorithm still recovers a significant amount of cluster information even for relatively low $\gamma$. As in the previous set of experiments, the recovery performance is substantially worse for the discontiguous case compared with the contiguous cases, again due to the highly constrained ensemble of possible partitions considered by the partition distance in the contiguous cases.

In Appendix B, we show the recovery performance results for the network distance between the true and inferred modes as we vary the mode separation $\gamma$. For the mode recovery, the results are even more robust to the changes in mode separation. This result is consistent with the recovery performance in Fig. 2, where the recovery performance of the partitions starts to worsen at lower noise levels than the recovery of the modes. Thus, small perturbations in the inferred clusters may not affect the inferred modes much, since misclassified networks likely have little in common with the rest of their cluster.

For our first example with empirical network data, we study a collection of worldwide import/export networks. The nodes represent countries and the edges encode trading relationships. The Food and Agriculture Organization of the United Nations (FAO) aggregates these data, and we use the trades made in 2010, as in de2015structural . Each network in the collection corresponds to a category of products, for example, bread, meat, or cigars. We ignore information about the intensity of trades and merely record the presence or absence of a trading relationship for each category of products. The resulting collection comprises 364 networks (layers) on the same set of 214 nodes with 874.6 edges (average degree of 8.2) on average, with some sparse networks having as little as one edge and the densest containing 6,529 edges. These networks are unordered, so we employ the discontiguous clustering method described in Sec. II.2. We run the algorithm multiple times with varying initial number of clusters $K_{0}$ to find the best optima, although as with the synthetic reconstruction examples the choice of $K_{0}$ has little impact on compression. The best compression we find results in 8 modes and achieves a compression ratio of $\eta(\mathcal{D})=0.562$, indicating that it is nearly twice as efficient to communicate the data when we use the modal networks and their clusters. In contrast, in de2015structural a clustering analysis of the same network layers using structural reducibility—a measure of how many layers can be aggregated to reduce pairwise information redundancies among the layers—yielded 182 final aggregated layers, which would poorly compress the data under our scheme and not provide a significant benefit in downstream analyses due to the large final number of clusters. Key properties of the configuration of modes and clusters inferred by our algorithm are illustrated in Fig. 4.

In Fig. 4a we show the number of edges in each inferred mode, which indicates that these modes vary substantially in density to reflect the key underlying structures in networks within their corresponding clusters. The sizes of the clusters, shown in Fig. 4b, also vary substantially, with the most populated cluster (cluster 4) containing nearly 7 times as many networks as the least populated cluster (cluster 6). Some striking geographical commonalities and differences in the structure of the modes can be seen due to the varying composition of their corresponding clusters of networks. Figures 4c and 4d show the differences and similarities respectively between the structure of the modes for clusters 5 and 7, which are chosen as example modes because of their modest densities and distinct distributions of product types (Fig. 4e). Edges that are in mode 7 but not in mode 5 are highlighted in blue, while edges in mode 5 but not in mode 7 are highlighted in red. Meanwhile, the shared edges common to both networks are shown in Fig. 4d in black. Mode 5, which contains a diversity of product types and a relatively large portion of grain and protein products, has a large number of edges connecting the Americas to Europe that are not present in mode 7. On the other hand, mode 7, which is primarily composed of networks representing the trade of fruits and vegetables, has many edges in the global south that are not present in mode 5. However, both modes share a common backbone of edges that are distributed globally.

We categorized the 364 products (the network layers being clustered) into 12 broader categories of product types, plotting their distributions within each cluster in Fig. 4e. There are a few interesting observations we can make about this figure. Nearly all of the dairy products are traded within networks in a single cluster (cluster 3), indicating a high degree of similarity in the trade patterns for dairy products across countries. A similar observation can be made for live animals, which are primarily traded in cluster 4. On the other hand, many of the other products (grains, proteins, sweets, fruits, vegetables, and drinks) are traded in reasonable proportion in all clusters, which may reflect the diversity of these products as well as their geographical sources, which can give rise to heterogeneous trading structures. The densities of the modes and the sizes of the clusters do not have a clear relationship, with cluster 6 containing the smallest number of networks but the densest mode, and clusters 4 and 5 having sparser modes and much larger clusters. This reflects a higher level of heterogeneity in the structure of the trading relationships captured in cluster 6, which requires a denser mode for optimal compression, while the converse is true for clusters 4 and 5.

We also identify substantial structural differences in the inferred modes. In Appendix C, we compute summary statistics (average degree, transitivity, and average betweenness) for the modes output in this experiment and the network layers in their corresponding clusters. The statistics vary much more across clusters than within the clusters, suggesting that the MDL optimal mode configuration exemplifies distinct network structures within the dataset. Because the within-cluster average value of each statistic and the corresponding value for the mode network are similar, our method provides an effective preprocessing step for network-level regression tasks.

We conclude our analysis with a study of a set of networks representing global marine fauna over the past 500 million years. We aggregate fossil occurrences of the shelled marine animals, including bryozoans, corals, brachiopods, mollusks, arthropods, and echinoderms, into a regular grid covering the Earth’s surface rojas_multiscale_2021 . From these data, we construct unweighted bipartite networks representing $90$ ordered time intervals in Earth’s history (geological stages): An edge between a genus and a grid cell indicates that the genus was observed in the grid cell during the network’s corresponding geological stage. We also construct the unipartite projections of these networks: An edge from one genus to another indicates that these two genera were present in the same grid cell during the stage corresponding to the network. In total, there were 18,297 genus nodes, 664 grid cell nodes, 67,371.5 edges on average for the 90 unipartite graphs (average degree of 7.4), and 1462.2 edges on average for the bipartite graphs (average degree of 0.08, corresponding to an average of roughly 10 percent of genera being present at each layer).

In Fig. 5, we show the results of applying our clustering method for contiguous network populations (Sec. II.3) to both the unipartite and bipartite populations representing the post-Cambrian fossil record. We find clusters that capture the known large-scale organization of marine diversity. Major groups of marine animals archived in the fossil record are organized into global-scale assemblages that sequentially dominated oceans and shifted across major biotic transitions. Overall, the bipartite and unipartite fossil record network representations both result in transitions concurrent with the major known geological perturbations in Earth’s history, including the so-called mass extinction events. However, differences in the clusters retrieved from the unipartite and bipartite representations of the underlying paleontological data highlight the impact of this choice on the observed macroevolutionary pattern eriksson_how_2021 .

We also use our methodology to assess the extent to which the standard division of the post-Cambrian rock record in the geological time scale and the well-known mass extinction events compress the assembled networks. Specifically, we evaluate the inverse compression ratio in Eq. (17) on three different partitions of the fossil record networks that are defined by clustering the assembled networks into geological eras (Paleozoic, Mesozoic, and Cenozoic), geological periods (Ordovician to Quaternary), and six time intervals between the five mass extinctions in Fig. 5, with planted modes constructed by placing the networks into each cluster using the algorithm described in Sec. II.2 and Appendix A.

Table 1 shows the results of these experiments. All three partitions compress the fossil record networks almost as much as the optimal partition, which represents a natural division based on major regularities. Accordingly, the planted partition based on mass extinctions is almost as good as this optimal partition because mass extinctions are concurrent with the major geological events shaping the history of marine life. In contrast, partitions based either on standard geological eras or periods are less optimal, likely because they represent, to some extent, arbitrary divisions that are maintained for historical reasons. Our results here provide a complementary perspective to the work in rojas_multiscale_2021 , where a multilayer network clustering algorithm was employed that clusters nodes within and across layers to reveal three major biotic transitions from the fossil data. In Appendix E we review this and other existing multiplex and network population clustering techniques, discussing the similarities and differences with our proposed methods.