Modularity in Multilayer Networks using Redundancy-based Resolution and Projection-based Inter-Layer Coupling
Alessia Amelio Giuseppe Mangioni Andrea Tagarelli




 

1]Alessia Amelio
2]Giuseppe Mangioni
1]Andrea TagarelliCorresponding author: Andrea Tagarelli. 

[1]Dept. of Computer Engineering, Modeling, Electronics and Systems Engineering (DIMES), University of Calabria, Rende (CS), Italy
 Email: {aamelio, tagarelli}@dimes.unical.it.

[2]Dept. of Electrical, Electronics and Computer Engineering, University of Catania, Catania, Italy
 Email: giuseppe.mangioni@dieei.unict.it.  


  Modularity in Multilayer Networks using 
Redundancy-based Resolution and
 Projection-based Inter-Layer CouplingAccepted for publication with IEEE Trans. on Network Science and Engineering, April 2019. DOI:  An abridged version of this work appeared in <cit.>.
    [
    
===================================================================================================================================================================================================================================================================






The generalized version of modularity for multilayer networks, a.k.a.  multislice modularity, is characterized by two model parameters, namely   resolution factor and  inter-layer coupling factor. The former corresponds to  a notion of layer-specific relevance, whereas the inter-layer coupling factor represents the strength of node connections across the network layers. Despite the potential of this approach, the setting of both parameters   can be arbitrarily selected, without considering specific characteristics from the topology of the multilayer network as well as from an available community  structure. Also, the multislice modularity is not designed to explicitly model order relations over the layers, which is of prior importance for dynamic networks. 

This paper aims to overcome the main limitations of the multislice modularity by introducing a new formulation of modularity for multilayer networks. We revise the role and semantics of both the resolution and inter-layer coupling factors based on information available from the within-layer and inter-layer structures of the multilayer communities. Also, our proposed multilayer modularity is general enough to consider orderings of network layers and their constraints on layer coupling.  Experiments  were carried out on synthetic and  real-world multilayer networks using   state-of-the-art approaches for multilayer community detection. The obtained results have shown the meaningfulness of the proposed modularity, revealing the effects of different combinations of the resolution and inter-layer coupling functions. This work also represents a starting point for the development of new optimization methods for community detection in multilayer networks. 






 


 



 




 

§ INTRODUCTION
 
Complex network systems, such as social networks, biological networks, and transportation networks, are inherently organized into communities, a.k.a. clusters or modules,  with dense internal links and sparse external links.  
  Since members of a community tend to generally share common
properties, revealing the community structure in a network can provide a better understanding of the overall functioning of the network. 

 

 The well-known    modularity <cit.> function was originally conceived to evaluate  a community structure in a network graph in terms of difference  between the  actual   number of edges linking nodes inside each community and the expected connectivity in the null model. Typically, the expected connectivity is expressed through a configuration graph model, having a certain degree distribution and randomized edges. Since this graph ignores any community structure, a large difference between actual connectivity and expected connectivity would indicate   the presence of a  community structure.   

Modularity has been widely utilized as objective function in several optimization methods designed for discovering communities  in networks <cit.>, including greedy agglomeration <cit.>,  spectral division <cit.>, simulated annealing <cit.>, 
  or extremal optimization <cit.>.
  
 
 
Traditional  approach to network analysis refers to the modeling of a real-world  system  as a single network of interacting entities. While this approach has been 
 
widely used to study a variety of applications, 

 there are plenty of scenarios for which this methodology appears strongly limiting <cit.>. 
In general, ties among entities could be induced by one or several types of  relations or interactions, or even be dependent on side-information 
based on specific dimensions or aspects of interest for the entities in the network. Within this view, multilayer networks
 <cit.> represent a powerful tool
to model systems interconnected by 
 multiple types of relations. In the multilayer network model, each type of connection is represented by a layer, and an entity may be present in 
different layers based on the type of relations it is connected to its neighbor  entities. 
 Just to mention one real example, nowadays online users usually have multiple accounts across different online social networks, and several online/offline relationships are likely to occur among the same group of individuals (e.g., following relations, like/comment interactions, working together, having lunch).
 
 It should be emphasized that  neglecting such a kind of complex 
organization by reducing the whole system to a single network (e.g., through some kind of projection,  or by aggregation), has been shown to be much less 
informative than the multilayer representation <cit.>.  
For the above mentioned reasons, multilayer networks are experiencing an increasing interest from the scientific community, leading to an explosion of 
scientific papers in many areas of science, thus becoming one of the most used tools for interdisciplinary research <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>.
 
Clearly, the problem of community detection in such multilayer networks takes a central role to unveil meaningful patterns of node groupings into communities,  by leveraging the various interaction modes that involve all the entity nodes in the network.   
To address these needs, modularity has been  extended to the   general case of  multilayer networks. In particular, Mucha et al. <cit.> extend the modularity function to   arbitrary multilayer networks (also called multislice in that work), by introducing two additional parameters w.r.t. classic modularity: a resolution parameter  and an inter-layer coupling factor. 
 The resolution parameter acts on the expected connectivity terms, thus controlling the effect on the size distribution of community due to the resolution limit known in modularity <cit.>. 
 The inter-layer coupling factor focuses on the links across layers and hence impacts on the strength of the inter-layer connections of entities in the network. 
 While being important to enhance the ability of modularity in evaluating a community structure, the two parameters introduced in the multislice modularity are nonetheless subjected to   arbitrary choices, which raise a number of issues in the application of this modularity function.  In particular,  the resolution parameter can be  arbitrarily   set for each layer, but it discards any structure information at graph or community level. Moreover, the inter-layer coupling terms do not differentiate among the selected layers, and  all pairs of layers can in principle be considered, which  makes no sense in certain scenarios such as modeling of time-evolving networks.  
 
 The above considerations prompted us to revisit the notion of modularity in multilayer networks, and in particular to introduce  novel aspects to take into account in both the resolution and inter-layer coupling definitions.  First, the layer-specific resolution factor is also made dependent on  each particular community. We notice that, since a high-quality multilayer community should embed high information content among its nodes,  the resolution of a specific layer to control the expected connectivity of a given community in the modularity function should be lower as the contribution of that layer to the information content of the community is higher. 
 By relating the information content of a multilayer community to the amount and  variety of types of links   internal to the community, we  provide a new definition of resolution factor based on the concept of redundancy of community. 
 Second, to determine the strength of coupling of nodes across layers, 
 we again  consider it at community level, such that for each pair of layers, the inter-layer coupling factor for nodes in a community depends on the relevance of the community projection on the two layers.    
 Moreover, we account for an available ordering of the layers, and relating constraints on their coupling validity. 
 

 Our main contributions are summarized as follows:


  * We propose a novel definition of multilayer modularity, in which we reconsider  the role and semantics of its two key terms, that is, the  resolution factor and inter-layer coupling factor. We conceive parameter-free unsupervised approaches for their computation, which leverage information from the within-layer and across-layer structures of the   communities in the multilayer network. 
Moreover, our formulation of multilayer modularity is general enough to account for an available ordering of the layers, therefore is also well-suited to deal with temporal multilayer networks. 

  * We provide  theoretical insights into properties of the proposed multilayer modularity. More specifically, we investigate the effect of increasing the number of communities in the behavior of the multilayer modularity, and we analytically derive the lower and upper bounds in the values of the multilayer modularity. 

  * We conduct an extensive experimental evaluation, primarily  to understand   how the proposed multilayer modularity behaves w.r.t. different settings regarding the resolution   and the inter-layer coupling terms.  Using 4  state-of-the-art methods for multilayer community detection (GL, PMM, LART,   and ), LFR synthetic multilayer networks and 10 real-world multilayer networks, 
results have shown the significance of our formulation and the different  expressiveness against the previously existing multislice    modularity.  


 

  

 

  
  


§ RELATED WORK
   

Community detection is a key-enabling task in network analysis and mining, with tons of methods developed in the last  ten years — please refer to <cit.> for surveys on this topic.  
 
In addition, 
   different metrics for community structure evaluation have been introduced. 
As discussed in Section <ref>, the most popular and widely accepted  measure is the so-called “modularity”, defined by Newman <cit.>.  Initially conceived for undirected networks, the modularity function has been subsequently extended to cover different cases. 
In <cit.><cit.>, modularity has been generalized to directed networks incorporating information contained in edge directions, while in <cit.>   modularity is also adapted to capture communities in weighted networks.  
To overcome the well-known modularity resolution limit <cit.>, 
in <cit.> and <cit.> modularity has been modified by incorporating a resolution parameter that helps reveal communities at different resolution scale.  
 A further step towards a generalization of the modularity refers to its extension to  signed networks <cit.>.  
 
  Also, to deal with bipartite networks,  modifications have been proposed in <cit.>. 
 
To uncover overlapping communities, in <cit.> the authors propose an extension to the modularity function that includes the notion of belonging (or membership) coefficient, which measures to which extent a node belongs to the various  communities. This approach is sometimes referred to   as fuzzy community discovering.   
  Finally, as introduced in Section <ref>, modularity has been generalized by <cit.> to capture communities in multislice  networks. 
  Such a version of the modularity function is detailed in the next section.  
 
 


§ BACKGROUND
 




 §.§ Modularity


Given an undirected graph G=(,ℰ), with n=|| nodes and m=|ℰ| edges,  let 𝒞 be a community structure over G. For any  v ∈,  we use   d(v) to denote the degree of v, and for any community C ∈𝒞, symbol d(C) to denote the degree of   C, i.e.,  d(C)=∑_v ∈ C d(v); also,  the  total degree of nodes over the entire graph,     d(), is defined as d()=∑_v ∈𝒱 d(v)=∑_C ∈𝒞 d(C)=2m.  
Moreover, we denote with d_int(C)  the internal degree of C, i.e., the portion of d(C) that corresponds to the number of edges linking nodes in C to other nodes in C. 
Newman and Girvan's modularity is defined as follows <cit.>:



    Q_NG(𝒞)=  ∑_C ∈𝒞d_int(C)/d() - ( d(C)/d())^2

 
 In the above equation, the first term is maximized when many edges are  contained in clusters, whereas the second term is minimized   by partitioning  the graph into many clusters with small total degrees. 
The value of Q_NG ranges within  -0.5 and 1.0 <cit.>;  
it is  minimum for any bipartite network with canonic clustering, and  maximum when the network is composed by disjoint cliques. 




 §.§ Multilayer network model

Let G_ℒ = (V_ℒ, E_ℒ, , ℒ) be a multilayer network graph, where  is a set of entities  and  ℒ= {L_1, …, L_ℓ} is a set of layers. Each layer represents a specific type of relation between entity nodes. 


Let V_ℒ⊆×ℒ be the set containing the entity-layer combinations, i.e., the occurrences of each   entity  in the   layers.   E_ℒ⊆ V_ℒ× V_ℒ is the set of undirected links connecting the entity-layer elements. 
 For every  L_i ∈ℒ, we define V_i = {v ∈ |     (v, L_i) ∈ V_ℒ}⊆𝒱 as the set of nodes in the graph of L_i, and E_i ⊆ V_i × V_i as the set of edges in L_i. Each entity must be present in at least one layer, i.e., ⋃_i = 1..ℓ V_i  =, but each  layer is not required to contain all  elements of . 
We assume that the inter-layer links only connect the same entity in different layers, however  each entity in one layer could be  linked to the same entity in a few or all other layers.  




 §.§ Multislice Modularity


Given a community structure   𝒞 identified over a multilayer network G_ℒ, the multislice modularity <cit.> of   𝒞 is defined as:
 
 

    Q_ms(𝒞)     =    1/d(V_ℒ)∑_u,v, 
     L_i,L_j[  ( A_uvL_i - γ_i d_L_i(u)d_L_i(v)/2|E_i|)  ×. 
       ×   . δ_L_i,L_j +  δ_u,vC_v,L_i,L_j] δ(g_u,L_i, g_v,L_j)



where 
d(V_ℒ) is the total degree of the multilayer network graph, 
d_L_i(u) denotes the degree of node u in layer L_i, 
A_uvL_i denotes a link between u and v in L_i,  

   2|E_i|  is the total degree of the graph of layer L_i, γ_i is the resolution parameter for layer L_i,   C_v,L_i,L_j quantifies the links of node v across  layers L_i,  L_j.  
Moreover, the Dirichlet terms have the following meanings: 
δ_L_i,L_j is equal to 1 if L_i=L_j and 0 otherwise, 

δ_u,v is equal to 1 if u=v and 0 otherwise (i.e., the inter-layer coupling is allowed only for nodes corresponding to the same entity), 
and 
δ(g_u,L_i, g_v,L_j) is equal to 1 if the community assignments of node u in L_i and node v in L_j are the same and 0 otherwise.  

 

Limitations of Q_ms.  
As previously  mentioned, a different   resolution parameter   γ_i can be associated with each layer to express its relevance weight; however,  in <cit.>, there is no  specification of  any principled way to set a layer-weighting scheme, possibly including information from the available multilayer community structure.  
Moreover, neither the inter-layer coupling term  (i.e., C_v,L_i,L_j) or any constraint on the layer comparability are clearly  defined; actually,  all nonzero inter-layer edges are set to a constant value ω, for all unordered pairs of layers.    
In general, both γ_i  and ω parameters can   assume any non-negative value, which further increases a clarity issue in the     properties of Q_ms. 
    


 


§ PROPOSED MULTILAYER MODULARITY

In this section, we propose a new definition  of   modularity for multilayer networks that aims to overcome all of the issues of Q_ms previously discussed. We pursue this goal by focusing on the role and semantics of the two key elements in multilayer modularity: the layer-specific resolution   and the inter-layer coupling.  

Our definitions of the two terms are independent on a-priori assumptions on the network and/or user-specified settings; by contrast, we  conceive  parameter-free unsupervised approaches for their computation, by leveraging information from the within-layer and inter-layer structures of the communities.  

 Our proposed   resolution factor is computed for pairs of layer and community, rather than  for each layer globally. Analogously, to define the inter-layer coupling term,  we   account for properties related to a community 
 on two layers; more in detail, 
 we evaluate the projections of a community  over any two comparable layers, i.e., the sets of nodes belonging to a community that lay on those layers.  
Remarkably, the comparability of layers is another key aspect of our definition of modularity:  we generalize   the inter-layer coupling term by admitting the existence of a partial order relation ≺_ℒ over the layers, in order  to properly represent  scenarios in which a particular ordering among layers is required. For instance, it may be the case that the network layers 
 have to be processed according to their natural order (e.g., lexicographic order of the network labels), or according to a temporal order; moreover, it may be required to compare   adjacent layers only, or each layer with any other  succeeding it in the ordering.  Figure <ref> provides an illustration of a multilayer network and the aforementioned key aspects we deal with in our proposed multilayer modularity, which is formally presented next.  

 




 


Let G_ℒ = (V_ℒ, E_ℒ, , ℒ) be  a multilayer network graph, and let  ≺_ℒ be  an optionally provided  partial order relation over the set of layers ℒ. 
Given a community structure 𝒞={C_1, …, C_k} as a partitioning of the multilayer graph G_ℒ, the   multilayer modularity  is defined as: 



    Q(𝒞)      =   1/d(V_ℒ)∑_C ∈𝒞[ ∑_L ∈ℒ(  d_L^int(C)  - γ(L,C) (d_L(C))^2/d(V_ℒ)   +    β∑_L' ∈𝒫(L)  IC(C, L, L') ) ]
 

