Block modelling in dynamic networks with non-homogeneous Poisson processes and exact ICL

Graph clustering (Schaeffer 2007) is probably one of the main exploratory tools used in network analysis as it provides data analysts with a high level summarized view of complex networks. One of the main paradigms for graph clustering is community search (Fortunato 2010): a community is a subset of nodes in a graph that are densely connected and have relatively few connections to nodes outside of the community. While this paradigm is very successful in many applications, it suffers from a main limitation: it cannot be used to detect other important structures that arise in graphs, such as bipartite structures, hubs, authorities and other patterns.

The alternative solution favoured in this paper is provided by block models (Lorrain and White 1971; White et al. 1976): in such a model, a cluster consists of nodes that share the same connectivity patterns to other clusters, regardless of the pattern itself (community, hub, bipartite, etc.). A popular probabilistic view on block models is provided by the stochastic block model (SBM, Holland et al. 1983; Wang and Wong 1987). The main idea is to assume that a hidden random variable is attached to each node. This variable contains the cluster membership information while connection probabilities between clusters are handled by the parameters of the model. The reader is send to Goldenberg et al. (2009) for a survey of probabilistic models for graphs and to Wasserman and Faust (1994), Ch.16, for an overview of the stochastic block models.

This paper focuses on dynamic graphs in the following sense: we assume that nodes of the graph are fixed and that interactions between them are directed and take place at a specific instant. In other words, we consider a directed multigraph (two nodes can be connected by more than one edge) in which each directed edge is labelled with an occurrence time. We are interested in extending the SBM to this type of graphs. More precisely, the proposed model is based on a counting process point of view of the interactions between nodes: we assume that the number of interactions between two nodes follows a non-homogeneous Poisson counting process (NHPP). As in a standard SBM, nodes are assumed to belong to clusters that do not change over time; thus, the temporal aspect is handled only via the non-homogeneity of the counting processes. Then, the block model hypothesis takes the following form: the intensity of the NHPP that counts interactions between two nodes depends only on the clusters of the nodes. In order to obtain a tractable inference, a segmentation of the time interval under study is introduced and the interactions are aggregated over the sub-intervals of the partition. Following Côme and Latouche (2015), the model is adjusted to the data via the maximization of the integrated classification likelihood (ICL Biernacki et al. 2000) in an exact form. As in Côme and Latouche (2015) (and Wyse et al. (2014) for latent block models), the maximization is done via a greedy search. This allows us to choose automatically the number of clusters in the block model.

When the number of sub-intervals is large, the model can suffer from a form of over-fitting as the ICL penalizes only a large number of clusters. Therefore, we introduce a variant, based on the model developed in Corneli et al. (2015), in which sub-intervals are clustered into classes of homogeneous intensities. Those clusters are accounted for in a new version of the ICL which prevents over-fitting.

This paper is structured as follows: in Sect. 2, we mention works related to the approach we propose, Sect. 3 presents the proposed temporal extension of the SBM and Sect. 4 derives the exact ICL for this model and presents the greedy search algorithm used to maximize the ICL. Section 5 gathers experimental results on simulated data and on real-world data. Section 6 concludes the paper.

Numerous extensions of the original SBM have already been proposed to deal with dynamic graphs. In this context, both nodes memberships to a cluster and interactions between nodes can be seen as stochastic processes. In Yang et al. (2011), for instance, authors introduce a Markov Chain to obtain the cluster of node in t given its cluster at time $t-1$. Xu and Hero (2013) as well as Xing et al. (2010) used a state space model to describe temporal changes at the level of the connectivity pattern. In the latter, the authors developed a method to retrieve overlapping clusters through time. In general, the proposed temporal variations of the SBM share a similar approach: the data set consists in a sequence of graphs rather than the more general structure we assume. Some papers remove those assumptions by considering continuous time models in which edges occur at specific instants (for instance, when someone sends an email). This is the case of e.g. Dubois et al. (2013) and of Guigourès et al. (2012, 2015). A temporal stochastic block model, related to the one presented in this paper, is independently developed by Matias et al. (2015). They assume that nodes in a network belong to clusters whose composition do not change over time and interactions are counted by a non-homogeneous Poisson process whose intensity only depends on the nodes clusters. In order to estimate (nonparametrically) the instantaneous intensity functions of the Poisson processes, they develop a variational EM algorithm to maximize an approximation of the likelihood.

