Improving temporal smoothness and snapshot quality in dynamic network community discovery using NOME algorithm

Dynamic network problem description

A dynamic network can be defined as a set of multiple static network snapshots G = {G¹, G², …, G^T}.G^t = {V^t, E^t}, where V^t is the set of nodes of G^t.E^t is the set of edges of G^t, denoting the set of edges at moment t (t ∈ [1, T]) consisting of two different nodes at that moment, and all the nodes and their connected edges form this network structure.

Dynamic network community discovery problem description

Dynamic network community detection is to find the community structure division of the network structure at each moment, it involves two main objectives in the process of finding, one is to ensure the quality of the community division in the current time step, and the other is to ensure that the community structure changes slowly in the continuous time, i.e., there is no huge change. Suppose that after dynamic network community detection, a network snapshot G^t of the community structure division is obtained at moment t, denoted as C = {C₁, C₂, …, C_k}, where $C_i^t$ denotes the community structure of the i- th community division, i ∈ [1, k]. In general, to better accomplish these two objectives, we use the Normalized Mutual Information (NMI) function to detect the similarity of two neighboring temporal community divisions and the modularity function Q to assess the quality of the current community division.

Normalized mutual information function (NMI)

NMI (Danon et al., 2005) is used to detect the similarity of community division between the current time sequence and the previous time sequence, and its formula is shown in Eq. (1). (1)${\displaystyle NMI\left(C^t,C^{t-1}\right)=\frac{-2\sum _{i=1}^{m^t}\sum _{j=1}^{m^{t-1}}{L}_{\mathit{ij}}\mathit{log}\left(\frac{L_{ij}N}{L_{i.}{L}_{.j}}\right)}{\sum _{i=1}^{m^t}{L}_{i.}\mathit{log}\left(\frac{L_{i.}}{N^t}\right)+\sum _{j=1}^{m^{t-1}}{L}_{.j}\mathit{log}\left(\frac{L_{.j}}{N^{t-1}}\right)}}$

where $C^t=\left\{C_1^t,C_2^t,\dots ,C_k^t\right\}$$C^{t-1}=\left\{C_1^{t-1},C_2^{t-1},\dots ,C_k^{t-1}\right\}$ denotes the community structure divided at the moment t (t-1) and L_ij denotes the number of nodes that are both in community $C_i^t$ and in community $C_j^{t-1}$. N is the number of nodes corresponding to C^t, m^t and m^t−1 denote the number of community divisions in $C_i^t$ and $C_j^{t-1}.L_{i.}$ and L_.j are the sum of elements in row i and column j of the confusion matrix L, respectively, NMI(C^t, C^t−1) = 1. If NMI(C^t, C^t−1) = 1, it means that C^t and C^t−1 have exactly the same community division result; conversely, 0 means that these two communities have completely different division results, so the closer the value of NMI is to 1 the more similar C^t and C^t−1 are on the surface.

Modularity function (Q)

The modularity function Q (Newman & Girvan, 2004) was proposed by Newman et al. to measure the quality of the division of the community structure, and a larger value of modularity indicates that the result is closer to the real community structure. It is defined as the expected value of the difference between the ratio of the number of edges inside the network community to the number of edges in the whole network and the ratio of the number of arbitrarily connected edges between nodes to the number of edges in the whole network under the same community structure, and its formula is shown in Eq. (2). (2)${\displaystyle Q=\frac{1}{2m}{\sum }_{i=1}^n{\sum }_{j=1}^n\left[A_{ij}-\frac{k_i{k}_j}{2m}\right]\delta \left(c_i,c_j\right)}$

where m denotes the number of edges, n is the number of network vertices, A_ij is the adjacency matrix that represents the relationship between edges before node i and node j, k_i and k_j are the degrees of nodes, and c_i is the community label of node i. If δ(c_i, c_j) = 1, node i and node j are in the same community, otherwise δ(c_i, c_j) = 0. The value range of module degree Q is [−0.5, 1].

Multi-objective optimization