where for any C ∈𝒞 and L ∈ℒ:


  * d_L(C) and d_L^int(C) are  the degree of  C and the internal degree of  C, respectively, by considering only edges  of layer L; 

  * γ(L,C) is the value of the  resolution function;  

 
  * IC(C, L, L') is the value of the inter-layer coupling function for any valid layer pairings with L;  
 
  * β∈{0,1} is a parameter to control the  exclusion /inclusion of inter-layer couplings; 
  and 

  * 𝒫(L) is  the set of  valid pairings   with   L defined as:

    𝒫(L)=
    {L' ∈ℒ |  L ≺_ℒ L' },    if≺_ℒis defined
    ℒ∖{L},    otherwise
 


 

 


 


Notably, unlike  the multislice modularity in  Eq. (<ref>), our proposed modularity originally introduces a resolution factor that varies with each community, and an inter-layer coupling scheme that might also depend on the layer ordering. Moreover,  Eq. (<ref>)  utilizes  the total degree d(V_ℒ) of the multilayer   graph  instead of the layer-specific degree (i.e., term  2|E_L_i|, for each L_i ∈ℒ). 
 The latter difference w.r.t. the multislice modularity  is also important because, as we shall later discuss more in detail,   the total degree of the multilayer graph includes the inter-layer couplings and it might be defined in different ways depending on the scheme of inter-layer coupling.  

 
In the following, we elaborate on the resolution functional term, γ(·,·), and the inter-layer coupling functional term,  IC(·, ·, ·). 









 §.§ Redundancy-based resolution factor
 
The layer-specific resolution factor intuitively expresses the  relevance of a particular layer to the calculation of the expected community connectivity in that layer. 
While this   can always reflect some predetermined scheme of relevance weighting of layers, we propose a more general definition that  accounts for the strength of the contribution that   a layer takes  in determining the internal connectivity for each community.  
The key assumption underlying our approach is   that, since a high quality community should envelope   high information content among its elements, the resolution   of a layer to control the expected connectivity of    a given community  should be lowered  as the layer's  contribution to the information content of the community increases. 

In this regard, the redundancy measure proposed in <cit.>  is particularly suited to  quantify  the variety of connections, such that it is higher if edges of more types (layers) connect each pair of nodes in the community.  
 
Let us  denote with P_1 the set of node pairs   connected in at least one layer in the graph, and with  P_2 the set of “redundant” pairs, i.e., the pairs of nodes connected in at least two layers. 
Given a community C, P_1^C and P_2^C denote the subset of P_1 and the subset of P_2, respectively, corresponding to nodes in   C. 
The redundancy of C, ρ(C),  expresses  the number of pairs in C with redundant connections, divided by  the  number of layers connecting the pairs. Formally <cit.>:  
 



    ρ(C) =   ∑_(v,u) ∈ P_2^CSL(v,u,ℒ)/|ℒ| × |P_1^C|,
 

with SL(v,u,ℒ)={L ∈ℒ |   (v,u,L) ∈ E_ℒ}.



Note that in the above formula, each of the sets SL refers to the layers on which two nodes in a redundant pair are linked. Upon this concept,  we define the set of supporting layers sup for each community C as:


    sup(C,ℒ)=⋃_(v,u) ∈ P_2^C  SL(v,u,ℒ).
 




Using the above defined sup(C,ℒ), we provide the following definition of redundancy-based resolution factor.  











Given a layer L and a community C,   the  redundancy-based resolution factor in Eq. (<ref>) is defined as:


    γ(L,C) = 2/1+ log_2(1+ nrp(L,C))
 

where 

nrp(L,C) = |{ s=SL(v,u,ℒ) ∈ sup(C,ℒ)   |  L ∈ s }|

expresses the number of times  layer L participates in redundant pairs.













Note that γ(L,C) ranges in (0, 1] ∪ [2]. In particular, it ranges   in (0, 1]  as long as L participates in at least one redundant pair, and it decreases as nrp(L,C) increases; moreover, as special case,  γ(L,C) is equal to 2  when nrp(L,C) = 0.   









Consider the network with  16 entities and 5 layers  shown in Fig. <ref>, and let us first focus on some specific cases for the computation of the γ(L,C) terms.    
  For instance, given community C_4 and  layers L_1 and L_5,  
  the corresponding values of redundancy-based resolution  are equal to 2, because L_1, L_5 never participate in redundant pairs for nodes of C_4. 
   By contrast,   layers L_2 and L_5 participate in one redundant pair for community C_1, which corresponds to the edge linking entities 2 and 3 for L_5 and the edge linking entities 1 and 2 for L_2; therefore,  the   values of redundancy-based resolution associated with C_1 on L_2 and L_5 are equal to 1. Also, the resolution for C_1 on L_1 takes a value lower than 1 since there  is  more than one redundant pair. 
   Table <ref> reports on the entire set of values for the resolution factor computed on the network of Fig. <ref>. 










 §.§ Projection-based inter-layer couplings

We propose a general and versatile approach  to quantify the strength of coupling of nodes in one layer with nodes on another layer. Our key idea is to determine the fraction of nodes belonging to a community onto a layer that  appears in the projection of the community on another layer, and express the relevance of this projection w.r.t. that pair of layers.  

Given a community C ∈𝒞 and layers L_i, L_j ∈ℒ, we will use symbols C^(i) and C^(j) to denote the projection of C onto the two layers, i.e., the set of nodes in C that lay on L_i and L_j, respectively. 

In the following, we define   two   approaches for measuring inter-layer couplings based on community projection. 

For any two layers L_i, L_j and community C, the first approach, we call symmetric, determines the relevance of inter-layer coupling of nodes belonging to C as proportional to the fraction of nodes shared between L_i and L_j that belong to C. 




Given a community C ∈𝒞 and layers L_i, L_j ∈ℒ,  the symmetric projection-based inter-layer coupling, denoted as IC_s(C,L_i, L_j) and referring to term IC  in Eq. (<ref>), is defined as the probability that C lays on L_i and L_j:


    IC_s(C,L_i, L_j)=[C  in  L_i, C  in  L_j] = |C^(i)∩ C^(j)|/|V_i ∩ V_j|




The above definition assumes that the two events “C in L_i” and  “C in L_j” are independent to each other, and it   does not consider that  the coupling might have a different meaning depending on the relevance a  community has on a particular layer in which it is located.  
By relevance of community, we simply mean here the fraction of nodes in a layer graph that belong to the community; therefore, the larger the community in a layer, the more relevant is w.r.t. that layer. 
However, we observe that more relevant community in a layer corresponds to less surprising projection from that layer to another. This would imply that the inter-layer  coupling for that community is less interesting w.r.t. projections of smaller communities, and hence the strength of the coupling might be lowered.  

We capture the above intuition by the following definition of   asymmetric projection-based inter-layer coupling. 

 
 
 
 
Given a community C ∈𝒞 and layers L_i, L_j ∈ℒ,  the asymmetric projection-based inter-layer coupling, denoted as IC_a(C,L_i, L_j) and referring to term IC  in Eq. (<ref>), is defined as the probability that C lays on L_j given that C lays on L_i:

 
    IC_a(C,L_i, L_j)     =    [C  in  L_j | C  in  L_i]  =  
        =    [C  in  L_i, C  in  L_j]/[C  in  L_i]
    =  
        =    |C^(i)∩ C^(j)|/|V_i ∩ V_j|×|V_i|/|C^(i)|

 

  



Dealing with layer ordering.  

Our formulation of multilayer modularity is general enough to account for an available ordering of the layers, according to a given partial order relation. 

The previously defined asymmetric inter-layer coupling is well suited to model  situations in which we might want to express the inter-layer coupling from a “source” layer to a “destination” layer. Given any two layers L_i, L_j, it may be the case that only comparison of  L_i to L_j, or vice versa, is allowed. 
This is clearly motivated when    there exist layer-coupling constraints, thus only some  of the layer couplings should  be considered in the computation of Q.  

In practical cases, we might assume that the layers can be naturally ordered to reflect a predetermined lexicographic order, which might be set, for instance, according to a progressive enumeration of layers or to a chronological order of the time-steps corresponding 
to the layers.  
That said,   we can consider two special cases of layer ordering:



  * Adjacent layer coupling: 
L_i ≺_ℒ L_j iff  j=i+1 according to a predetermined natural order.  

  * Succeeding-layer coupling:   L_i ≺_ℒ L_j iff  j>i according to a predetermined natural order. 




Note that the adjacent layer coupling scheme requires ℓ - 1 pairs to consider, while the succeeding-layer coupling scheme involves the comparison between a layer and its subsequent ones, i.e., (ℓ^2-ℓ)/2 pairs. 



  Moreover, it should be noted that the availability of a layer ordering enables two variants of the asymmetric projection-based inter-coupling given in Eq. (<ref>).  
For any two layers L_i, L_j ∈ℒ, such that L_i ≺_ℒ L_j holds, we refer to as  inner the direct evaluation of IC_a(C, L_i, L_j), and as outer  the case in which L_i and L_j are switched, i.e., IC_a(C, L_j, L_i). 




In the inner case, the strength of coupling is higher as the projection of C on the source layer (i.e., the   preceding one  in the order) is less relevant;   vice versa,  the outer case 
weights more the coupling as the projection on the destination layer  (i.e., the   subsequent one in the order) is less relevant.




Consider again the example network of Fig. <ref>. The asymmetric coupling for the projection of community C_2 from L_2 to L_3 is IC_a(C_2, L_2, L_3) = (2/5) × (9/4)=9/10=0.9 in the inner case, and  IC_a(C_2, L_3, L_2)= (2/5) × (6/3)=4/5=0.8 in the outer case. 



We hereinafter use symbols IC_ia and IC_oa to distinguish between the inner asymmetric and the outer asymmetric evaluation cases. 
 





Time-evolving multilayer networks.  
So far we have assumed that when comparing  any two layers L_i, L_j, with L_i ≺_ℒ L_j,  it does not matter the number of layers between L_i and L_j. 
Intuitively, we might want to penalize the strength of the  coupling  as more “distant” L_j is from L_i. This is often the case in time-sliced networks, whereby we want to understand how community structures evolve over time. 

In light of the above remarks, we define a refinement of the asymmetric projection-based inter-layer coupling, by introducing a multiplicative factor that smoothly decreases the value of the IC_a function as the temporal distance between L_i and L_j increases. 




Given a community C ∈𝒞 and layers L_i, L_j ∈ℒ, such that  L_i ≺_ℒ L_j, the time-aware asymmetric projection-based inter-layer coupling, denoted as IC_a^t(C,L_i, L_j), is defined as  

 
    IC_a^t(C,L_i, L_j)=IC_a(C,L_i, L_j) ×2/1+ log_2(1+j-i)

 
 
 
 Note that the second term in the above equation is 1 for the adjacent layer coupling scheme, thus making no penalization effect when only consecutive   layers are considered. 




 





 
Referring again to the example in Fig. <ref>, 
in Table <ref> we summarize the mean and standard deviation values for the different variants of the inter-layer coupling factor. 
 
 One remark is  that the values for communities C_2 and C_3 are higher than those corresponding to communities C_1 and C_4. 

This is mainly due to the representativity of  C_2 and C_3 in all layers.  The lowest values are obtained for community C_4, which is in fact less represented than other communities (only in layers L_1 and L_5). For instance, let us focus on this community. The mean inter-layer coupling factor IC_oa^Suc for community C_4 is 0.11, since: |C_4^(5)∩ C_4^(1)|=2, |V_5 ∩ V_1|=10 (which is exactly the size of L_5 without node 16), |V_5|=11 and |C_4^(5)|=2; this  determines an inter-layer coupling factor of 1.1, which is divided by the admissible pairings of layers, i.e., 10. On the contrary, the mean IC_oa^Adj for C_4 is equal to zero, because the projection of this community is always empty when the adjacent coupling scheme is used. 
 
Finally, Table <ref> reports the multilayer modularity values, including the com­munity-specific contributions. Note that, regardless of the settings of γ and IC factors,   

 communities C_2 and C_3 obtain the highest  values of modularity, which is mainly determined  since   they are disconnected from the rest of the graph at layer L_3. 
 In general, it should be noted that the contribution given by each community is consistent w.r.t. the various settings of γ and IC factors.  
 Also, it is interesting to note that discarding the inter-layer couplings (IC=0, which corresponds to the 9th and 13th columns) can lead to   
  values of community-modularity and global modularity that  tend  to be much higher than the corresponding cases with IC≠ 0.  This overestimation can also occur, though to a lesser extent, when fixing γ=1 (13th column) vs. redundancy-based γ (9th column). 
Also, it is worth noting that using the redundancy-based resolution factor γ with unordered layers (2nd, 3rd and 4th columns) increases the community as well as global modularity  vs. the same cases with γ=1 (10th, 11th, and 12th columns).





 







  §.§.§ Relations between the resolution and inter-layer coupling factors

 


Both factors take into consideration the network context, however they differ in that  
γ(L, C) considers a “global” multiplex context, whereas IC considers a “local” multiplex context. 
Intuitively, γ(L, C)  is defined for each valid   layer-community pair according to the status of the links among nodes in   community C that lays on L versus their status on the other layers.  By contrast,   IC considers the status of the same community from one layer to  another comparable layer.  

In terms of numerical comparison,  when the size of the community structure  tends to the number of nodes of the network,  γ(L, C) tends to increase (i.e., to the maximum value of 2) while   IC tends to decrease (i.e., to zero).



 

 §.§ Properties of the proposed  multilayer modularity

We provide theoretical insights on   Q, focusing on the effect of increase in the size of the community structure and on the analytical derivation of the range of values of Q.  

 


  §.§.§ Effect of increase in the number of communities

 



We discuss the effect of increasing k (i.e., decreasing the average size of communities in 𝒞) by distinguishing three configurations of Q: (i) symmetric inter-layer coupling, (ii) asymmetric inter-layer coupling, and (iii) ordered layers.

In the first case, Q tends to have a monotonic decreasing trend.   
  
 This is easily explained by the combination of three contingencies. The first one is an average decrease in the internal degree d_L^int. The second contingency is an   increase in the redundancy-based resolution factor γ: in fact, smaller communities correspond to lower probability of observing redundant pairs within  communities over   different layers; this decreases the logarithmic term in the resolution factor,   which will progressively  tend to 2 (maximum value). The third contingency is a decrease in the inter-layer coupling factor IC_s, since the size of community intersection  becomes increasingly smaller as the community size decreases.
 
By contrast, when equipped with the asymmetric projection-based inter-layer coupling IC_a, Q tends to  differ from  a monotonic decreasing trend because of the bias term |V_i|/|C^(i)|, which increases with communities of smaller size. 


In the third case (i.e., ordered layers), Q can again follow an increasing or decreasing trend. Recall that  the term d(V_ℒ)    includes the contribution of the inter-layer edges, which obviously are fewer   when the  layer couplings  are order-dependent. A decrease in the number of inter-layer couplings    also makes the decrease in the actual connectivity term (i.e., d_L^int(C)/d(V_ℒ))  slower as k increases, since d(V_ℒ) is smaller than in the unordered-layer  contingency. 
Consequently, the inter-layer coupling term could compensate the actual connectivity term, which will result in increasing the value of Q. 
Finally, considering time-aware asymmetric  inter-layer coupling, Q is more likely to  follow a  decreasing  trend because of the effect due to the smoothing term 2/1+ log_2(1+j-i), which penalizes IC_a for any two  no time-consecutive layers. Consequently, since  the inter-layer coupling factor IC_a^t is smaller than   IC_a, Q could  monotonically decrease despite  the bias term |V_i|/|C^(i)| in IC_a^t. 






  §.§.§ Lower and upper bounds




To determine the range of values of the basic modularity in simple graphs, the theoretical frameworks   previously studied in <cit.> and <cit.> define two canonical structures to support the analytical computation of the minimum and maximum value of the modularity, respectively. More specifically, the former work proved that  any bipartite graph  with the canonical two-way clustering   obtains the minimum value of modularity, whereas the latter work proved that the maximum modularity is reached in a graph composed of  disjoint  cliques. 

 
 
 Following the lead of the above works, here  we provide  theoretical results about   the analytical derivation of the lower bound and upper bound   of our proposed multilayer modularity.   
 
 
 
 Given a multilayer network G_ℒ = (V_ℒ, E_ℒ, , ℒ), with n=||, ℓ = |ℒ|, and a community structure 𝒞 for G_ℒ, 
 the lower bound of Q  is as follows:   


    Q(𝒞) =   -n^2ℓ/(nℓ +2p)^2 + 4(1+η)p/n^2(nℓ + 2p),


with η=0 for IC_s and η=1 for IC_a, 
and  p= ∑_L ∈ℒ|𝒫(L)| is   the total number of valid layer-pairings.
 
 
 Proof.    Proof is reported in the Appendix. 


 
 Given a multilayer network G_ℒ = (V_ℒ, E_ℒ, , ℒ), with n=||, ℓ = |ℒ|, and a community structure 𝒞 for G_ℒ, 
 the upper bound of Q  is as follows:   


    Q(𝒞)  
     =     2 [ 1/2(n/2-1)ℓ/(n/2-1)ℓ+p-γℓ(1/2(n/2-1)/(n/2-1)ℓ+p)^2 + . 
    
     +    .  (1+η)p/n^2(n/2-1)ℓ+n^2p]


with η=0 for IC_s and η=1 for IC_a, 
and  p= ∑_L ∈ℒ|𝒫(L)| is   the total number of valid layer-pairings.

 
 Proof.    Proof is reported in the Appendix. 

  

Note that,  in the special case for β=0, i.e.,  the inter-layer coupling factor is discarded, the lower bound of Q is 


    Q(𝒞) =-1/ℓ.


Analogously, the upper bound of Q is rewritten as:
 
 
    Q(𝒞) =   2ℓ-γ/2ℓ,


with γ=2(1+log_2(1+n/2(n/2-1)/2))^-1.
 









 


§ EVALUATION METHODOLOGY



We discuss here the evaluation networks (Sect. <ref>), the multilayer community detection methods  (Sect. <ref>),  and the experimental settings (Sect. <ref>).
 
 


 §.§ Datasets

 

 
Our selection of network datasets was motivated to fulfill the reproducibility requirement: in fact, all of our evaluation datasets, including both real-world networks and synthetic generators, are publicly available. Moreover, we also took the opportunity of diversifying the choice of real-world networks by considering various domains that are profitably modeled as multilayer networks. 
 

 

  §.§.§ Real-world network datasets

We considered 10 real-world multilayer network datasets. 




AUCS  <cit.>  
describes relationships among university employees:   work together, lunch together, off-line friendship, friendship  on Facebook, and coauthorship. 






EU-Air transport  network <cit.>  
(EU-Air, for short) 
represents European airport connections considering  different  airlines. 


FAO Trade network (FAO-Trade) <cit.> represents different types of trade relationships among countries, obtained from FAO (Food and Agriculture Organization of the United Nations).  


 (stands for FriendFeed, Twitter, and YouTube) <cit.>  was built by exploiting the feature of FriendFeed as social media aggregator to align registered users who were also members of Twitter and YouTube. 


Flickr refers to the  dataset studied in <cit.>. 

We used the corresponding    timestamped interaction network whose links express “who puts a favorite-marking to a photo of whom”.   
We extracted the layers on a month-basis and aggregated every six (or more)  months. 


 (stands for GitHub, StackOverflow and Twitter) <cit.> is another cross-platform network   
where  edges express followships on Twitter and GitHub, and  “who answers to whom” relations on   StackOverflow.


Higgs-Twitter <cit.>  
represents  friendship,   reply,   mention, and retweet relations among Twitter users. 


London   transport network <cit.>  
 (London, for short) models three types of connections of train stations in London: underground lines, overground, and DLR.   


ObamaInIsrael2013 <cit.>
(Obama, for short) models retweet, mention, and reply relations of users of  Twitter during   Obama's visit to Israel in 2013.


7thGraders  <cit.>  
(VC-Graders, for short) represents  students involved in    friendship,  work together, and affinity relations in the class. 



  Table <ref> reports for each dataset, the size of set , the number of edges in all layers, 
 and the average coverage of node set (i.e.,  1/|ℒ| ∑_L_i ∈ℒ (|V_i|/|𝒱|)). 
The table also shows basic, monoplex structural statistics (degree, average path length, and clustering coefficient)  for the layer graphs of each dataset.   
    
    


 
 
 
 
 
 
 
 

  §.§.§ Synthetic network datasets

 
Besides the real-world network data, we generated four synthetic multilayer network datasets. Our goal was the evaluation of the multilayer modularity Q on different network models.
 
 Two out of the four networks are composed of 2 layers and 256 entities. In one   network, hereinafter referred to as ER-ER,  the two layers are Erdös-Rényi (ER) random graphs. In the second network, dubbed LFR-ER, the first layer is generated by the Lancichinetti-Fortunato-Radicchi (LFR) benchmark, while the second layer is an Erdös-Rényi random graph. The other two networks are composed of 4 layers and 128 nodes. Both networks are characterized by two Erdös-Rényi  layers  and two layers built as Girvan-Newman (GN) graphs, but they differ in the layer ordering:  GN-ER-GN-ER in the first network, and GN-ER-ER-GN in the second network.  

Moreover, mainly for purposes of efficiency evaluation, 

 we  generated a set of synthetic multilayer networks using the Lancichinetti-Fortunato-Radicchi (LFR) benchmark. In particular, single-layer network datasets were provided by the LFR benchmark using a variable number of nodes with steps of 128 until 1024. Also, the maximum and average node degrees were set to 16, and the mixing coefficient μ was set to 0.1.   Each network dataset was characterized by four communities.  
From each of such networks, a multilayer network   was created by replicating the LFR single-layer     from 2 to 10. 





 §.§ Community detection methods


 
We resorted to   state-of-the-art methods for  community  detection in multilayer networks, which belong to the two  major  approaches, namely  aggregation and direct methods. The former  detect a community structure separately for each network layer, after that  an aggregation mechanism is used to obtain the final community structure, while the latter directly work on the multilayer graph by    optimizing a multilayer quality-assessment criterion. 
(Note that while it is possible to flatten the multilayer graph in order to apply on it any conventional  community detection algorithm, this approach can be too simplistic, since, e.g.,  it would not permit to investigate about the temporal evolution of communities.) 

 
As   exemplary methods of the   aggregation approach, we used  Principal Modularity Maximization (PMM) <cit.> and Enhanced Modularity-driven Ensemble-based Multilayer Community Detection () <cit.>. 
 PMM aims to find a concise representation of features from the various layers (dimensions) through   
 structural feature extraction and cross-dimension integration.  
Features from each dimension are first extracted via modularity maximization, then concatenated and subjected to PCA to select the top eigenvectors, which represent possible community partitions. Using these eigenvectors, a low-dimensional embedding is computed  to   capture  the principal patterns across all the dimensions of the network, finally a   k-means on this embedding is carried out to discover a community structure.  

 
 is a parameter-free extension of  the method proposed in <cit.>. Given an ensemble of community structures available for a multilayer network, optimizes a consensus objective function  to discover a consensus solution with maximum modularity, subject   to the constraint of being searched over a hypothetical space of consensus community structures that are valid w.r.t. the input ensemble and topologically bounded by two baseline solutions.  To detect the initial cluster memberships of nodes,  utilizes a   consensus or co-association matrix, which stores the fraction of clusterings in which any two nodes are assigned to the same cluster. To filter out noisy, irrelevant co-association, a user-specified threshold must be specified.  
Besides introducing flexibility in community assignments of nodes during the modularity optimization,  overcomes the limitation of setting such a parameter of minimum co-association,   by providing a parameter-free identification of consensus clusters   based on generative models for graph pruning. 

 
As for the direct methods, we resorted to Generalized Louvain (GL) <cit.> and Locally Adaptive Random Transitions (LART) <cit.>. 

GL extends the classic Louvain method using  multislice  modularity, so that each node-layer tuple is assigned separately to a community.  Majority voting is adopted to decide the final assignment of an entity node to the  community that contains the majority of its layer-specific instances.  

LART is a   random-walk based method. It first runs a different random walk for each layer, then a dissimilarity measure between nodes is obtained leveraging the per-layer transition probabilities. A hierarchical clustering method is used   to produce a hierarchy of  communities which is eventually cut at the level corresponding to the best value of multislice modularity.  
 

It should be emphasized that  we selected the above methods because, while having  different characteristics, they all use modularity either as optimization criterion (GL, PMM   and ) or as evaluation criterion to produce the final community structure (LART).  







Note also that  PMM requires the desired number of communities (k) as input. 
Due to different size of our evaluation datasets, we devised several configurations of variation of parameter k in PMM, by reasonably adapting each of the configuration range and increment step   to the network size.  Concerning , we used  the marginal likelihood filter (MLF) to perform parameter-free detection of the number of communities <cit.>. 
 
 

It should be noted that the selected methods actually discover different community structures, thus supporting our choice in terms of diversity of evaluation scenarios for  the two competing modularity measures under study. 
 Table <ref> reports the number of communities of the solutions found by the various methods on the  real-world network datasets. (The number of communities k in PMM is selected according to the solution with highest modularity value.)  
We found  that GL tends to discover a high number of communities for larger networks (i.e.,   Flickr, Higgs-Twitter,  ,    and Obama), and 
the size distribution   of these communities (results not shown)   is highly right-skewed  
on the larger networks, while it is moderately left-skewed on the remaining datasets.  A similar result can be observed in for the different networks, although in Higgs-Twitter and  (resp. )  the number of communities is much lower (resp. higher) than in GL. 
By contrast,  the best performances of PMM  usually correspond to a low and quite stable number of communities.  Also,   LART generally tends to produce much more communities than the other methods, on the networks for which it is able to discover communities. 

As a final general remark, we used the original implementations of the selected methods, based on the source code made available by the respective authors. We emphasize that it is beyond the goals of this work to make any performance improvement in the community detection methods under study, which hence are considered here with no intent of comparative evaluation and with all their limitations. (This justifies, in particular, the inability of LART in terminating the task for some network datasets.) 
 





 §.§ Experimental settings



We carried out GL, PMM, LART  and   methods on each of the network datasets and measured, for each community structure solution, our proposed multilayer modularity (Q) as well as the  Mucha et al.'s multislice modularity (Q_ms).  

 We   evaluated Q using the redundancy-based resolution factor γ(L,C) with either  the symmetric   (IC_s) or the asymmetric (IC_a) projection-based inter-layer coupling.   

 We  also considered   cases corresponding to   ordered   layers, using either the adjacent-layer scheme or the succeeding-layer scheme, and for both schemes considering inner (IC_ia) as well as outer   (IC_oa) asymmetric coupling.  
We further evaluated the case of temporal ordering, using   the time-aware asymmetric projection-based inter-layer coupling. 
Yet, we considered the particular setting  of uniform resolution (i.e., γ(L,C)=1, for each layer L and community C).  
 

As for Q_ms, we   devised two settings: the first by varying  
ω  within [0..2]  while fixing γ=1,  

 the second by varying γ and ω=1-γ <cit.>.   
 


 
 





§ RESULTS



We organize our main experimental results into two parts, depending on whether layer ordering  was considered in the evaluation networks. 
Experiments were  carried out on an Intel Core i7-3960X CPU @3.30GHz, 64GB RAM machine. 
 




 §.§ Evaluation with unordered layers
 


 




  §.§.§ Synthetic network datasets

 
Table <ref> reports the multilayer modularity Q, multislice modularity Q_ms and number of communities obtained by the GL solution on the four synthetic networks. 

 
One first  remark is that using the redundancy-based resolution factor γ always leads to higher Q w.r.t. the cases corresponding to γ fixed to 1.  In particular,  we observe gains up to 0.1 on ER-ER, 0.07 on LFR-ER, and 0.12 on GN-ER-GN-ER and GN-ER-ER-GN. 

Another  remark is that the fully combination of resolution and inter-layer coupling factors (i.e., rightmost two columns) tends to lower the value of Q w.r.t. the cases corresponding to varying γ with β=0 (i.e., third last column); moreover, the asymmetric inter-layer coupling results in a higher Q w.r.t. the symmetric setting of IC.  This would hint that when the normalization term in the Q equation accounts for  the inter-layer couplings, this results in lowering the value of  Q, which is turn smoother when the asymmetric setting is used.  
 

Comparing Q and Q_ms, it should be noted that the two measures behave consistently on ER-ER vs. LFR-ER, i.e., the presence of a layer with a (LFR) modular structure actually leads to an increase in both modularity measures w.r.t. ER-ER. 
By contrast,  Q tends increase faster  than Q_ms    on the two GN-ER networks: this can be explained since a higher number of layers  (as occurs for the two GN-ER networks than for the ER-ER and LFR-ER networks) has a higher effect on  the inter-layer coupling factor IC, which is not present in  Q_ms. 
 
 
  






  §.§.§ Real-world network datasets




 Tables <ref>–<ref>  and 
Fig. <ref> report Q  measurements on the community structure solutions obtained by the various  community detection methods. 

Concerning GL (Table <ref>), 
 we observe that with the exception of  on which effects on Q are equivalent,  
 using IC_a leads to higher  Q  than  IC_s.  
On average over all networks, using  IC_a yields an increment of  13.4% and 14.6% (with γ fixed to 1) w.r.t. the value of Q corresponding to IC_s.  
 This   higher performance of Q due to IC_a supports our initial hypothesis on the opportunity of   asymmetric inter-layer coupling.     
It is also interesting to note that, when  fixing γ to 1, Q decreases w.r.t. the setting  with  redundancy-based resolution γ(L,C)  — 
decrement of 11% and 12%   using IC_a and IC_s, respectively.   

 
 
 


 
Table <ref> shows results obtained from LART solutions. 

 (Due to memory-resource and efficiency issues shown by the currently available implementation of LART, we are able to report  results only on some networks). 
  
   We observe that  the relative performance difference  between IC_s and IC_a settings is consistent with results found in the GL evaluation; this difference is even extreme (0.98 or 0.99) on  EU-Air and London, which is likely due also to the different sizes of community structures detected by the two methods (cf. Sect. <ref>).  

 
 Table <ref> shows results obtained by solutions. Also in this case, 
 using IC_a generally leads to better   Q  than  IC_s, regardless of the setting of γ. In particular, the observed increase is higher  in VC-Graders and London  (0.2),  followed by AUCS (0.18) and Obama (0.13). Moreover, when  fixing γ to 1, in most cases Q decreases (0.01-0.06) w.r.t. the setting  with  redundancy-based resolution γ(L,C).
 
 


 
  
 
  


 


 Figure <ref> shows how Q varies in function of the number (k) of clusters given as input to   PMM. 
 
 One major remark is that  Q tends to decrease as k increases. This holds consistently   for the configuration of Q with symmetric inter-layer coupling; in fact, as discussed in Sect. <ref>, the decrease of Q for increasing k  depends on a combination of  decrease of the internal degree d_L^int, decrease of the symmetric inter-layer coupling factor IC_s, and     increase of the redundancy-based resolution factor γ(L, C).
 
 Moreover, values of Q corresponding to IC_a tend to be close to the ones obtained for IC_s on the large networks, while on the smaller ones, IC_a trends  are above IC_s, by diverging for high k in some cases; in particular,   in London 
modularity for IC_a follows a rapidly, roughly linear  increasing trend with k; even more evident is the divergence of the    IC_s and IC_a trends for AUCS.  Again, as we previously discussed in   Sect. <ref>, this is due to the bias term |V_i|/|C^(i)| of IC_a, which increases with communities of smaller size. 

 
Note that, from an inspection of  the behavior of Q for higher regimes of k, we also found that  Q values eventually tend to stabilize below 1.    

 As concerns  the   setting with γ fixed to 1 (results not shown), while the trends of Q for IC_a and for IC_s do not change significantly, the values are typically lower than those obtained with redundancy-based resolution, which is again consistent with results observed for GL, LART and evaluations.  
 
 

  


 

 

Correlation analysis.  
We investigated whether any correlation may exist at community-level between the value of Q   and selected   statistics based on   structural characteristics of the input network. 
For this purpose, we focused on the average path length, clustering coefficient, redundancy, and node- and edge-set coverage,  for each community in an evaluation network; note that the latter two statistics are computed as, given a community C, the fraction of nodes (resp. edges) in a layer L_i that belong to  C, averaged over all layers in the network.  




Figure <ref> shows the correlation between each of the above structural characteristics and the values of  Q, with redundancy-based resolution factor γ(L,C) and IC_s, on the solution found on selected networks by GL, and PMM; for the latter,   k was chosen as that corresponding to the best modularity performance. 
 Note also that  the correlation results obtained by Q with γ=1 and IC_s, γ=1 and IC_a, γ(L,C) only, and combination of  γ(L,C) and IC_a, do not show significant differences, hence their presentation is discarded.   
  
Looking at the three plots in the figure, we observe a mid-high positive correlation of Q with the topological measures  in most cases.  More in detail, in Fig. <ref> (a) the modularity of the solution found by GL on EU-Air shows an average correlation of 0.85 with the other measures. Also, an average correlation of 0.95 and 0.96 is obtained between Q  and respectively node-set  and edge-set coverage on . For AUCS, Q has a positive correlation of 0.76 with   clustering coefficient and a negative correlation with the other measures. For VC-Graders, Q shows a positive correlation with all measures except   with   redundancy. For London and , the correlation is up to 0.5.   For FAO-Trade, Q shows a higher correlation up to 1 with node-set and edge-set coverage, and a lower correlation up to 0.5 with average path length, clustering coefficient and redundancy.
 
Considering  Fig. <ref> (b), the multilayer modularity Q of the solution found by PMM shows an average correlation of 0.99 with   clustering coefficient for EU-Air, of 1 and 0.92 with   average path length, node-set and edge-set coverage for  AUCS and , respectively, and of 1 with   redundancy for VC-Graders.
 On the contrary, a correlation of -1 is obtained between Q and the clustering coefficient and redundancy for AUCS, and between Q and all measures except the redundancy for VC-Graders.  For FAO-Trade, Q shows a positive correlation up to 1 with average path length, node-set and edge-set coverage, and a negative correlation up to -1 with the redundancy. A weakly negative correlation is shown between Q and clustering coefficient. In the other cases, the correlation ranges between -0.5 and 0.5.
  
  Finally, considering Fig. <ref> (c), the multilayer modularity Q of the solution found by shows a very high   correlation   with node-set and edge-set coverage in all networks. Also, Q shows a correlation with the average path length which is up to 0.5 in all networks, with the only exception of EU-Air and GH-SO-TW. For redundancy and clustering coefficient,  Q obtains a high correlation   with clustering coefficient and redundancy in FAO-Trade and with clustering coefficient in VC-Graders. 
 
 


    Figure <ref> shows the correlation between various settings of  Q and the previously analyzed set of statistics for solutions obtained by LART. Looking at the plots,  Q obtains the highest correlation with the edge-set coverage, followed by the node-set coverage, clustering coefficient and redundancy.  
    Overall, results by LART confirm the trends observed for GL and PMM, with even higher tendency to positive correlation in general. Remarkably, this particularly  holds when Q involves the inter-layer coupling terms, with IC_a leading to higher correlation than IC_s.
     
   


 


 §.§ Evaluation with ordered layers
 
In this section we focus on   evaluation scenarios that  correspond to the specification of an ordering of the set of layers. 
 We will present results on  
the real-world networks EU-Air and Flickr. The former was chosen because of its highest dimensionality (i.e., number of layers) over all datasets, the latter is a time-evolving multilayer network and   was chosen for evaluating the time-aware asymmetric inter-layer coupling. 
 
 





 
 
  





 

 


Table <ref> summarizes results by   GL, LART and on EU-Air, corresponding to adjacent and succeeding-layer coupling.  
We observe that, regardless of the setting of the resolution factor, values of Q with  succeeding-layer coupling are higher than the corresponding ones for the adjacent layer coupling scheme. This   suggests that  the impact on the inter-layer coupling term is higher when all ordered pairs of layers are taken into account, than when only adjacent pairs are considered. In this regard,  recall that the total degree of the multilayer graph, which normalizes the inter-layer coupling term as well, is properly computed according to the actual number of inter-layer couplings considered,  depending on whether   adjacent or   succeeding-layer scheme was selected. 

 The above result is also confirmed by PMM, as shown in Fig. <ref>, where the plots for the succeeding-layer scheme superiorly bound those for the adjacent scheme, over the various k. 
 Note also that, while results on EU-Air are shown only for  the ascendent layer ordering, by inverting this order we will have    a switch between results corresponding to the inner asymmetric case with results corresponding to the outer asymmetric case.  

   Moreover, Fig. <ref>(b) compares  the effect of  asymmetric inter-layer coupling on Flickr with and without time-awareness, for PMM solutions. Here we observe that  both IC^t_ia and IC^t_oa plots are above those corresponding to IC_ia and IC_oa.   This indicates that considering a smoothing   term for the temporal distance between layers
   (Eq. (<ref>))   leads to an increase in modularity. 
   This general result is also confirmed by GL, LART and (results not shown); for instance, GL achieved on Flickr  modularity 
 0.462 for IC^t_ia,  0.468 for IC^t_oa, and 0.460 for IC^t_s, which compared to results shown in Table <ref>  represent increments in Q of    43%. Similarly, obtained on Flickr  modularity 
 0.975 for IC^t_ia, IC^t_oa, and IC^t_s, which is higher than the corresponding values reported in Table <ref> for IC_a and IC_s, respectively.  
 
 
 
  
Correlation analysis.  
Analogously to correlation analysis performed for the unordered case,  
 we compare different settings of Q with selected  statistics on topological properties.  Figures <ref> 
  show results on   EU-Air   
   obtained by GL, LART, PMM and , respectively.  Again, for PMM, k was set to the number of communities corresponding to the best modularity performance achieved by the method.  
 


As a general remark  Q is always non-negatively  correlated with all topological measures. More specifically,  the correlation is highly positive with all measures, when GL is used,  and with all measures but average path length and clustering coefficient, when LART and are used; for PMM, correlation is very high with clustering coefficient, and mid-low with the other measures. 
 

 When equipped with   succeeding-layer coupling, 
 correlation is higher than in the   adjacent-layer setting with average path length (up to +0.14), node-set coverage (up to +0.02) and redundancy (up to +0.05) for the solution found by GL and , and with average path length (up to +0.07), node-set coverage (up to +0.11) and edge-set coverage (up to +0.11) for the solution found by PMM. 
 
  We also found that  the layer ordering does not provide meaningful variations on the correlation values — plots regarding descendent layer ordering are reported in   the Appendix. 
  
  
 



 
                        
                        
                        
                        





                        
 


 §.§ Analysis of Q_ms and comparison with Q
 


We discuss here performance results obtained by the community detection algorithms with   Q_ms as assessment criterion. We will refer to the default setting of unordered set of layers as stated in <cit.>. 
  
Using GL,  Q_ms tends to decrease as γ increases (while ω decreases, as it was varied with γ as ω=1-γ). This occurs monotonically in most datasets, within positive ranges (e.g.,    from 0.636 to 0.384 on ,  from 0.525 to 0.391 on ) or including negative modularity for higher γ (e.g., 
from 0.645 to -0.05 on Flickr, from 0.854 to -4 on AUCS). 
Remarkably, the simultaneous effect of  γ and  ω=1-γ  on  Q_ms leads on some datasets (Obama, EU-Air, London) to  a drastic degradation of modularity (down to much negative values) for some     γ >1, followed by a rapid increase to modularity of 1 as γ increases closely to 2.   

 Analogous considerations hold for LART, PMM and . For the latter method, the trend of drastic degradation of modularity followed by a rapid increase is only visible for EU-Air. For PMM, the plots of Fig. <ref> show results by varying k, on the real-world datasets.
  Surprisingly, it appears that Q_ms is  relatively less sensitive to the variation in the community structure than our Q. 
 
  This is particularly visible in AUCS, London (not shown),  and Obama 
   where Q_ms shows no variation for increasing k. Also, it is worth noting that, for specific values of γ, Q_ms may have an abrupt decrease  with very low peaks, as it happens   for 
    Obama.  By contrast,  for FAO-Trade, a value around 0.8 for γ induces high modularity   which is stable for the different k values.
  
  

When varying ω within [0..2], with γ=1, Q_ms tends to monotonically increase as ω increases.   This holds  consistently on all datasets (Fig. <ref> shows plots for some of them) with the only exception of FAO-Trade, where Q_ms is stable at 1 for ω > 0.8. Variations are always on positive intervals (e.g.,   from 0.248 to 0.621 on Flickr, from 0.305 to 0.541 on , from 0.136 to 0.356 on Higgs-Twitter). 



0.5em


  §.§.§ Comparison between Q and Q_ms: qualitative evaluation on the solutions generated by the community detection methods
   

In the light of the above analysis, 

 a few interesting remarks arise by observing the different behavior of Q and Q_ms over the same  community structure solutions, in function of the resolution and inter-layer coupling factors. From a qualitative viewpoint, the effect  of ω on  Q_ms turns out to be opposite, in most cases, to the effect of our IC terms on Q: that is, accounting more for the inter-layer couplings leads to increase    Q_ms, while this does not necessarily happen in Q. 
Less straightforward is comparing the use of a constant resolution for all layers, as done in Q_ms, and the use of variable (i.e., for each pair of layer and community)  resolutions, as done in Q. In this regard, we have previously observed that the use of a varying redundancy-based resolution factor improves Q w.r.t. the setting γ=1. By coupling this general remark with the results (not shown) of an inspection of the values of γ(L,C)   in the computation of Q on the various network datasets (which confirmed   that γ(L,C) values span over its range,  in practice), we can conclude that a more appropriate   consideration of the term modeling the expected connectivity of community is realized in our Q w.r.t. keeping the resolution as constant for all layers in Q_ms.


  


0.5em


  §.§.§ Comparison between Q and Q_ms: evaluation on AUCS  ground-truth  communities
   

Let us now consider a further   stage of evaluation of Q vs. Q_ms, which is complementary to the previous comparative analysis, with the specific purpose of assessing their behavior w.r.t. a ground-truth community structure.  To this aim, we resorted again to the AUCS data: in their original work <cit.>, the authors filled a gap in the literature (actually, still largely open) corresponding to a lack of benchmarks for understanding multilayer/multiplex networks. In that work,  the authors also provided a ground-truth multiplex community structure for AUCS, which reflects  affiliation of the university employees/students to research groups. Please refer to <cit.> for a detailed description of how this ground-truth was obtained. 

 

[ph!]
Comparison between specific settings of  Q and Q_ms for each of the non-singleton  ground-truth communities  of   AUCS  

 
   #nodes/    avg.      avg.    clust.   2c||Q_ms   6c||Q   3c|Q_ms

 6-16
    #edges    degree    path     coeff.    2c||γ=1   2c|γ=1   4c||γ(L,C)   γ=1   2c|γ= γ(L,C)

 6-16
            length         ω=1   ω=2   IC_s   IC_ia   IC_s   IC_ia   IC^Adj_ia   IC^Suc_ia    ω= IC_s    ω=1    ω= IC_s

 C_1   4/8	   	3.000   	1.081   	0.523   0.076   0.083   0.037   0.052   0.038   0.054   0.073   0.069   0.064   0.078   0.067
C_2   11/29	   4.853   	1.564   	0.635    0.202   0.239   0.118   0.131   0.135   0.148   0.223   0.229   0.145   0.226   0.181
C_3   6/4	   1.210   	1.148   	0.120   0.052   0.061   0.020   0.035   0.020   0.036   0.044   0.039   0.038   0.053   0.040
C_4   5/6   2.105   	1.038   	0.413	   0.062   0.069   0.029   0.043   0.029   0.044   0.057   0.054   0.051   0.063   0.053
C_5   3/4   2.100   	1.133   	0.600	   0.054   0.066   0.018   0.034   0.018   0.035   0.041   0.034   0.036   0.056   0.038
C_6   6/8	   2.430   	1.208   	0.428   0.083   0.095   0.036   0.054   0.037   0.055   0.072   0.067   0.065   0.085   0.068
C_7   7/13   3.570   	1.551   	0.659   0.123   0.141   0.057   0.074   0.058   0.075   0.103   0.101   0.096   0.128   0.103

5|r|Global modularity        0.651   0.754   0.314   0.424   0.337   0.447   0.615   0.593   0.495   0.689   0.550


  



In this context of evaluation, we analyzed  again the behaviors of  Q and  Q_ms under particular settings of the resolution and inter-layer coupling factors, while keeping fixed the community structure to a reference one corresponding to  ground-truth  knowledge. 
 
 Before going into details of such analysis, let us first provide some remarks on the trends of Q_ms   under its two   settings previously considered in Sect. <ref>:   varying γ within [0..2] with ω=1-γ, and   varying ω within [0..2] with γ=1.  
  As we observe from the results shown in 
Fig. <ref>, Q_ms can vary significantly depending on the setting of γ, ω: when γ=1, ω∈ [0..2], 
Q_ms monotonically increases with ω, varying within a relatively small range (i.e., from 0.38 to 0.76); however,  when γ∈ [0..2],  ω=1-γ, Q_ms follows an inverse trend, with a more rapid decrease for γ > 1.5, and overall wider range (from 0.768 to -1.09).  Note that these considerations on the trends are consistent with the previous analysis for Q_ms computed on the solution found by PMM for AUCS (cf. Figs. <ref>–<ref> (e)). 
 We shall come back later on such a high parameter-sensitivity of Q_ms.    
 

Table <ref> shows the global and community-specific values of Q and Q_ms for particular combinations of their corresponding    γ and inter-layer coupling factors (IC and ω, respectively). 
Columns 2 to 5 report basic structural statistics for each of the communities, while the rest of the table is organized into three subtables: the first refers to  Q_ms results, the second to our Q, and the third again to Q_ms with Q-biased settings of γ, ω. In the latter case, we wanted to make Q_ms “closer to”  Q by setting  its parameters to the values of γ(L,C) and IC_s, respectively, averaged over the different communities and layers.  

Moreover, note that ground-truth communities in AUCS are 14 in total, however we only reported  results for those (7) that contain more than one node, in order to  avoid  cluttering the table with roughly constant, zero-close modularity values that correspond to the singleton communities. 



Looking at the table, communities C_2 and C_7 (resp. C_5, C_3, and C_4) correspond to the highest (resp. lowest) modularity values, for either modularity measure, under each of the parameter-combinations considered. 
 In general, beyond the differences in the respective values of modularity (note that  Q_ms values still differ  from the corresponding Q ones even for the Q-biased settings), the two measures appear to  behave  similarly  over the various communities.  
  To confirm this intuition, we evaluated the Pearson correlation of different pairings of Q and Q_ms community-specific values. 
  Results show indeed almost perfect correlation (always above 0.98).  


One aspect of evaluation that we also investigated is whether disrupting the multiplex network by layer may have effect on the comparison between 
 Q and  Q_ms, which resembles a sort of layer-oriented  resiliency analysis. 
  More specifically, based on the structural impact due to the various layers in AUCS <cit.>, we considered  the following alternative configurations:  
  (i) we removed layer co-authorship   (i.e., the smallest and less connected of all layers), 
  (ii) we removed layers co-authorship and leisure    (i.e., the ones  having the lowest number of edges), 
    (iii) we removed layer work   (i.e., the one with the most edges), 
  (iv) we retained layers work and lunch   only.   
 



 For each of such multiplex-disruption configurations, we replicated the above analysis corresponding to the full multiplex. 
  Results  (not shown)  indicated no particular differences in terms of rankings of community modularity obtained   by  Q_ms and Q, respectively; however, we also observed a general decrease in Pearson correlation of the pairings of Q_ms and Q community-specific values, although the correlation remained still high, in particular always above 0.96 (e.g., when removing  layer work, correlation was 0.965 between Q_ms with γ=1, ω= IC_s and Q with γ=1, IC_s, and 0.98 between Q_ms with γ= γ(L,C), ω= IC_s and Q with varying γ(L,C), IC_s). 
 

 
To sum up, in this ground-truth evaluation,  Q_ms and Q exhibited consistently  similar behaviors at community level for specific settings of their respective parameters of resolution and inter-layer coupling. However, it should be emphasized that such a similarity between the two modularity values was actually achieved for either canonical settings of γ, ω in Q_ms (i.e., γ=1 and ω∈{1, 2}) or Q-biased settings of γ, ω (i.e., γ= γ(L,C), ω= IC_s). In general, Q_ms has shown to be highly sensitive to the settings of its parameters, whereas by contrast, our Q has the  key advantage of  automatically determining the resolution and inter-layer coupling factors   based on the structural information of the communities in the multilayer network.    

 



 
 



 
 






 

 §.§ Efficiency evaluation

We analyzed the computation time of Q for the different combinations of redun­dancy-based resolution factor γ and inter-layer coupling factors IC. 
The 3-D plots in Figs. <ref>–<ref> display the   time vs. the number of layers and the number of nodes per layer.   For this analysis, we referred to  the solutions found by GL for the  multilayer networks  generated through the LFR benchmark  (cf. Section <ref>).  



As expected, the computation time increases with both the number of layers and nodes per layer, with the latter being less evident when setting γ fixed to 1. Also, while it is  obvious that the computation time is higher when using variable (i.e., redundancy-based) resolution factor than in the case γ=1 (with percentage increase of 50%, for the maximum number of layers and nodes per layer), we observe   much less fluctuations in the plot surfaces  than in the case of resolution factor fixed to 1, regardless of the setting of inter-layer coupling factors. 
 The reader is also referred to the Appendix for further results by varying the inter-layer coupling settings. 

 
 
 







 




§ CONCLUSION


We proposed a new definition of modularity for multilayer networks. Motivated by the  opportunity  of  revising the multislice modularity proposed in <cit.>,  
 we conceived  alternative  notions of layer resolution and inter-layer coupling, which are key-enabling for  generalizations of modularity for multilayer networks.  

Using four state-of-the-art methods for multilayer community detection, synthetic multilayer networks and ten real-world multilayer networks, we provided empirical evidence of the significance of our proposed modularity. 
 
Our work paves the way for the development of new optimization methods of community detection in multilayer networks  which, by embedding our multilayer modularity, can discover     community structures having the interesting   properties relating to the proposed per-layer/community redundancy-based resolution factors and projection-based inter-layer coupling schemes.   
 In this respect, we point out that our multilayer modularity is able to cope with communities that are overlapping at entity level,  which eventually  reflect   the different roles that the same entity can play when occurring in two or more layers of the network. 
  Within this view, one benefit of adopting the multilayer network model is that the problem of computing soft community-memberships of entities can be   translated into a simpler problem of identification of crisp community-memberships of nodes within each layer.  
 Nonetheless,   a further  interesting direction would be to evaluate our multilayer modularity in contexts of node-overlapping communities. In this case, however,  one challenge to face is whether and to  what extent an   overlapping-aware multilayer modularity should be able to measure the community  overlaps    within each layer and/or across the layers. 
 Along this direction, it would be interesting to study an integration   of our multilayer modularity into recently developed works that propose probabilistic representations or stochastic generative models for overlapping community detection  in   multilayer networks <cit.>. 

 
 
 





§ ACKNOWLEDGEMENTS

 This research work has been partly supported by the Start-(H)Open POR Grant No. J28C17000380006 funded by the Calabria Region Administration, and by the NextShop PON Grant No. F/050374/01-03/X32 funded by the Italian Ministry for Economic Development. 
   
 


10

url@samestyle








AmelioT17
A. Amelio and A. Tagarelli, “Revisiting Resolution and Inter-Layer Coupling
  Factors in Modularity for Multilayer Networks,” in Proc. IEEE/ACM
  ASONAM, 2017, pp. 266–273.

Newman04
M. E. J. Newman and M. Girvan, “Finding and evaluating community structure in
  networks,” Phys. Rev. E, vol. 69, no. 2, p. 026113, 2004.

Newman04b
M. E. J. Newman, “Fast algorithm for detecting community structure in
  networks,” Phys. Rev. E, vol. 69, 2004.

BrandesDGGHNW08
U. Brandes, D. Delling, M. Gaertler, R. Görke, M. Hoefer, Z. Nikoloski,
  and D. Wagner, “On modularity clustering,” IEEE Trans. Knowl. Data
  Eng., vol. 20, no. 2, pp. 172–188, 2008.

ChenKS14
M. Chen, K. Kuzmin, and B. K. Szymanski, “Community detection via maximization
  of modularity and its variants,” IEEE Trans. Comput. Social Syst.,
  vol. 1, no. 1, pp. 46–65, 2014.

Clauset04
A. Clauset, M. E. J. Newman, and C. Moore, “Finding community structure in
  very large networks,” Phys. Rev. E, vol. 70, 2004.

Newman05
M. E. J. Newman, “Modularity and community structure in networks,”
  Proc. Natl. Acad. Sci., pp. 8577––8582,
  2005.

WhiteS05
S. White and P. Smyth, “A spectral clustering approach to finding communities
  in graph,” in Proc. SDM, 2005.

Reichardt06
J. Reichardt and S. Bornholdt, “Statistical mechanics of community
  detection,” Phys. Rev. E, vol. 74, 2006.

Duch05
J. Duch and A. Arenas, “Community detection in complex networks using extremal
  optimization,” Phys. Rev. E, vol. 72, 2005.

BOCCALETTI20141
S. Boccaletti, G. Bianconi, R. Criado, C. I. del Genio, J. G.-G. nes,
  M. Romance, I. S. na Nadal, Z. Wang, and M. Zanin, “The structure and
  dynamics of multilayer networks,” Phys. Rep., vol. 544, no. 1,
  pp. 1–122, 2014.

Magnanibook
M. E. Dickison, M. Magnani, and L. Rossi, Multilayer social
  networks.1em plus 0.5em minus 0.4emCambridge University Press,
  2016.

Kivela+14
M. Kivela, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A.
  Porter, “Multilayer networks,” J. Complex Netw., vol. 2,
  no. 3, pp. 203–271, 2014.

e20120909
L. G. A. Alves, G. Mangioni, F. A. Rodrigues, P. Panzarasa, and Y. Moreno,
  “Unfolding the complexity of the global value chain: Strength and entropy in
  the single-layer, multiplex, and multi-layer international trade networks,”
  Entropy, vol. 20, no. 12, 2018.

Mucha10
P. J. Mucha, T. Richardson, K. Macon, M. A. Porter, and J.-P. Onnela,
  “Community structure in time-dependent, multiscale, and multiplex
  networks,” Science, vol. 328, no. 5980, pp. 876–878, 2010.

coscia2011
M. Coscia, F. Giannotti, and D. Pedreschi, “A classification for community
  discovery methods in complex networks,” Stat. Anal. Data
  Min., vol. 4, no. 5, pp. 512–546, 2011.

KunchevaM15
Z. Kuncheva and G. Montana, “Community detection in multiplex networks using
  locally adaptive random walks,” in Proc. IEEE/ACM ASONAM, 2015, pp.
  1308–1315.

ZhangWLY16
H. Zhang, C. Wang, J. Lai, and P. S. Yu, “Modularity in Complex Multilayer
  Networks with Multiple Aspects: A Static Perspective,” CoRR, vol.
  abs/1605.06190, 2016.

Tagar17
A. Tagarelli, A. Amelio, and F. Gullo, “Ensemble-based community detection in
  multilayer networks,” Data Min. Knowl. Discov., vol. 31, no. 5, pp.
  1506–1543, 2017.

Tagar18
D. Mandaglio, A. Amelio, and A. Tagarelli, “Consensus Community Detection in
  Multilayer Networks Using Parameter-Free Graph Pruning,” in Proc. PAKDD, 2018, pp. 193–205.

torreggiani2018identifying
S. Torreggiani, G. Mangioni, M. J. Puma, and G. Fagiolo, “Identifying the
  community structure of the food-trade international multi-network,”
  Environ. Res. Lett., vol. 13, no. 5, p. 054026, 2018.

mangioni2018multilayer
G. Mangioni, G. Jurman, and M. De Domenico, “Multilayer flows in molecular
  networks identify biological modules in the human proteome,” IEEE
  Trans. Netw. Sci. Eng., 2018.

alves2019nested
L. G. Alves, G. Mangioni, I. Cingolani, F. A. Rodrigues, P. Panzarasa, and
  Y. Moreno, “The nested structural organization of the worldwide trade
  multi-layer network,” Sci. Rep., vol. 9, no. 1, p. 2866, 2019.

FortunatoB07
S. Fortunato and M. Barthelemy, “Resolution limit in community detection,”
  Proc. Natl. Acad. Sci., vol. 104, no. 36, 2007.

FORTUNATO201075
S. Fortunato, “Community detection in graphs,” Phys. Rep., vol.
  486, no. 3, pp. 75–174, 2010.

2010Tang
L. Tang and H. Liu, Community Detection and Mining in Social Media, ser.
  Synthesis Lectures on Data Mining and Knowledge Discovery.1em plus
  0.5em minus 0.4emMorgan & Claypool Publishers, 2010.

FORTUNATO20161
S. Fortunato and D. Hric, “Community detection in networks: A user guide,”
  Phys. Rep., vol. 659, pp. 1–44, 2016.

leicht2008
E. A. Leicht and M. E. J. Newman, “Community structure in directed networks,”
  Phys. Rev. Lett., vol. 100, p. 118703, Mar 2008.

arenas2007
A. Arenas, J. Duch, A. Fernández, and S. Gómez, “Size reduction of
  complex networks preserving modularity,” New J. Phys.,
  vol. 9, no. 6, p. 176, 2007.

newman2004
M. E. J. Newman, “Analysis of weighted networks,” Phys. Rev. E,
  vol. 70, p. 056131, Nov 2004.

XIANG20124995
J. Xiang and K. Hu, “Limitation of multi-resolution methods in community
  detection,” Physica A Stat. Mech. Appl.,
  vol. 391, no. 20, pp. 4995–5003, 2012.

Zhang2009
J. Zhang, K. Zhang, X. ke Xu, C. K. Tse, and M. Small, “Seeding the kernels in
  graphs: toward multi-resolution community analysis,” New J. Phys., vol. 11, no. 11, p. 113003, 2009.

arenas2007-2
A. Arenas, A. Fernández, and S. Gómez, “Analysis of the structure of
  complex networks at different resolution levels,” New J. Phys., vol. 10, no. 5, p. 053039, 2008.

XIANG2015127
J. Xiang, Y.-N. Tang, Y.-Y. Gao, Y. Zhang, K. Deng, X.-K. Xu, and K. Hu,
  “Multi-resolution community detection based on generalized self-loop
  rescaling strategy,” Physica A Stat. Mech. Appl., vol. 432, pp. 127–139, 2015.

gomez2009
S. Gómez, P. Jensen, and A. Arenas, “Analysis of community structure in
  networks of correlated data,” Phys. Rev. E, vol. 80, p. 016114, 2009.

traag2009
V. A. Traag and J. Bruggeman, “Community detection in networks with positive
  and negative links,” Phys. Rev. E, vol. 80, p. 036115, Sep 2009.

barber2007
M. J. Barber, “Modularity and community detection in bipartite networks,”
  Phys. Rev. E, vol. 76, p. 066102, 2007.

barber2008
Y. Xu, L. Chen, and S. Zou, “Community detection from bipartite networks,” in
  Proc. IEEE WISA, 2013, pp. 249–254.

guimera2007
R. Guimerà, M. Sales-Pardo, and L. A. N. Amaral, “Module identification in
  bipartite and directed networks,” Phys. Rev. E, vol. 76, p. 036102,
  2007.

nicosia2009
V. Nicosia, G. Mangioni, V. Carchiolo, and M. Malgeri, “Extending the
  definition of modularity to directed graphs with overlapping communities,”
  J. Stat. Mech., vol. 2009, no. 03, p. P03024, 2009.

Berlingerio2011
M. Berlingerio, M. Coscia, and F. Giannotti, “Finding and characterizing
  communities in multidimensional networks,” in Proc. IEEE/ACM ASONAM,
  2011, pp. 490–494.

KimL15
J. Kim and J. Lee, “Community detection in multi-layer graphs: A survey,”
  SIGMOD Record, vol. 44, no. 3, pp. 37–48, 2015.

Rossi2015
L. Rossi and M. Magnani, “Towards effective visual analytics on multiplex and
  multilayer networks,” Chaos, Solitons & Fractals, vol. 72, pp.
  68–76, 2015.

Domenico2015
M. De Domenico, V. Nicosia, A. Arenas, and V. Latora, “Structural
  reducibility of multilayer networks,” Nat. Commun., vol. 6,
  p. 6864, 2015.

cha2009www
M. Cha, A. Mislove, and K. P. Gummadi, “A Measurement-driven Analysis of
  Information Propagation in the Flickr Social Network,” in Proc. WWW,
  2009, pp. 721–730.

kdweb2015
G. Silvestri, J. Yang, A. Bozzon, and A. Tagarelli, “Linking accounts across
  social networks: the case of stackoverflow, github and twitter,” in
  Proc. Work. on Knowl. Discov. on the Web, 2015, pp. 41–52.

Omodei2015
E. Omodei, M. De Domenico, and A. Arenas, “Characterizing interactions in
  online social networks during exceptional events,” Frontiers in
  Physics, vol. 3, p. 59, 2015.

TangWL09
L. Tang, X. Wang, and H. Liu, “Uncovering groups via heterogeneous interaction
  analysis,” in Proc. IEEE ICDM, 2009, pp. 503–512.

Stanley2016
N. Stanley, S. Shai, D. Taylor, and P. J. Mucha, “Clustering network
  layers with the strata multilayer stochastic block model,” IEEE
  Trans. Netw. Sci. Eng., vol. 3, no. 2, pp. 95–105, April 2016.

Ali2019
H. T. Ali, S. Liu, Y. Yilmaz, R. Couillet, I. Rajapakse, and A. O.
  Hero, “Latent heterogeneous multilayer community detection,” in
  Proc. ICASSP, 2019.


 










§ APPENDIX







§ ANALYTICAL DERIVATION OF LOWER AND UPPER BOUNDS





 Proof of Proposition 1 (Lower bound of Q).    
Let us assume that each of the ℓ layer graphs in G_ℒ has the form of a bipartite graph  K_i(a, b), with i=1, …, ℓ, and sets a and b induce a partitioning of the set of nodes in two communities denoted as C_1 and C_2, respectively, so that 𝒞={C_1, C_2} with |C_1|=|C_2|=n/2, and no internal links are drawn between nodes of the same community (because of the bipartite assumption).

To begin with, consider the reduction of Q to its simplest form, i.e., 
  γ(L,C) fixed to 1 for any L,C and β=0. Therefore, the contribution of community C_1 to Q is: 









    Q(C_1) =  -∑_L ∈ℒ(d_L(C_1)/d(V_ℒ))^2.


Since d_L(C_1)=n^2/4 and d(V_ℒ)=n^2/4(2ℓ)=n^2 ℓ/2, Q(C_1) is calculated as:


    Q(C_1) = -∑_L ∈ℒ(n^2/4/n^2ℓ/2)^2=-1/4ℓ.


The same above holds for C_2. Therefore, the minimum value of Q when γ(L,C)=1 and β=0 is as follows:  

    Q(𝒞) = -1/4ℓ2=-1/2ℓ.


Let us     now consider Q with its  redundancy-based resolution factor while keeping β=0. 






 Since the internal degree of community C_1 is 0,  there are no redundant pairs for any community and layer, and hence  γ(L,C)=2. Consequently,  the contribution of  C_1 to  Q is: 


    Q(C_1) =  -2 ∑_L ∈ℒ(d_L(C_1)/d(V_ℒ))^2 = 
    -2 ∑_L ∈ℒ(1/2ℓ)^2 = -1/2ℓ.


The same above holds for C_2. Therefore, the lower bound of Q when γ(L,C) is variable and β=0 is as follows:  


    Q(𝒞) =-1/2ℓ2=-1/ℓ.



In the general form of Q with both resolution and inter-layer coupling factors (i.e., varying γ and β=1), the contribution of C_1 to Q is:


    Q(C_1) =  ∑_L ∈ℒ[-2(d_L(C_1)/d(V_ℒ))^2 + ∑_L' ∈𝒫(L)IC(C_1, L, L')/d(V_ℒ)].


 In the above formula, let us indicate the terms  X(C_1)=∑_L ∈ℒ -2(d_L(C_1)/d(V_ℒ))^2, and Y(C_1)=∑_L ∈ℒ, L' ∈𝒫(L)IC(C_1, L, L')/d(V_ℒ). 