We consider a fixed set of N nodes, $\{1,\ldots ,N\}$, that can interact as frequently as wanted during the time interval [0, T]. Interactions are directed from one node to another and are assumed to be instantaneous.¹ A natural mathematical model for this type of interactions is provided by counting processes on [0, T]. Indeed a counting process is a stochastic process with values that are non-negative integers increasing through time: the value at time t can be seen as the number of interactions that took place from 0 to t. Then, the classical adjacency matrix $(X_{ij})_{1\le i, j\le N}$ of static graphs is replaced by a $N\times N$ collection of counting processes, $(X_{ij}(t))_{1\le i, j\le N}$, where $X_{ij}(t)$ is the counting process that gives the number of interactions from node i to node j. We still call ${\mathbf {X}}=(X_{ij}(t))_{1\le i, j\le N}$ the adjacency matrix of this dynamical graph.

We introduce in this section a generative model for adjacency matrices of dynamical graphs that is inspired by the classical stochastic block model (SBM).

In this section, experiments on both synthetic and real data are provided. All running times are measured on a twelve cores Intel Xeon server with 92 GB of main memory running a GNU Linux operating system, the greedy algorithm described in Sect. 4.4 being implemented in C++. A Euclidean hierarchical clustering algorithm was used to initialize the labels, and $K_{\max }$ was set to N/2.

In the following, we call TSBM the temporal SBM we propose and we refer to the optimization algorithm described in the previous section as greedy ICL.

5.1 Simulated data

5.1.1 First scenario

We start by investigating how the proposed approach can be used to efficiently estimate the vector ${\mathbf {z}}$ of labels in situations where the standard SBM fails. Thus, we simulate interactions between 50 (N) nodes, grouped in two hidden clusters ${\mathcal {A}}_1$ and ${\mathcal {A}}_2$, over 100 (U) time intervals of unitary length. The generative model considered for the simulations depends on two time clusters ${\mathcal {C}}_1$ and ${\mathcal {C}}_2$ containing a certain number of time intervals $I_1,\ldots I_U$. If $I_u$ is in ${\mathcal {C}}_1$, then $Y_{ij}^{I_u}$ is drawn from a Poisson distribution ${\mathcal {P}}(P_{z_i z_j})$. Otherwise, $Y_{ij}^{I_u}$ is drawn from a Poisson distribution ${\mathcal {P}}(Q_{z_i z_j})$. The matrices P and Q are given by

$$\begin{aligned} P= \begin{pmatrix} \psi &{} 1 \\ 1 &{} \psi \\ \end{pmatrix} \quad {\text {and}}\quad Q= \begin{pmatrix} 1 &{} \psi \\ \psi &{} 1 \\ \end{pmatrix}, \end{aligned}$$

where $\psi$ is a free parameter in $[1, \infty )$. When this parameter is equal to 1, we are in a degenerate case and there is not any structure to detect: all the nodes are placed in the same, unique cluster. The higher $\psi$, the stronger the $contrast$ between the interactions pattern inside and outside the cluster. In this paragraph, $\psi$ is set equal to 2 and the proportions of the clusters are set equal ($\varvec{\omega }=(1/2, 1/2)$). The number of time intervals assigned to each time cluster is assumed to be equal to U/2. In the following, we consider

$$\begin{aligned} {\mathcal {C}}_1&:=\{I_1, \ldots , I_{25}\}\cup \{I_{51},\ldots ,I_{75}\}, \\ {\mathcal {C}}_2&:=\{I_{26}, \ldots , I_{50}\}\cup \{I_{76},\ldots ,I_{100}\}. \end{aligned}$$

This generative model defines two integrated intensity functions (IIFs), say $\varLambda _1(t)$ and $\varLambda _2(t)$. The former is the IIF of the Poisson processes counting interactions between nodes sharing the same cluster, the latter is the IIF of the Poisson processes counting interactions between vertices in different clusters. These IIFs are shown in Fig. 1a.

Open image in new window