The multi-objective optimization algorithm can optimize two or more conflicting objectives at the same time. We define the multi-objective problem of community partitioning at time t as (3)${\displaystyle \left\{\begin{array}{c}minF\left(C^t\right)=\left[f_1\left(C^t\right),f_2\left(C^t\right),\dots ,f_m\left(C^t\right)\right]\hfill \\ s.t.C^t\in \Omega \hfill \end{array}\right.}$

where f_i(C^t) is the i-th (i =1,2, …,m) objective function, m is the number of objective functions, C^t is the community structure at time step t, and Ω is the feasible domain of the set of all community structures at time step t.

Suppose $C_1^{\mathrm{t}}$ and $C_2^{\mathrm{t}}$ are two feasible solutions to this dynamic network optimization problem. We define the dominance relation of the solution of community partitioning at time t as Eq. (4), and $C_1^{\mathrm{t}}$ is said to dominate $C_2^{\mathrm{t}}$ when and only if all the conditions of formula (4) are satisfied. (4)${\displaystyle \left\{\begin{array}{c}\forall i=1,2,\dots ,m;\hfill \\ f_i\left(C_1^t\right)\le f_i\left(C_2^t\right)\wedge \exists j\in \left[1,m\right];\hfill \\ f_j\left(C_1^t\right)<f_j\left(C_2^t\right)\hfill \end{array}\right.}$

If a solution is not dominated by any other solution, the solution is called an non-dominated solution. All non-dominated solutions form a Pareto optimal set (PS), and the set of PS mapped to the target space is called the Pareto front (PF).

Coding method

In performing the community discovery process, the community of each individual is divided using coding. The main current encoding methods are the string encoding method and the bit-adjacent encoding method (Gong et al., 2014). In this article, a combination of string encoding and genomic matrix encoding is used for each individual.

Genome matrix coding refers to the representation of information about the community and node-to-node information in the form of a genome. The individuals for which the community finds a solution to the problem are encoded as n × n matrix genomes of order n, where n equals the number of nodes —v—. e_ij is used in this article to represent the relationship between node v_i and node v_j, while the value of e_ij proves whether v_i and v_j belong to an edge within a community or between communities. where e_ij =  − 1, indicating that v_i and v_j are connected to each other; e_ij = 0, indicating that there is no connection between the two nodes; and e_ij = 1, indicating that v_i and v_j are nodes belonging to the same community.

Figure 1A shows an original network and the original community, which can be encoded into the form of a genome matrix using the genome matrix according to the connection relationship between the network nodes, and then the genes in the matrix were updated after the reorganization (crossover, mutation, and BNO assignment) operation, and finally, the structure of the divided community was obtained by decoding. Figure 1B shows the divided network, and nodes of the same color indicate that they are in the same community. In the reconstructed matrix how to decode as string encoding method, we divide the sub-matrix genes into the same community based on the relationship of inner connection in the matrix until all the inner connection genes that are divided. As the yellow and green areas in the matrix before decoding are the two largest inner-connected sub-matrices we find, and also satisfy all inner-connected genes are divided, so the decoding operation is completed by converting all the divided sub-matrices to string encoding.

Proposed algorithm

This section proposes an improved algorithm NOME to discover dynamic community structures based on the work on evolutionary community detection in dynamic social networks proposed by Liu et al. (2019). A decomposition-based MOEA/D algorithm framework is used to improve the quality of community detection by migrating nodes with respect to their occupancy in neighboring communities. In addition, in order to obtain a better time step for the division of the network structure, NOME uses a rationalized selection scheme to obtain the final solution based on the final Pareto optimal solution set at the current moment. Algorithm 1 shows the overall flow of NOME.

Boundary node occupancy assignment

To improve the quality of community checking solutions in dynamic networks, we propose the community boundary node occupancy (BNO) assignment method, which can better classify nodes into the correct communities. This operation emphasizes the interconnection of nodes within and between communities, and the BNO is specifically calculated as (5)${\displaystyle O_c^i=\left\{\begin{array}{c}\frac{B_c^i}{N_c-1},i\in c;\hfill \\ 0,N_c=1;\hfill \\ \frac{B_c^i}{N_c},others\hfill \end{array}\right.}$

where ${\mathrm{O}}_{\mathrm{c}}^i$ means to calculate the occupancy rate of the current i-th node in the neighbor community c. c represents the c th community in the node i neighbor community. $B_c^i$ represents the total number of neighbor nodes of node i in the c th neighbor community, and N_c represents the total number of nodes in the c th community. If node i ∈ c, then the node itself needs to be subtracted when calculating the total number of nodes in the c-th community. If N_c = 1, then the occupancy of node i in the c-th neighboring community is 0. The value range of ${\mathrm{O}}_{\mathrm{c}}^i$ is between [0,1].

Figure 2 illustrates the process of BNO through an example. Where the intra-community and inter-community connections are represented by solid and dashed lines, respectively, and nodes of the same color indicate being in the same community. Figure 2A shows the structure of the divided network community. Node v3 is a boundary node of community C1. v3′s neighboring nodes are {v1, v2, v4, v5, v6, v7}, where nodes v1 and v2 belong to community C1, v4 and v5 belong to community C2, and v6 and v7 belong to community C3, so its neighboring communities are C2 and C3. For a boundary node, we need to calculate its occupancy in each neighboring community. The calculation by Eq. (5) shows that since node v3 belongs to community C1, it needs to subtract the node itself when calculating the total number of community nodes, so its total number of nodes is 3. And there are two neighbor nodes of v3 in community C1, which are v1 and v2, so the final calculation gets the occupancy of v3 in community C1 as 2/3. Similarly, the occupancy of v3 in community C2 is 1 and in community C3 is 1/2. Finally, through the node occupancy allocation method, node v3 will eventually be assigned to community C2, as in Fig. 2B.

There are three possible cases in the BNO allocation process, and we use different allocation schemes for different cases. The first one is that the occupancy rate of the community to which the boundary node belongs is higher than the occupancy rate of its neighbor community, at which point the node will not be assigned out. The second case is when there are multiple communities with the same maximum occupancy among all neighboring communities of the boundary node, we randomly assign the node to one of them. If one of the communities belongs to the node community, then we will determine whether the node is assigned to the other maximum occupancy community by the probability Pm. The third case is when the occupancy rate of the community to which the boundary node belongs is lower than the occupancy rate of its neighboring communities, in which case the node will be assigned to the neighboring community with the highest occupancy rate.

Initialization rules

Dynamic networks require us to perform initialization operations on the algorithm when dividing the network at each time step. The specific initialization operations are as follows.

Initialize the weight vector

We use a uniformly distributed weight vector generation method (Ji, Zhang & Zhou, 2020) where the population size is denoted by N and the N weight vectors can be expressed as λ¹, …, λ^N. Therefore, the population size and weight vectors can be set by expressing them as (6)${\displaystyle N=C_{H+m-1}^{m-1}}$

where N represents the population size and the number of weight vectors, m represents the number of optimization objectives, and H is the parameter to be set to solve for the desired number of populations by setting different parameters.

Initialize neighborhood

The T weight vectors closest to each weight vector are found by computing the Euclidean distance between any two weight vectors. The neighbors of each weight vector obtained are stored in B (i) for all i =1 to N, B(i) = i₁, i₂, …, i_T.

Initialize population

The Physarum-based network model PNM (Gao et al., 2018b) is used to initialize the population, PNM is an improvement on the mathematical model named PM (Tero, Kobayashi & Nakagaki, 2007), whose main improvement is the modification of single-entry single-exit to single-entry multiple-exit in the selection scheme for choosing whether the vertices are entrances or exits. The core mechanism of PNM is the feedback relationship between the cytoplasmic flux and the conductivity of the tube based on Poiseuille’s law. The method can improve the efficiency of the community detection algorithm by identifying the inter-community edges from the intra-community edges in the network and further enhancing the positive feedback of the solving process in the algorithm. Therefore using this method for initializing populations can give a better structure of the network community at first.

Initialize reference points

By calculating the fitness value of the initialized population, the maximum value of fitness of each target is used as the reference point Z∗. In this article, we mainly use two target module degrees Q and NMI, so Z∗ can be expressed as $Z\ast =\left\{NM{I}_{\max },Q_{\max }\right\}$.

Reorganization rules

The partitioning of dynamic network communities is solved in the decomposition-based MOEA/D framework, whose reorganization part uses a genomic matrix-based crossover and variation strategy combined with our proposed method of BNO assignment to find community solutions with high quality and high temporal smoothness.

Update rules

In the dynamic network community partitioning process, we use the relevant update strategy in the decomposition-based multi-objective optimization algorithm in order to retain a better community structure. In network snapshots with time steps greater than 2, the Tchebycheff method is used to update the neighbors of each sub-problem as a way to eliminate the poorer community structure. For the neighbor x^j(j ∈ B(i)) of the i-th sub-problem $y_i^{{}^{\prime }}$ after reorganization, if ${\mathrm{g}}^{te}\left(y_i^{{}^{\prime }}|{\lambda }^j,z\right)\le {\mathrm{g}}^{te}\left(x^j|{\lambda }^j,z\right)$, then the solution of that sub-problem will replace its neighbor. Then all solutions of $y_i^{{}^{\prime }}$ are added to the external population EP, and all dominated solutions are removed. Finally, the maximum modularity Q and NMI values are found and the reference point $Z\ast =new\left\{NM{I}_{\max },Q_{\max }\right\}$ is updated.

Rationalized selection strategy

The result we obtain in the partitioning process of dynamic network communities based on multi-objective optimization is a set of Pareto optimal solution sets consisting of multiple non-dominated solutions. We know that in the process of dynamic network community discovery, we need to choose a non-dominated solution from the optimal set of solutions at the current time step as the final pole at the current moment and also as a reference for the next moment. Most researchers choose the community structure with the highest snapshot quality, but this can easily lose the community structure with better time cost. Therefore, in order to accommodate such a situation, we propose to lose less snapshot quality in exchange for a better time cost. First we find the solution with the highest snapshot quality, then we use that solution as a reference to find other solutions in the range of 0.01 from it, and finally we select the solution with the best time cost from all the solutions in that range as the final solution.

SYN-FIX datasets

The SYN-FIX datasets has 128 nodes in the network containing four communities, each community has 32 nodes at each time step, the average degree of each network is set to 16, and the edge connecting each node to a node that is not in the same community is zout. The size of zout determines the fuzziness of the network, and the larger the zout value the fuzzier the community structure. To simulate the evolutionary properties of dynamic networks, three nodes are randomly selected from each community to move to other communities at each time step, starting from the second time step. In this data, experiments were performed for two dynamic networks with zout = 3 and zout = 5, respectively, with a time step of 10.

As the results in Fig. 3 show, on the synthetic datasets SYN-FIX, algorithm NOME has the same NMI value of size 1 in Fig. 3A zout = 3 as algorithms ECD and ePMCL at all moments, both outperforming DYNMOGA and FacetNet. In Fig. 3B zout = 5, algorithms NOME and ePMCL have NMI values of 1 at all moments. The specific values for the five algorithms are given in Table 2, where DYNMOGA and ECD have NMI values of 1 for only two moments. Therefore NOME and ePMCL outperform other algorithms overall on this datasets.

SYN-VAR datasets

The network of the SYN-VAR datasets consists of 256 nodes divided into four communities, each consisting of 64 nodes. Compared with the SYN-FIX datasets, SYN-VAR simulates the creation and disappearance of nodes, while the meaning of the parameter zout in SYN-VAR is the same as that of SYN-FIX. Experiments are also conducted in the SYN-VAR datasets for two dynamic networks with zout = 3 and zout = 5, respectively. In order to generate a dynamic network with a time step of 10, eight nodes are randomly selected from the four communities in each of the first five time steps to synthesize a new community consisting of 32 nodes as the newly generated nodes and communities in the next time step. In the subsequent five time steps, the nodes will return to their original communities, so the number of communities in this dynamic network with 10 time steps is 4, 5, 6, 7, 8, 8, 7, 6, 5 and 4. In addition, the average degree of each node in the community is set to half the community size, and at each time step, 16 nodes are randomly removed from the network, while 16 new nodes are added to the network.

Figure 4 shows the results of the five algorithms on the synthetic datasets SYN-VAR. In Fig. 4A zout = 3, algorithm ePMCL has slightly higher NMI values than NOME and ECD at moments 3 and 5, and all three are better than the other two algorithms overall. In Fig. 4B zout = 5, combined with Table 3, it is clear that NOME has the optimal NMI value for eight moments, ECD has the optimal NMI value for four moments, and ePMCL has the optimal NMI value for nine moments. Therefore, the ePMCL and NOME algorithms outperform the other algorithms overall in the synthetic datasets SYN-VAR.

Cellphone calls datasets

The Cellphone Calls datasets, developed by VAST2008 Mini-Challenge3, consists of cellphone call records between members of the Paraiso movement, which cover 10 days of calls from 400 cellphones in June 2006. To build a dynamic network, in this datasets, a snapshot of the network at each time step consists of all phone call records during a day, where a node represents a phone and an edge between nodes represents a call between two phones that are available. Thus this data contains 10 network snapshots.

Enron mail datasets

The Enron Mail datasets is derived from a collection of emails from 1999 to 2002 from a U.S. company called Enron Corporation. This experiment only selected email data from 2001 with potentially unusual emails, which contained 252,756 emails from 151 employees. It is divided by month into 12 copies representing 12 network snapshots, each snapshot consists of 151 nodes representing 151 employees. Also the node-to-node edges represent users with emails to and from each other.

Table 4 shows the mean and standard deviation of NMI and modularity returned by NOME, DECS, and ECD on the real dataset Cellphone Calls. It can be found that the mean value of NMI and the mean value of modularity obtained by our proposed algorithm during the evolution process on continuous time steps are much higher than other algorithms. In addition, the smaller the standard deviation, the more stable the algorithm. The standard deviation values of our proposed algorithm are mostly better than other algorithms, among which NMI accounts for seven out of 10-time scales and modularity accounts for 6. Therefore, our proposed method has better accuracy and stability than other community detection algorithms on this dataset.

Figure 5 shows the average NMI values obtained for NOME, DECS, ePMCL and ECD on Fig. 5A Cellphone Calls and Fig. 5B Enron Mail datasets. As we observed, our algorithm NOME has the highest NMI values at each time step on both real datasets, while ePMCL performs the worst. The gap between the remaining ECD and DECS algorithms is not very big, but there is still a certain gap compared with our proposed algorithm. The results show that NOME can significantly improve the similarity between the time steps of the community structure after using the strategy of rationalization rules, and also reveal that the time cost of the evolving community structure on the Cellphone Calls and Enron Mail datasets is better than other algorithms.

Tables 5 and 6, respectively, show the mean, standard deviation, and difference of NMI and the mean, standard deviation, and difference of modularity Q returned by the three algorithms NOME, DECS, and ECD on the real data set Enron Mail. Among them, DE-—NMI— represents the difference between NOME and DECS, EC-—NMI— represents the difference between NOME and ECD; DE-—Q— represents the difference between NOME and DECS, and EC-—Q— represents the difference between NOME and ECD. To further compare the pros and cons of the three different algorithms on the Enron dataset, we integrated the difference between their NMI and Q into Table 7 for comparison.

Figure 6 shows the modularity values of NOME, DECS and ECD on the real datasets. From Fig. 6A, we can see that the modularity values of the NOME algorithm are much higher than the other two algorithms at each moment on the Cellphone Calls datasets, so the performance of the NOME algorithm on the Cellphone Calls datasets is much better than the other algorithms. However, in Fig. 6B the NOME algorithm has slightly lower modularity values than the other two algorithms on the Enron Mail datasets for all moments except the first moment. To better evaluate the effectiveness of the algorithms, we compare the three algorithms on the Enron Mail datasets in terms of the difference between NMI and modularity, and we expect a smaller modularity in exchange for a larger boost in NMI on this datasets.

From Table 7, we can see that in the comparison between NOME and DECS, there are nine best results for NMI and only two best results for modularity, and in the comparison between NOME and ECD, there are 10 best results for NMI and only one best result for modularity. Therefore, the results show that the strategy of border node occupancy distribution proposed by us can better assign nodes to the correct community, to improve the quality of community division. Although affected by the rationalized selection strategy, some sacrifices may have better community results, but overall it can be concluded that the NOME algorithm has better accuracy in dynamic community discovery than other advanced algorithms.