If we denote with p= ∑_L ∈ℒ|𝒫(L)|  the total number of valid layer-pairings, then   d(V_ℒ)=n^2ℓ/2+np, and X(C_1) can be rewritten  as follows:


    X(C_1)=-2∑_L ∈ℒ(n^2/4/n^2ℓ/2+ np)^2=-1/2n^2ℓ/(nℓ + 2p)^2.


For the reduction of the term Y(C_1), we consider the two cases of symmetric inter-layer coupling and asymmetric inter-layer coupling. 
In the first case, the minimum value for IC_s(C_1,L_i, L_j) is equal to |C_1^(i)∩ C_1^(j)|/|V_i ∩ V_j|=1/n.  

Accordingly, Y(C_1) is reduced as follows:








    Y(C_1)  =  ∑_L ∈ℒ∑_L' ∈𝒫(L)1/n/n^2ℓ/2+ np = 2p/n^2(nℓ+2p)


In the second case, IC_a(C_1,L_i, L_j) = |C_1^(i)∩ C_1^(j)|/|V_i ∩ V_j|×|V_i|/|C_1^(i)|=1/n×n/n/2=2/n. 

Accordingly, Y(C_1) is reduced as follows:


    Y(C_1)=∑_L ∈ℒ∑_L' ∈𝒫(L)2/n/n^2ℓ/2+ np=4p/n^2 (nℓ   +  2p).