A tensor Y, with dimensions $N \times N \times U$, is drawn. Its (i, j, u) component is the sampled number of interactions from node i to node j over the time interval $I_u$. Moreover, sampled interactions are aggregated over the whole time horizon to obtain an adjacency matrix. In other words, each tensor is integrated over its third dimension. We compared the greedy ICL algorithm with the Gibbs sampling approach introduced by Nouedoui and Latouche (2013). The former was run on the tensor Y (providing estimates in 11.86 s on average) the latter on the corresponding adjacency matrix. This experiment was repeated 50 times, and estimates of random vector ${\mathbf {z}}$ were provided at each iteration. Each estimate $\hat{{\mathbf {z}}}$ is compared with the true ${\mathbf {z}}$ and an adjusted rand index (ARI Rand 1971) is computed. This index takes values between zero and one, where one corresponds to the perfect clustering (up to label switching).

Remark 5

The true structure is always recovered by the TSBM: 50 unitary values of the ARI are obtained. Conversely, the standard SBM never succeeds in recovering any hidden structures present in the data (50 null ARIs are obtained). This can easily be explained since the time clusters have opposite interaction patterns, making them hard to uncover when aggregating over time.

Relying on an efficient estimate of ${\mathbf {z}}$, the two integrated intensity functions can be estimated through the estimator in Eq. (15). Results are shown in Fig. 1b, where the estimated functions (coloured dots) overlap the real functions 1a.

Over-fitting We now illustrate how the model discussed so far fails in recovering the true vector ${\mathbf {z}}$ when the number of time intervals (and hence of free parameters) grows. We consider the same generative model of the previous paragraph, with a lower $\psi$:

$$\begin{aligned} P= \begin{pmatrix} 1.4 &{} 1 \\ 1 &{} 1.4 \\ \end{pmatrix} \quad {\text {and}}\quad Q= \begin{pmatrix} 1 &{} 1.4 \\ 1.4 &{} 1 \\ \end{pmatrix}. \end{aligned}$$

Despite the lower contrast (from 2 to 1.4 in P and Q), with $U=100$ and time sub-intervals of unitary length, the TSBM model still always recovers the true vector ${\mathbf {z}}$. Now, we consider a finer partition of [0, 100] by setting $U=1000$ and $\varDelta _u=0.1$ as well as scaling the intensity matrices as follows

$$\begin{aligned} {\tilde{P}}:= \begin{pmatrix} 0.14 &{} 0.1 \\ 0.1 &{} 0.14 \\ \end{pmatrix} \quad {\text {and}}\quad {\tilde{Q}}= \begin{pmatrix} 0.1 &{} 0.14 \\ 0.14 &{} 0.1 \\ \end{pmatrix}. \end{aligned}$$

Moreover, we set

$$\begin{aligned} {\mathcal {C}}_1:=\{I_1, \ldots , I_{250}\} \cup \{I_{501}, \ldots , I_{750}\} \end{aligned}$$

and ${\mathcal {C}}_2$ is the complement of ${\mathcal {C}}_1$, as previously. Finally, we sampled 50 dynamic graphs over the interval [0, 100] from the corresponding generative model. Thus, each graph is characterized by a sampled tensor Y.

Unfortunately, the model is not robust to such changes. Indeed, when running the greedy ICL algorithm on each sampled tensor Y, the algorithm does not see any community structure and all nodes are placed in the same cluster. This leads to a null ARI, for each estimation. As mentioned in Sect. 4.4, the ICL penalizes the number of parameters and since the tensor $\pi$ has dimension $K \times K \times U$, for a fixed K, when moving from the larger decomposition ($U=100$) to the finer one ($U=1000$), the number of free parameters in the model is approximatively⁴ multiplied by 10. The increase we observe in the likelihood, when increasing the number of clusters of nodes from $K=1$ to $K=2$, is not sufficient to compensate the penalty due to the high number of parameters, and hence, the ICL decreases. Therefore, the maximum is taken for $K=1$ and a single cluster is detected.

Model B allows to tackle this issue. When allowing the integrated intensity functions $\varLambda _1(t)$ and $\varLambda _2(t)$ to grow at the same rate on each interval $I_u$ belonging to the same time cluster ${\mathcal {C}}_d$, we basically reduce the third dimension of the tensor $\pi$ from U to D.

