2.1 Data matrix and partition matrix

Consider a set of nodes indexed by $[n] = \{1,\dots,n\}$ , and an $n$ -by- $m$ data matrix $M$ in which row $i$ represents data associated with node $i$ , and entry $M_{ij}$ represents the value of the $j$ th data item associated with node $i$ . In a basic setting, $M$ equals the adjacency matrix $A$ (with $m=n$ ) of a single graph, and the rows correspond to adjacency profiles of the nodes. Alternatively, network data could be summarized by a distance matrix $D$ (with $m=n$ ), in which case the rows of the data matrix correspond to distance profiles. Indeed, matrices $D$ and $A$ provide equivalent representations of the network topology.Footnote ¹ In this article, we take a more general approach and allow the number $m$ of data items associated with a node to be arbitrary. This flexibility allows to model multilayer networks by concatenating, say, several adjacency or distance matrices side by side into a single data matrix. Furthermore, any data associated with nodes can easily be concatenated to the data matrix as extra columns. In this more general case, each row of $M$ corresponds to a data profile of a node.

Using the data matrix, we partition the node set into $k$ disjoint sets called communities. Such a partition can be represented as an $n$ -by- $k$ partition matrix $R$ with entries

\begin{equation*} R_{iu} \ = \ \begin {cases} 1 &\quad \text {if node $i$ is in community $u$}, \\ 0 &\quad \text {else}. \end {cases} \end{equation*}

In applications, some entries of the data matrix $M$ may be undefined or unobserved. These situations are handled by equipping $M$ with an $n$ -by- $m$ indicator matrix $C$ in which $C_{ij} =1$ indicated a valid entry, and $C_{ij}=0$ indicates an undefined or unobserved entry. The corresponding matrix elements of $M$ , which are not defined or are missing, are replaced by dummy values, which is chosen to be $0$ in all our sample cases. We demonstrate how this works in Section 5.

Our aim is to group nodes into few communities in such a way that description length of $M$ is minimized based on a probabilistic model for the matrix elements of $M$ . In other words, we try to find an optimal number of communities $k$ and corresponding optimal partitioning $\mathcal{U}_k=\{U_1,U_2,\dots,U_k\}$ of $[n]$ to achieve this. In case of non-negative integer-valued $M$ , rows in a community $s \in [k]$ are modeled by a sequence of $m$ independent Poisson-distributed random variables denoted as $(X^s_1,X^s_2,\cdots,X^s_m)$ . As a result, the probabilistic model constitutes of $n \times m$ independent Poisson random variables with $k \times m$ parameters, the expectations of the corresponding variables. Such a choice is adopted from (Reittu et al., Reference Reittu, Bazsó, Weiss, Fränti, Brown, Loog, Escolano and Pelillo2014). Communities are selected in order to minimize the magnitude