The  above expressions for X and Y hold for C_2. Therefore, the lower bound of Q in its general form is as follows:   


    Q(𝒞) =     2(X(C_1)+Y(C_1)) = 
    
     =     -n^2ℓ/(nℓ +2p)^2 + 4(1+η)p/n^2(nℓ + 2p),


with η=0 for IC_s and η=1 for IC_a.

It should be noted that Eq. <ref> is a special case of  Eq. <ref> with β=0 and d(V_ℒ) discarding the contribution given by the inter-layer edges (i.e., d(V_ℒ) = n^2ℓ/2).































































 
 Proof of Proposition 2 (Upper bound of Q). 
Let us assume that each of the ℓ layer graphs in G_ℒ has community structure 𝒞={C_1, C_2} such that |C_1|=|C_2|=n/2 and each community is a clique with n/2(n/2-1)/2 edges. Moreover, there are no external edges connecting the communities, therefore  d_L(C)=d_L^int(C). 
 Note that by uniformly distributing the n nodes into the two communities, it can easily be shown that the  maximum of Q is higher.

Analogously to the  analysis of the minimum value of Q, let us first consider the setting γ(L,C)=1, for any L,C, and β=0. 





 
  Therefore, the contribution of community C_1 to Q is: 
   



    Q(C_1) =   ∑_L ∈ℒd_L^int(C_1)/d(V_ℒ)  -  (d_L^int(C_1)/d(V_ℒ))^2.

 