The greedy ICL algorithm for Model B was run on each sampled tensor Y, providing estimates of ${\mathbf {z}}$ and ${\mathbf {y}}$ in 2.38 min, on average. A hierarchical clustering algorithm was used to initialize the time labels ${\mathbf {y}}$, and the initial number of time clusters was set to $D_{\max }=\sqrt{U}$. In an attempt to avoid convergence to local maxima, ten estimates are built for each tensor, and the estimate leading to the best ICL is finally retained. The adjusted rand index is used to evaluate the clustering, as previously, and the results are presented as box plots in Fig. 2. Note that the results were obtained through the optimization strategy TN. The other two strategies described in Sect. 4.4, namely the NT strategy and the M strategy, led to similar results in terms of final ICL and ARIs.

Open image in new window

5.1.2 Second scenario

Since the node clusters are fixed over time, the TSBM model can be seen as an alternative to a standard SBM to estimate the label vector ${\mathbf {z}}$. The previous scenario shows that the TSBM can recover the true vector ${\mathbf {z}}$ in situations where the SBM fails. In this paragraph, we show how the TSBM and the SBM can sometimes have similar performances.

We considered dynamic graphs with 50 (N) nodes and 50 (U) time intervals

$$\begin{aligned} I_1, \ldots , I_{50}. \end{aligned}$$

These time intervals are grouped in two time clusters ${\mathcal {C}}_1$ and ${\mathcal {C}}_2$, the former containing the first 25 time intervals, the latter the last 25 time intervals. If $I_u$ is in ${\mathcal {C}}_1$, then $Y_{ij}^{I_u}$ is drawn from a Poisson distribution ${\mathcal {P}}(P_{z_i z_j})$. Otherwise, $Y_{ij}^{I_u}$ is drawn from a Poisson distribution ${\mathcal {P}}(2P_{z_i z_j})$. The P matrix is given by

$$\begin{aligned} P= \begin{pmatrix} \psi &{} 2 \\ 2 &{} \psi \end{pmatrix} \end{aligned}$$

and $\psi$ is a free parameter in $[2, +\infty )$. Hence, we have two different integrated intensity functions, say $\varLambda _1(t)$ and $\varLambda _2(t)$ with the same roles as in the previous section. These two functions are plotted in Fig. 3a, for a value of $\psi =4$.

Open image in new window

We investigated six values for the parameter $\psi$

$$\begin{aligned} \{2.1, 2.2, 2.3, 2.4,2.5,2.6\}. \end{aligned}$$

For each value of $\psi$, we sampled 50 tensors Y, of dimension $(50 \times 50 \times 50)$, according to the generative model considered. Interactions are aggregated over the time interval [0, 50] to obtain adjacency matrices. We ran the greedy ICL algorithm on each tensor and the Gibbs sampling (SBM) algorithm on each adjacency matrix. For the greedy ICL algorithm, estimates of vector ${\mathbf {z}}$ were obtained in a mean running time of 5.52 s. As previously, to avoid convergence to local maxima, ten different estimates are built for each tensor, the one leading to the highest ICL being retained. The results are presented as box plots in Fig. 4. Although the SBM leads to slightly better clustering results for small values of $\psi$ (2.2, 2.3) and the TSBM for higher values of $\psi$ (2.5, 2.6), we observe that the two models have quite similar performances (in terms of accuracy) in this scenario.

Open image in new window

To provide some intuitions about the scalability (see next paragraph) of the proposed approach we repeated the previous experiment by setting $K=3$ clusters, corresponding to the following connectivity matrix:

$$\begin{aligned} P= \begin{pmatrix} \psi &{} 2 &{} 2\\ 2 &{} \psi &{} 2\\ 2 &{} 2 &{} \psi \end{pmatrix}. \end{aligned}$$

The assignment of the time intervals to the time clusters is unchanged as well as the connectivity pattern on each time cluster are unchanged. The contrast parameter $\psi$ takes values in the set $\{2, 2.5,2.10, \ldots , 2.8 \}$ and 50 dynamic graphs were sampled, according to the described settings, for each value of $\psi$. We ran the TSBM on each dynamic graph obtaining 50 estimates of the labels vector ${\mathbf {z}}$ (one for each $\psi$) and box and whiskers plots for each group of ARIs are shown in Fig. 5. By comparing this figure with Fig. 4a, we can see that the model needs a slight higher contrast to fully recover the true structure. Actually, when increasing the number of clusters without increasing the number of nodes, the size of each cluster decreases (on average) and since the estimator of ${\mathbf {z}}$ we are using is related to the ML estimator, we can imagine a slower convergence to the true value of ${\mathbf {z}}$.