\begin{equation*} L(\mathcal {U}_k) \,{:\!=}\, - \sum _{s \in [k]} \sum _{d\in U_s } \sum _{j \in [m]} \log _2\!(\mathbb {P}(X^s_j = M_{dj}). \end{equation*}

According to classical information theory (e.g., Cover & Thomas, Reference Cover and Thomas2006), there exists a binary code, the Shannon code, for encoding $M$ with code length $= \lceil L \rceil$ . We use $L$ as a cost function of the partition $\mathcal{U}_k$ . The quality of the found communities is the compression ratio

\begin{equation*} r(\mathcal {U}_k)=\frac {L(M)}{L(\mathcal {U}_k)}, \end{equation*}

in which $L(M)=\sum _{i,j: M_{ij}\gt 0} \log M_{ij}$ is the number of bits needed to represent data matrix $M$ as a string of integers.

2.2 Likelihood function

We employ a statistical latent-variable model in which all observable entries (those with $C_{ji}=1$ ) of the data matrix $M$ are conditionally independent and Poisson-distributed random variables given the community structure. Furthermore, all rows of $M$ corresponding to the same community are identically distributed. This statistical model is parameterized by an $n$ -by- $k$ partition matrix $R$ and an $m$ -by- $k$ expectation matrix $\Lambda$ of Poisson variables and corresponds to likelihood function

(1) \begin{equation} f(M | \Lambda,R) \ = \ \prod _{j=1}^{n} \prod _{i=1}^{m} \prod _{u=1}^k \operatorname{Poi}(M_{ji} | \Lambda _{i u})^{C_{ji} R_{ju}}, \end{equation}

where $\operatorname{Poi}(x | \lambda ) = e^{-\lambda } \frac{\lambda ^x}{x!}$ is the probability mass function of a Poisson distribution with mean $\lambda$ . The corresponding log-likelihood can be written as

\begin{equation*} \label {log} \log f(M | \Lambda,R) \ = \ \sum _{i=1}^m \sum _{j=1}^n \sum _{u=1}^k R_{ju} C_{ji} \Big ( M_{ji} \log \Lambda _{iu} - \Lambda _{iu} \Big ) - \text {const}(M), \end{equation*}

where $\text{const}(M) = \sum _{j=1}^n \sum _{i=1}^m C_{ji} \log \!( M_{ji}!)$ does not depend on the model parameters and can be ignored. The above model is structurally similar to the SBM in which the data matrix has $m=n$ columns and corresponds to the adjacency matrix, and the Poisson distributions are replaced by Bernoulli distributions (Holland et al., Reference Holland, Laskey and Leinhardt1983; Zhang & Zhou, Reference Zhang and Zhou2016). In contrast to SBMs, the above model allows more flexibility in choosing data matrices with an arbitrary number of columns $m$ . Note also that only the rows are grouped in blocks, all columns are treated as separate. In this sense, the model has maximal number of variables with respect to the number of columns.

2.3 Maximum likelihood estimation

Having observed a data matrix $M$ , maximum likelihood estimation searches for a partition matrix $R$ of $[n]$ for which the function in Equation (2.2) is maximized. For any fixed $R$ , the $\Lambda$ -parameters are set equal to

(2) \begin{equation} \hat \Lambda _{iv}(R) \ = \ \frac{\sum _{j=1}^n M_{ji}R_{jv}C_{ji}}{\sum _{j=1}^n R_{jv}C_{ji}}, \end{equation}

which is the observed average of the $i$ th data item in community $v$ . As a consequence, a maximum likelihood estimate of $R$ is obtained by maximizing the profile log-likelihood $f(M \, | \, \hat \Lambda (R), R)$ , or equivalently minimizing

(3) \begin{equation} L(R) \ = \ \sum _{i=1}^m \sum _{j=1}^n \sum _{v=1}^k R_{jv}C_{ji} \left ( \hat \Lambda _{iv}(R) - M_{ji} \log \hat \Lambda _{iv}(R) \right ) \end{equation}

subject to $n$ -by- $k$ partition matrices $R$ , in which $\hat \Lambda _{iv}(R)$ is given by (2). We note that the above function can be written as $L(R) = \sum _{j=1}^n \ell _{j Z_j}(R)$ , in which

(4) \begin{equation} \ell _{jv}(R) \ = \ \sum _{i=1}^m C_{ji}\left ( \hat \Lambda _{iv}(R) - M_{ji} \log \hat \Lambda _{iv}(R) \right ) \end{equation}

is a normalized minus log-likelihood of the data vector of node $j$ , given that $j$ is placed in community $v$ and the rate parameters are equal to $\hat \Lambda (R)$ . We note that $L(R)$ , up to an additive constant, equals to the description length of the data matrix $M$ in the sense of Shannon coding.

The same algorithm can be used also in case of data matrix with positive real values, which was already shown in (Reittu et al., Reference Reittu, Bazsó, Weiss, Fränti, Brown, Loog, Escolano and Pelillo2014). In this case, each matrix element $M_{ij}$ is treated as a parameter (the expectation) of a Poisson distribution. The $\Lambda$ -matrix is computed according to Equation (2) for each partition $R$ . Equation (3) also remains intact, and $L$ equals to Kullback–Leibler divergence between the corresponding Poisson distributions, the original with parameters from data matrix $M$ , and those with parameters from Equation (2). The task is to minimize $L$ which can be done with the same algorithm as in the integer case.

2.4 Regular decomposition algorithm

Minimizing the cost function in Equation (3) with respect to $R$ is a hard nonlinear discrete optimization problem with an exponentially large input space of the order of $\Theta (k^n)$ , making exhaustive search computationally infeasible. This is why we suggest solving the problem using a greedy Algorithm 1 which is an expectation maximization type algorithm which alternates between updating $\Lambda$ according to Equation (2) for a fixed $R$ (E-step) and updating the partition $R$ by greedily updating the community of each node, one by one, to minimize $\ell _{jv}(R)$ in Equation (4) for a fixed $\hat \Lambda$ (M-step). Starting from a uniformly random initial assignment $R_0$ , the algorithm finds a local optimum as a limit of the greedy algorithm. Running the greedy algorithm for several initial random states $R_0$ , and selecting the community assignment with smallest cost in Equation (3) as the final output.

Figure 1. Pseudocode for the regular decomposition algorithm according to Reittu et al. (Reference Reittu, Leskelä, Räty and Fiorucci2018).

The runtime of Algorithm 1 is $O(s_{\textrm{max}} t_{\textrm{max}} k (n+m))$ , where $s_{\textrm{max}}$ is the number of random initializations and $t_{\textrm{max}}$ is the number of iterations per optimization round. Especially, the runtime is linear in $n$ and $m$ for bounded $k,s_{\textrm{max}},t_{\textrm{max}}$ and is hence scalable for large data sets.

2.5 Boosted regular decomposition

Setting up the data matrix for Algorithm 1 may require costly preprocessing, for example computing distances between node pairs (see Section 2.7). The matrix $M$ may also be simply too large to be treated as a whole, a situation frequently encountered in the realm of so-called big data. In such cases, the algorithm can be boosted by replacing the data matrix $M$ by a submatrix $M_{VW}$ with a row set $V \subset [n]$ of size $n_0$ and a column set $W \subset [m]$ of size $m_0$ , and running Algorithm 1 with the submatrix $M_{VW}$ as input. This results in a partition matrix $R^*$ of the node set $V$ . A community assignment of the remaining node set is then computed in a subsequent classification phase where the community index of each node $j \in [n] \setminus V$ is chosen to be $v(j) \in [k]$ :

(5) \begin{eqnarray} v(j)=\text{arg min}_{v'\in [k]}\ell _{jv'}(R^*),\quad \ell _{jv'}(R^*) \ = \ \sum _{i \in W} C_{ji}\left ( \hat \Lambda _{iv'}(R^*) - M_{ji} \log \hat \Lambda _{iv'}(R^*) \right ), \end{eqnarray}

where

\begin{equation*} \hat \Lambda _{iv}(R^*) \ = \ \frac {\sum _{j \in V} M_{ji}R_{jv}^*C_{ji}} {\sum _{j \in V} R_{jv}^*C_{ji}} \end{equation*}

is the observed average value of data item $i \in W$ among the nodes of $V$ classified into community $v \in [k]$ according to $R^*$ .

The runtime of Algorithm 1 applied to the submatrix $M_{VW}$ is $O(s_{\textrm{max}} t_{\textrm{max}} k (n_0+m_0))$ , and the runtime of the subsequent classification phase is $O(km_0 n)$ . Hence, the boosted RD algorithm has complexity $O(s_{\textrm{max}} t_{\textrm{max}} k (n_0+m_0) + k m_0 n)$ . The feasibility of this boosting approach requires that the row set $V$ is large enough to contain nodes from all communities, and the column set $W$ is a sufficiently informative collection of data items. A simple way of selecting $V$ and $W$ is by random sampling. This approach was developed in (Reittu et al., Reference Reittu, Leskelä, Räty and Fiorucci2018, Reference Reittu, Norros, Räty, Bolla and Bazsó2019) in which sufficient sample sizes were estimated and convergence proved in some model cases.

2.6 Estimating the number of communities

Algorithm 1 requires the number of communities $k$ as an input parameter. However, in most situations this parameter is not a priori known and needs to be estimated from the observed data. The problem of estimating the number of communities can be approached by recasting the maximum likelihood problem in terms of the MDL principle (Rissanen, Reference Rissanen1983; Grünwald, Reference Grünwald2007) where the goal is to select a model which allows a minimum coding length for both the data and the model. MDL adheres to the principle of Occam’s razor in which the best hypothesis follows the best compression of data, hence justifying the selection of MDL for this task.

When restricting to the model described in Section 2.2, then the $R$ -dependent part of the coding length equals $L(R)$ given by (3), and an MDL-optimal partition $R^*$ for a given $k$ corresponds to the minimal coding length

\begin{equation*} R^* = \text {arg min}_{R} L(R). \end{equation*}

It is not hard to see that $L(R^*)$ is monotonously decreasing as a function of $k$ , and a balancing term, the model complexity, is added to select the model that best explains the observed data. The model complexity is the length of a code that uniquely describes the mode itself. In all of our experiments, $L(R^*)$ (the negative log-likelihood) as a function of $k$ becomes essentially a constant above some value $k^*$ . Such an elbow point $k^*$ is used as an estimate of $k$ in the experiments in this article, see also (Ketchen & Shook, Reference Ketchen and Shook1996). In general, it might be necessary to have a more sophisticated method using a model complexity term (Reittu et al., Reference Reittu, Bazsó and Norros2017a; Norros et al., Reference Norros, Reittu and Bazsó2022; Peixoto, Reference Peixoto2012). However, in examples we are using it suffices to use a simplified version of the MDL principle based on the elbow point.

2.7 Using distances as data items in multiplex networks

In a multiplex network consisting of $s$ directed graphs with a common node set, each graph represents one layer. The distance matrices of the layers are denoted by $D_1, \ldots, D_s$ . If layer $r$ contains no path from node $i$ to node $j$ , we declare the corresponding entry as missing by setting $(C_r)_{ij} = 0$ , and we may define $(D_r)_{ij} = 0$ without loss of generality. As a result, we obtain $s$ indicator matrices $C_1, \dots, C_s$ . When layers are undirected, data about path lengths are encoded in a concatenated data matrix

(6) \begin{equation} M \ = \ [D_1, D_2, \dots, D_s], \end{equation}

and the corresponding indicator matrix is $C = [C_1, C_2, \dots, C_s]$ . In case of directed layers, data about directed path lengths are encoded in matrix

(7) \begin{equation} M = \left[D_1, D_1^T, D_2, D_2^T, \dots, D_s, D_s^T\right], \end{equation}

where row $i$ of matrix $D_r$ (resp. $D_r^T$ ) contains the shortest directed path lengths from $i$ to other nodes (resp. from other nodes to $i$ ). The corresponding indicator matrix is denoted $C = \left[C_1, C_1^T, C_2, C_2^T, \dots, C_s, C_s^T\right].$ The data matrix $M$ and the indicator matrix $C$ are then given as input to Algorithm 1, and the optimal number of communities is determined as in Section 2.6.

Computing the distance matrix in an unweighted directed graph with $n$ nodes and $e$ links has complexity $O(n(n+e))$ using breadth-first search (Bang-Jensen & Gutin, Reference Bang-Jensen and Gutin2009). The RD algorithm based on distances can be boosted by computing distances only for a restricted set of reference nodes $W \subset [n]$ of size $m_0$ , resulting in an $n$ -by- $m_0$ distance matrix $D$ in which $D_{ij}$ equals the distance from node $i \in [n]$ to node $j \in W$ , and $D$ can be computed using breadth-first search in $O(m_0(n+e))$ time. The same complexity bound is valid for concatenated data matrices of form (6)–(7) in multilayer networks with a bounded number of layers. Then, we may apply the boosted RD algorithm (Section 2.5) with $V=[n]$ and $W$ as above. The total complexity of the preprocessing step (restricted distance matrix computations) and boosted RD algorithm then equals $O(s_{\textrm{max}} t_{\textrm{max}} k (n+m_0) + k m_0 n + m_0(n+e))$ . For bounded number of communities $k$ and bounded iteration parameters $s_{\textrm{max}},t_{\textrm{max}}$ , this bound is $O(m_0(n+e))$ . Especially, by selecting $m_0$ constant, we obtain a scalable algorithm for sparse massive networks with $e=O(n)$ , capable of identifying communities in linear time with respect to the number of nodes $n$ .

3.1 Poissonian power-law graph

A simple random graph model can be induced from a random graph process described in (Norros & Reittu, Reference Norros and Reittu2006; van der Hofstad, Reference van der Hofstad2017) which we call a Poissonian power-law graph. We first sample node attributes $\lambda _1, \dots, \lambda _n$ independently at random from a probability distribution on the non-negative reals and thereafter connect each unordered node pair $ij$ by a link with probability $1-\exp \{-\lambda _i\lambda _j/\lambda \}$ , independently of other node pairs, where $\lambda = \sum _i \lambda _i$ . When the node attributes are distributed according to a power-law with density exponent $\tau \in (2,3)$ , we obtain a generator of a sparse random graph where the degree distribution has a finite mean and infinite variance.

In such power-law graphs, quite a rich topological structure spontaneously arises in the limit of large graph size and with probability tending to one. The nodes are categorized into sets called tiers according to their degree, such groping we call “soft hierarchy” (Norros & Reittu, Reference Norros and Reittu2006; Reittu & Norros, Reference Reittu and Norros2004; Norros & Reittu, Reference Norros and Reittu2008b,a). The top tier $V_0$ is formed, asymptotically as $n\rightarrow \infty$ , by nodes with degrees in the range $(n^{1/2}, \infty )$ , and the other tiers $V_k$ are formed by nodes with degrees in range $(n^{\beta _k}, n^{\beta _{k-1}}]$ , in which $\beta _0 = 1/2$ and $\beta _k=(\tau -2)^k/(\tau -1)+o(1)$ for $k \ge 1$ . For large values of $n$ , it is known that with high probability, the top tier $V_0$ is fully connected, and further, every node in $V_k$ has a link to $V_{k-1}$ for all $k$ up to order $\log \log n$ (Norros & Reittu, Reference Norros and Reittu2006; van der Hofstad & Hooghiemstra, Reference van der Hofstad and Hooghiemstra2008; Reittu & Norros, Reference Reittu and Norros2002). The subgraph induced by the union of the tiers $V_k$ for $0 \le k \le \log \log n$ is called the core network. Most of the nodes have small degrees and, as a result, are outside the core. However, any node is at a very short distance, compared to $\log \log n$ , from the bottom tier. This explains ultra-short distances in the graph.

Our aim is to show, through experiments on synthetic and real data, that our version of the RD algorithm can identify soft hierarchical structures by using network distances and node degrees together as a data matrix. The soft hierarchy is a compact description of the organization of shortest paths between most of the nodes in power-law graphs of the type we are interested in. RD compresses the distance matrix, and that is why it is likely that a soft hierarchy shows up in the resulting short description of the matrix. Degrees are also essential in describing the soft hierarchy, and that is why their inclusion in data matrix should help. This is why we argue that a degree-augmented and distance-based community structure matches qualitatively with a kind of rough soft hierarchy in synthetic and real power-law graphs. Notably, in expectation $\mathbb{E}{\lvert \cup _k V_k\rvert }/n\rightarrow 0$ , see (Norros & Reittu, Reference Norros and Reittu2006; Reittu & Norros, Reference Reittu and Norros2002) as $n \to \infty$ , which means that communities associated with the layers are very small and which means that such communities are undetectable for typical community detection algorithms assuming comparable community sizes.

This illustrates how our method should be used. At first, there should be an intuitive understanding of which data are essential for the problem to be solved. In the current example, the problem is to find the soft hierarchy as a community structure. In other problems, the data matrix could be completely different.

3.2 Regular decomposition using distances

We generated a power-law graph with ten thousand nodes using the model in Section 3.1 with node attributes drawn from a power-law distribution with density exponent $\tau =2.5$ , and we extracted its largest connected component as input for subsequent analysis. As a result, we have a graph with $n=7775$ nodes and adjacency matrix shown in Figure 2(a). Next, we computed the distance matrix of this graph. We identified communities of the graph by applying Algorithm 1 using the distance matrix as data matrix. By experimenting with different values of the number of communities $k$ , we found that the cost function in Section 2.6 saturates at $k=5$ . This value is identified as the most informative number of communities.

Figure 2. The adjacency matrix of a synthetic power-law graph ( $n=7775$ ) with rows and columns organized according to (a) random order of nodes, (b) communities identified by distance matrix, (c) communities identified by distance matrix and degrees (right).

The block structure of the adjacency matrix induced by the identified communities is shown in Figure 2(b). A clear block structure is revealed with one large block with relatively high density. All five identified communities are rather large. Hence, the identified community structure differs remarkably from the theoretical tier structure (Section 3.1) where the top layers are dense and small. This is a natural consequence of partitioning the graph using distance profiles because the small-degree neighbors of high-degree nodes are likely to have similar distance profiles with each other.

The degrees of the graph are plotted in Figure 3(a). The linear shape in the log-log plot on the left is typical for power-law graphs. The degrees of nodes in the five identified communities are shown in Figure 3(b). We see that all high-degree nodes are in the same community, but this community also contains nodes of smaller degree and is quite large. Furthermore, Figure 4(a) displays the graph with the five identified communities, plotted using Mathematica’s CommunityGraphPlot tool. We see that nodes in the red community are central for connecting nodes in the graph. The ring-like communities around the center appear to roughly play a role similar to tiers in a soft hierarchy, in that most shortest paths between peripheral nodes use the rings to reach the red center.

Figure 3. Degrees of nodes in a synthetic power-law graph ( $n=7775$ ) sorted from largest to smallest within each community (indicated by color). (a) Full graph viewed as one community. (b) Nodes organized into communities identified by distance matrix. (c) Nodes organized into communities identified by distance matrix and degrees.

Figure 4. Topology of a synthetic power-law graph ( $n=7775$ ) colored according to the community structure identified by the regular decomposition algorithm with (a) distance matrix, (b) distance matrix and degrees.

3.3 Regular decomposition using distances and degrees

We continue experimenting with partitioning the same graph sample as in the previous section. Instead of using the distance matrix with 7775 columns as in Section 3.2, we will now use a data matrix with only 101 columns, consisting of 100 randomly sampled columns of the distance matrix and 1 additional column containing the node degrees. The aim of this experiment is twofold. First, we wish to investigate how adding the degrees to the data matrix affects the inferred community structure. Second, we will demonstrate that computing distances to a relatively small set of reference nodes suffices to well characterize the distance profiles of most nodes.

The adjacency matrix organized according to five identified communities using the modified data matrix is shown in Figure 2(c). The main difference from the previous case shown in the middle of the same plot is a small central community with high-degree nodes only. This can be seen in Figure 3(c) which presents the node degrees grouped by communities. The blue community contains all high-degree nodes, and its degree sequence does not overlap with other communities. Nodes in the blue community may hence be thought as tier 1 nodes. As a result, the community structure is qualitatively closer to the theoretical tiers of the power-law graph. According to the theory, there are only $\Theta (\!\log \log n)$ tiers, we may expect only a couple of layers in our sample with $\log \log n = 2.19$ .

Figure 4(b) visualizes the identified communities from a topological point of view. The small and dense red community in the center corresponds qualitatively to the top tier of the theoretical power-law graph structure in Section 3.1. The second and third largest communities can be seen as tier 2 and tier 3 communities, and the remaining communities form the periphery of low-degree nodes. Incorporating degrees in the data matrix can hence substantially change the community structure and in our case align the communities better with a soft hierarchy of nodes.

4.1 Regular decomposition using distances

We partition the AS graph into communities by Algorithm 1 using $100$ randomly sampled columns of the distance matrix as data matrix. This appears sufficient for this type of network where we expect that the distances from a typical reference node depend heavily on the position of the node in the network hierarchy. Using the method in Section 2.6, we found that the most informative number of communities is $k=10$ . Figure 5(b) displays the adjacency matrix of the AS graph organized by the identified community structure, showing that the communities are all rather large and of comparable size. The link densities inside and between communities are all low and comparable to the overall link density. The degrees of nodes grouped into communities are displayed in Figure 6(b).

The results for the AS graph have similarities with the synthetic power-law graph in Section 3.2. For instance, although all high-degree nodes are in the same community, this community also contains many low-degree nodes and thus has a low internal density. This can be seen in Figure 6(b) where the degree sequences of different communities overlap. As a result, using only graph distances as the data matrix, we were not able to identify the tiers of the AS graph. For instance, the smallest community has much more nodes than there are tier 1 nodes³.

4.2 Regular decomposition using distances and degrees

We repeat the community identification experiment of the AS graph by augmenting the 100-column data matrix used in Section 4.1 with 1 column containing the node degrees. Again, $k=10$ is identified as the most informative number of communities. Figure 5(c) displays the adjacency matrix of the graph organized according to the identified communities. The resulting community structure substantially differs from the one in the previous section. The smallest two communities in Figure 5(c) are dense and approximately correspond to tier 1 and tier 2 subnetworks of the AS graph. The degree sequences of the communities are shown in Figure 6(c). There are three rather small and dense communities which contain all high-degree nodes, but no low-degree nodes.

To assess the quality of the identified communities, we determined the AS identities in the three smallest and densest communities using a list of AS networks³ and an AS lookup tool.Footnote ⁴ Our aim is to verify that the three identified communities are close to tier 1–3 subnetworks and contain the most important AS. The smallest identified community (Figure 7) contains $36$ nodes, out of which $23$ are tier 1 and the rest are tier 2. The members of tier 2 are important telecom carriers. For example, “Hurricane Electric”³ has a very high degree ( $7061$ ), which explains why it is included in the smallest community by the RD algorithm. The discrepancy between tier 1 and the smallest identified community may be due to the fact that the AS graph topology deviates from a simple power-law graph (e.g., Chen et al., Reference Chen, Chang, Govindan, Jamin, Shenker and Willinger2002). The second smallest identified community is shown in Figure 8(a). The largest degrees are around $200$ . In this community, the ten nodes of highest overall degree do not belong to tier $1$ , but nevertheless correspond to networks operated by major companies such as British Telecom (UK) and Microsoft (US). The third smallest community, displayed in Figure 8(b), is much sparser, and its highest degrees are around $20$ . Top-degree members are Telefonica Data S.A. (BR), Orano (US), and Harvard University (US).

Figure 7. Subgraph of the AS graph induced by the smallest community identified by regular decomposition using distances and degrees.

Figure 8. Subgraphs of the AS graph induced by the second (a) and third (b) smallest communities identified by regular decomposition using distances and degrees. Ten nodes of highest overall degree in the second smallest community are highlighted in red: 1. Net Access Corporation, 2. Microsoft, 3. London Interconnection Point, 4. BT, 5. Internet Initiative Japan, 6. Frontier Communications of America, 7. MCI Communications Services, 8. INAP, 9. TransTeleCom, and 10. Telstra.

We conclude that augmenting the graph distance matrix by a column containing the node degrees allows to identify much more meaningful communities, compared to only using the distance matrix. The RD method was able to identify central carriers in the top tiers with good accuracy from a large data set. In particular, we discovered a dense tier 1-rich subnetwork. The suggested method could be used even for extremely large graphs encountered in areas such as biology and social networks, where it might be impossible to acquire the entire graph for analysis. Our methods need only a limited sample of shortest paths between a set of sampled nodes and node degrees. Community detection on such a sample results in a model which can be used to classify any other node outside the sample. It has the potential to rapidly detect soft hierarchies in massive networks.

5.1 Regular decomposition using layerwise distances

Using Algorithm 1 and the method in Section 2.6, we discovered $k=6$ as the most informative number of communities and identified the corresponding communities. The resulting data matrix, organized by the identified communities, is shown at the bottom of Figure 10. The smallest community has 9 nodes consisting of airports in the Middle East (4), Europe (2), East Asia (1), and South America (1). The second smallest community has 13 nodes consisting of airports in East Asia (5), Africa (3), South-East Asia (2), Pacific (2), and North America (1). The remaining four communities are of comparable size: One of them has 27 airports in France and 85 in USA, and the other three have most airports located in the USA.

6.1 Community detection

6.1.1 Internet autonomous systems graph

As stated in the introduction, our aim is to have a community detection method which identifies communities which reflect various aspects of data associated with the network and their role in the graph topology. There are of course many existing powerful community detection methods which can do this in some particular cases. We illustrate this point by finding communities in the AS graph (Section 4) using two popular methods: modularity maximization and SBM fitting.

Modularity maximization (Newman, Reference Newman2006) aims to find densely connected communities which have as little as possible links between the communities. Figure 12(a) displays the adjacency matrix organized according to the 25 identified communities, and Figure 13 illustrates the subgraphs induced by the communities. As expected, the community structure determined by modularity maximization is substantially different from the community structures identified by RD in Figure 5(b–c). The subgraphs in Figure 13 do not respect the tiered structure found with RD. For instance, the high-degree nodes forming tier 1 are embedded in very large communities, which can be inferred from Figure 13.

Figure 11. Communities of the world airline graph identified by (a) regular decomposition using layerwise distances, (b) SBM fitting with a layer-aggregated adjacency matrix.

Figure 12. Adjacency matrix of the AS graph organized according to (a) 25 communities identified by modularity maximization, (b) 10 communities identified by SBM fitting.

Figure 13. Subgraphs of the AS graph induced by the $25$ communities identified by modularity maximization.

SBM fitting (Zhao et al., Reference Zhao, Levina and Zhu2012) uses the adjacency matrix to find communities in which relations inside and between the communities are like those in classical random graphs. Information-theoretic model fitting can be used to find the communities. We used the RD method for this. For computational tractability, instead of the full adjacency matrix we restricted to the largest connected component of the subgraph induced by a uniform random sample of 10,000 nodes. After identifying the communities of the restricted graph, the parameters of the model were used to classify all nodes of the graph. A smaller sample is not feasible, because the resulting induced subgraph would hardly have any links. This is in stark contrast to RD using graph distances, where a sample of 100 nodes suffices. Figure 12(b) shows the resulting community structure, where a dense tier 1-like community appears, but a deeper tier hierarchy seems quite weakly represented, as manifested by a large unstructured block of low-degree nodes. In contrast, a much finer community structure is identified by RD based on distances and degrees in Figure 5(c). A typical subgraph induced by a pair of communities identified by SBM fitting is a well-connected graph in which one part is like a neighborhood of the other community, see Figure 14. On the other hand, in a community structure identified by distance profiles, such pairs are typically not connected—this is simply because many distances are larger than one.

Figure 14. A subgraph of the AS graph induced by two distinct communities (red, blue) identified by SBM fitting. The nodes of the blue community are almost entirely low-degree neighbors of the nodes in the red community.

For a more quantitative comparison, we computed the PageRank centrality of the network nodes (with damping factor 0.8) and plotted in Figure 15 the PageRank distributions within communities identified by three methods: modularity maximization, RD using distances, and RD using distances and degrees. We see that the latter is the only method able to separate nodes of high PageRank into a common community. This illustrates the ability of our method in identifying community structures associated with the topological roles of nodes in the network.

Figure 15. PageRank centrality of AS graph nodes in communities identified using three alternative methods and sorted from smallest to largest within each community (indicated by color). (a) Modularity maximization, (b) RD using distances, (c) RD using distances and degrees. Only the last method can group the most central nodes into separate communities, indicated by blue and orange circles in (c).

6.1.2 World airline graph

We compare the community structure of the world airline graph determined in Section 5.1 with a more customary SBM fitting applied to the adjacency matrix of the undirected graph obtained by collapsing the layers and ignoring link directions. We impose the same number $k=6$ of communities as previously. The resulting community structure, visualized in Figure 11(b), bears some similarity with the one in Figure 11(a), but mostly in the periphery of the network which forms the largest community in both cases. The peripheral communities are highlighted in Figure 16. Furthermore, Figure 17 displays the overlap matrix of communities. The overlap is weak, and the community structures are quite different. If the community structures were similar, there should be at each row exactly one strong element in this matrix and those strong elements would be all in different columns.

Figure 16. Largest community (highlighted in red) of the world airline graph identified using (a) regular decomposition using layerwise distances, (b) SBM fitting with a layer-aggregated adjacency matrix.

Figure 17. Overlap of communities determined by SBM fitting using layer-aggregated graph (row index) and by regular decomposition using layerwise distances (column index). The $(i,j)$ -entry is computed as the number of common nodes in community pair $(i,j)$ divided by the number of nodes in the larger community in the pair. Only the pair $(3,6)$ indicates a substantial value around $0.5$ . This pair corresponds to the peripheral communities highlighted in Figure 16.

6.2 Data compression

A natural measure of our method is the compression rate of the data matrix used in community detection. This experiment is done on the distance matrix for the AS graph. As the uncompressed description length, we used Mathematica’s internal memory requirement for storing a distance matrix (ByteCount), which for the AS graph equals $L_0 \approx 1.7 \times 10^{10}$ bits. First, as a standard compression, we applied Mathematica’s built-in implementation of the zlibFootnote ⁶ algorithm. Second, we computed the description length $L$ for the distance matrix by applying RD. This consists of the Shannon code length, which is minus base-2 logarithm of the probability of the distance matrix in the probabilistic model induced by the community structure, and the code lengths of the parameters. We used the leading part of the prefix code lengths of integers (Grünwald, Reference Grünwald2007). For a positive integer $m$ , the code length of the integer is $ \lceil \log _2 m \rceil + \lceil \log _2\log _2 m \rceil +\cdots$ , in which $\log _2$ is iterated as long as the result remains positive, after which the sum is truncated. The parameters we need to encode are the expectations of $k\times n$ Poisson variables, $k=10$ , and $n=22963$ , and the partition into communities which can be represented as an $n$ -vector with coordinates in $\{1,\dots,10\}$ . Third, we computed similar code length for community structure found using node degrees along with the distance matrix. The compression ratio $L/L_0$ is displayed in Table 1. We see that zlib compresses the original data by around $9$ times while RD using distances does a better job with around $12$ times compression. Augmenting the RD by also using degrees leads to a slight further improvement.

Table 1. Compression ratio for the AS graph and the world airline graph using three methods: zlib compression algorithm, RD using distances, and RD using distances and degrees

We also repeated this experiment for the world airline graph for which the uncompressed layerwise distance matrix requires $L_0 \approx 5.5 \times 10^8$ bits. The zlib algorithm compresses this data by around $80$ times, while RD on the layerwise distance matrix is able to compress $235$ times, almost three times more than zlib, see Table 1.