Because d_L^int(C_1)=n/2(n/2-1) and d_L()=n/2(n/2-1)2ℓ, the above expression is rewritten as:
 
 
 

    Q(C_1) =   ∑_L ∈ℒn/2(n/2-1)/n/2(n/2-1)2ℓ  -  (n/2(n/2-1)/n/2(n/2-1)2ℓ)^2 = 1/2-1/4ℓ.


The same above holds for C_2. Therefore, the maximum of Q when γ(L,C)=1 and β=0 is as follows:[If C_1 and C_2 would have  n-1 and 1 nodes, respectively, the resulting maximum value of Q will be ℓ-1/ℓ, hence lower than what we obtain in Eq. <ref>.]  


    Q(𝒞) =   2(1/2-1/4ℓ)=1-1/2ℓ=2ℓ-1/2ℓ.


Consider now the setting with redundancy-based γ(L, C) and β=0. The contribution of C_1 to Q is calculated as: 







    Q(C_1) =   ∑_L ∈ℒn/2(n/2-1)/n/2(n/2-1)2ℓ  -  γ(L, C_1)(n/2(n/2-1)/n/2(n/2-1)2ℓ)^2.


Since nrp(L, C_1)=n/2(n/2-1)/2, it follows that  γ(L, C_1) is   equal to 2(1+log_2(1+n/2(n/2-1)/2))^-1 for each layer and community. 

 Note that the above constant quantity, hereinafter denoted as   γ, tends to be ≪ 1, and  it is  smaller for higher number of nodes n. 