Open image in new window

5.1.3 Scalability

A full scalability analysis of the proposed algorithm as well as the convergence properties of the proposed estimators are outside the scope of this paper. Nonetheless, in “Appendix” we provide details about the computational complexity of the greedy ICL algorithm. Future works could certainly be devoted to improve both the algorithm efficiency and scalability through the use of more sophisticated data structures.

5.2 Real data

The data set used in this section was collected during the ACM Hypertext conference held in Turin, from 29 June 2009 to 1 July 2009. We focus on the first conference day (24 h) and consider a dynamic network with 113 (N) nodes (conference attendees) and 96 (U) time intervals (the consecutive quarter-hours in the period: 8 a.m. of June 29th–7.59 a.m. of June 30th). The network edges are the proximity face to face interactions between the conference attendees. An interaction is monitored when two attendees are face to face, nearer than 1.5 m for a time period of at least 20 s.⁵ The data set we considered consists of several lines similar to the following one

ID1	ID2	Time interval (15 min)	Number of interactions
52	26	5	16

It means that conference attendees 52 and 26, between 9 a.m. and 9.15 a.m., have spoken for $16 \times 20\,\hbox {s} \approx 5\,\hbox {min}\,30\,\hbox {s}$.

We set $K_{\max }=20$ and the vector ${\mathbf {z}}$ was initialized randomly: each node was assigned to a cluster following a multinomial distribution. The greedy algorithm was run ten times on the considered dataset, each time with a different initialization and estimates of ${\mathbf {z}}$ and K were provided in 13.81 s, on average. The final values of the ICL can be observed as box plots in Fig. 6.

Open image in new window

The estimates associated with the highest ICL correspond to 5 node clusters. In Fig. 7, we focus on the cluster ${\mathcal {A}}_4$, containing 48 nodes. In Fig. 7a, we plotted the time cumulated interactions inside the cluster. As it can be seen that the connectivity pattern for this cluster is very representative of the entire graph: between 13 p.m. and 14 p.m. and 18 p.m. and 19.30 p.m. there are significant increases in the interactions intensity. The estimated integrated intensity function (IIF) for interactions inside this cluster is shown in Fig. 7b. The function has a higher slope on those time intervals where attendees in the cluster are more likely to have interactions. The vertical red lines delimit two important times of social gathering:⁶

Open image in new window

• 13.00–15.00—lunch break.

• 18.00–19.00—wine and cheese reception.

We conclude this section by illustrating how Model B can be used to assign time intervals on which interactions have similar intensity to the same time cluster. We run the greedy ICL algorithm for Model B on the data set by using the optimization strategy M described at the end of Sect. 4.4 (other strategies lead in this case to similar results) and $D_{\max }$ was set equal to 20. The time clustering provided by the greedy ICL algorithm is shown in Fig. 8. On the left-hand side, the aggregated interactions for each quarter-hour during the first day are reported. On the right-hand side, interactions taking place into those time intervals assigned to the same time cluster have the same form/colour. Two important things should be noticed:

1. 1.

   The obtained clustering seems meaningful: the three time intervals with the highest interactions level are placed in the same cluster (blue), apart from all the others. More in general, each cluster is associated with a certain intensity level, so time intervals in the same cluster, not necessarily adjacent, share the same global interactivity pattern.

    
2. 2.

   There are not constraints on the number of abruptly changes connected with these five time clusters. In other words, time clusters do not need to be adjacent and this is the real difference between the approach considered in this paper (time clustering) and a pure segmentation one.

    

Open image in new window

We proposed a non-stationary extension of the stochastic block model (SBM) allowing us to cluster nodes of a network is situations where the classical SBM fails. The approach we chose consists in partitioning the time interval over which interactions are studied into sub-intervals of fixed length. Those intervals provide aggregated interaction counts that are increments of non-homogeneous Poisson processes (NHPPs). In a SBM inspired perspective, nodes are clustered in such a way that aggregated interaction counts are homogeneous over clusters. We derived an exact integrated classification likelihood (ICL) for such a model and proposed to maximize it through a greedy search strategy. Finally, a nonparametric maximum likelihood estimator was developed to estimate the integrated intensity functions of the NHPPs counting interactions between nodes. The experiments we carried out on artificial and real-world networks highlight the capacity of the model to capture non-stationary structures in dynamic graphs.