The contribution of C_1 to Q is rewritten as:
 

    Q(C_1)=1/2-γ1/4ℓ.


 
The same above holds for C_2. Therefore, the maximum  of Q with redundancy-based γ(L,C) and β=0 is as follows:  


    Q(𝒞) =   2[1/2-γ1/4ℓ]=1-γ1/2ℓ=2ℓ-γ/2ℓ.


 Note that the above quantity is as much closer to 1 as  n and ℓ are higher. 



In the general setting of Q with redundancy-based   γ(L, C) and β=1, the contribution of C_1 to Q is: 


    Q(C_1) =    ∑_L ∈ℒ[  d_L^int(C_1)/d(V_ℒ)  - γ(L, C) (d_L^int(C_1)/d(V_ℒ))^2  + .
    
     +    . ∑_L' ∈𝒫(L)IC(C_1, L, L')/d(V_ℒ)].


 In the above formula, let us indicate the terms  X(C_1)=∑_L ∈ℒ[  d_L^int(C_1)/d(V_ℒ)  - γ(L, C) (d_L^int(C_1)/d(V_ℒ))^2], 
 and Y(C_1)=∑_L ∈ℒ ∑_L' ∈𝒫(L)IC(C_1, L, L')/d(V_ℒ).  

 Because d_L^int(C_1)=n/2(n/2-1), d_L()=n/2(n/2-1)2ℓ+np, and γ(L, C_1)=2(1+log_2(1+ n/2(n/2-1)/2))^-1 for each layer and community, we obtain:


    X(C_1)=   ∑_L ∈ℒn/2(n/2-1)/n/2(n/2-1)2ℓ+np  -  γ(n/2(n/2-1)/n/2(n/2-1)2ℓ+np)^2= 
    
     =    1/2(n/2-1)ℓ/(n/2-1)ℓ+p-γℓ(1/2(n/2-1)/(n/2-1)ℓ+p)^2,



    Y(C_1)=∑_L ∈ℒ∑_L' ∈𝒫(L)IC(C_1, L, L')/n/2(n/2-1)2ℓ+np.


For the reduction of the term Y(C_1), we again consider the two cases of symmetric inter-layer coupling and asymmetric inter-layer coupling. 
In the first case, IC_s(C_1,L_i, L_j) = |C_1^(i)∩ C_1^(j)|/|V_i ∩ V_j|=1/n. Therefore:
 

    Y(C_1)=∑_L ∈ℒ∑_L' ∈𝒫(L)1/n/n(n/2-1)ℓ+np=p/n^2(n/2-1)ℓ+n^2p.


In the second case, IC_a(C_1,L_i, L_j) = |C_1^(i)∩ C_1^(j)|/|V_i ∩ V_j|×|V_i|/|C_1^(i)|=1/n×n/n/2=2/n. Therefore:


    Y(C_1)=∑_L ∈ℒ∑_L' ∈𝒫(L)2/n/n/2(n/2-1)2ℓ+np=2p/n^2(n/2-1)ℓ+n^2p.


The same above holds for C_2. Therefore,  the maximum value of Q  is as follows:


    Q(𝒞) =     2(X(C_1)+Y(C_1)) = 
    
     =     2 [ 1/2(n/2-1)ℓ/(n/2-1)ℓ+p-γℓ(1/2(n/2-1)/(n/2-1)ℓ+p)^2 + . 
    
     +    .  (1+η)p/n^2(n/2-1)ℓ+n^2p]


with η=0 for IC_s and η=1 for IC_a.

  It is worth noting that Eq. <ref> is a special case of  Eq. <ref>  with  β=0 and d(V_ℒ) discarding the contribution given by the inter-layer edges.





















































 
 

 

§ EVALUATION WITH ORDERED LAYERS

 
 Figure <ref> provides further details on  correlation analysis using descendent layer ordering.
 
 
 
 




  




 



§ EFFICIENCY RESULTS
 
Figures <ref>–<ref> show the computation time of Q for the different combinations of  γ and  IC factors.


