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

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

Marco Corneli Pierre Latouche Fabrice Rossi Université Paris 1 Panthéon-Sorbonne - Laboratoire SAMM
90 rue de Tolbiac, F-75634 Paris Cedex 13 - France
Abstract

We develop a model in which interactions between nodes of a dynamic network are counted by non homogeneous Poisson processes. In a block modelling perspective, nodes belong to hidden clusters (whose number is unknown) and the intensity functions of the counting processes only depend on the clusters of nodes. In order to make inference tractable we move to discrete time by partitioning the entire time horizon in which interactions are observed in fixed-length time sub-intervals. First, we derive an exact integrated classification likelihood criterion and maximize it relying on a greedy search approach. This allows to estimate the memberships to clusters and the number of clusters simultaneously. Then a maximum-likelihood estimator is developed to estimate non parametrically the integrated intensities. We discuss the over-fitting problems of the model and propose a regularized version solving these issues. Experiments on real and simulated data are carried out in order to assess the proposed methodology.

keywords:
Dynamic network, Stochastic block model, exact ICL, Non homogeneous Poisson Process.

1 Introduction

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 multi-graph (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 take 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 chose 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 Section 2 we mention works related to the approach we propose, Section 3 presents the proposed temporal extension of the SBM, Section 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.

2 Related Works

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 given its cluster at time . Xu and Hero III (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 (non-parametrically) the instantaneous intensity functions of the Poisson processes, they develop a variational EM algorithm to maximize an approximation of the likelihood.

3 The model

We consider a fixed set of nodes, , that can interact as frequently as wanted during the time interval . Interactions are directed from one node to another and are assumed to be instantaneous111In practice, the starting time of an interaction with a duration will be considered.. A natural mathematical model for this type of interactions is provided by counting processes on . Indeed a counting process is a stochastic process with values that are non negative integers increasing through time: the value at time can be seen as the number of interactions that took place from to . Then the classical adjacency matrix of static graphs is replaced by a collection of counting processes, , where is the counting process that gives the number of interactions from node to node . We still call 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).

3.1 Non-homogeneous Poisson counting process

We first chose a simple form for : we assume that this process is a non-homogeneous Poisson counting process (NHPP) with instantaneous intensity given by the function from to , . For , it then holds

(1)

where is the (non negative) number of interactions from to that took place during . (We assume that .)

3.2 Block modelling

The main idea of the SBM (Holland et al., 1983; Wang and Wong, 1987) is to assume that nodes have some (hidden) characteristics that solely explain their interactions, in a stochastic sense. In our context this means that rather than having pairwise intensity functions , those functions are shared by nodes that have the same characteristics.

In more technical terms, we assume the nodes are grouped in clusters () and introduce a hidden cluster membership random vector such that

The random component is assumed to follow a multinomial distribution with parameter vector such that

In addition, the are assumed to be independent (knowing ) and thus

(2)

where denotes the cardinal of . Notice that this part of the model is exactly identical to what is done in the classical SBM.

In a second step, we assume that given , the counting processes are independent and in addition that the intensity function depends only on and . In order to keep notations tight we denote the common intensity function and we will not use directly the pairwise intensity functions . We denote the matrix valued intensity function .

Combining all the assumptions, we have for

(3)

3.3 Discrete time version

In order to make inference tractable, we move from the continuous time model to a discrete time one. This is done via a partition of the interval based on a set of instants

that defines intervals (with arbitrary length ). The purpose of the partition is to aggregate the interaction. Let us denote

(4)

In words, measures the increment, over the time interval , of the Poisson process counting interactions from to . We denote by the random vector

Thanks to the independence of the increments of a Poisson process, we get the following joint density:

(5)

The variations of inside an interval have no effect on the distribution of . This allows us to use the integrated intensity function defined on by

In addition, we denote by the increment of the integrated intensity function over

Then equation (5) becomes

(6)

with .

Using the block model assumptions, we have in addition

(7)

where we have used the fact that (which leads to , etc.).

Considering the network as a whole, we can introduce two tensors of order 3. is a random tensor whose element is the random variable and is the tensor whose element is . can be seen as an aggregated (or discrete time version) of the adjacency process while can be seen as summary of the matrix valued intensity function .

The conditional independence assumption of the block model leads to

(8)

To simplify the rest of the paper, we will use the following notations

The joint distribution of , given z and , is

(9)

where

is the total number of interactions from cluster to cluster (possibly equal to ) and with

3.4 A constrained version

As will be shown in Section 4.4, the model presented thus far is prone to over fitting when the number of sub-intervals is large compared to . Additional constraints on the intensity functions are needed in this situation.

Let us consider a fixed pair of clusters . So far, the increments are allowed to differ on each over the considered partition. A constraint can be introduced by assigning the time intervals to different time clusters and assuming that increments are identical for all the intervals belonging to the same time cluster. Formally, we introduce clusters () and a hidden random vector , labelling memberships

Each is assume to follow a multinomial distribution depending on parameter

and in addition the are assumed to be independent, leading to

(10)

The random variable is now assumed to follow the conditional distribution

(11)

Notice that the new Poisson parameter replaces in the unconstrained version. The joint distribution of , given and , can easily be obtained

(12)

where

Remark 1.

The introduction of this hidden vector is not the only way to impose regularity constraints to the integrated function . For example, a segmentation constraint could be imposed by forcing each temporal cluster to contain only adjacent time intervals.

3.4.1 Summary

We have defined two generative models:

Model A

the model has two meta parameters, the number of clusters and the parameters of a multinomial distribution on . The hidden variable is generated by the multivariate multinomial distribution of equation (2). Then the model has a tensor of parameters . Given and , the model generates a tensor of interaction counts using equation (3.3).

Model B

is a constrained version of model A. In addition to the meta parameters and of model A, it has two meta parameters, the number of clusters of time sub-intervals and the parameters of a multinomial distribution on . The hidden variable is generated by the multivariate multinomial distribution of equation (10). Model B has a tensor of parameters . Given , and , the model generates a tensor of interaction counts using equation (12).

Unless specified otherwise “the model” is used for model A.

4 Estimation

4.1 Non parametric estimation of integrated intensities

In this Section we assume that is known. No hypothesis has been formulated about the shape of the functions and the increments of these functions over the partition introduced can be estimated by maximum likelihood (ML), thanks to equation (3.3)

where denotes those terms not depending on . It immediately follows

(13)

where denotes the ML estimator of . In words, can be estimated by ML as the total number of interactions on the sub-graph corresponding to the connections from cluster to cluster , over the time interval , divided by the number of nodes on this sub-graph. Once the tensor has been estimated, we have a point-wise, non parametric estimator of , for every , defined by

(14)

Thanks to the properties of the ML estimator, together with the linearity of (14), we know that is an unbiased and convergent estimator of .

Remark 2.

Estimator (14) at times , can be viewed as an extension to random graphs and mixture models of the non parametric estimator proposed in Leemis (1991). In that article, -trajectories of independent NHPPs, sharing the same intensity function, are observed and the proposed estimator is basically obtained via method of moments.

In all the experiments, we consider the following step-wise linear estimator of

(15)

which is a linear combination of estimators in equation (14) on the interval . This is a consistent and unbiased estimator of at times only.

When considering model B, equations (13) and (14) are replaced by

(16)
(17)

Equation (15) remains unchanged, but an important difference between the constrained model and the unconstrained one should be understood: in the former, each interval corresponds to a different slope for the function whereas in the latter we only have different slopes, one for each time cluster.

4.2 Icl

Since the vector , as well as the number of clusters are unknown, estimator (13) cannot be used directly. Hence we propose a two step procedure consisting in

  1. providing estimates of and ,

  2. using these estimates to implement (13) and (14).

To accomplish the first task, the same approach followed in Côme and Latouche (2015) is adopted: we directly maximize the the joint integrated log-likelihood of complete data (ICL), relying on a greedy search over the labels and number of clusters. To perform such a maximization, we need the ICL to have an explicit form. This can be achieved by introducing conjugated prior distributions on the model parameters. The ICL can be written as

(18)

This exact quantity is approximated by the well known ICL criterion (Biernacki et al., 2000). This criterion, obtained through Laplace and Stirling approximations of the joint density on the left hand side of equation (18), is used as a model selection tool, since it penalizes models with a high number of parameters. In the following, we refer to the joint log-density in equation (18) as to the exact ICL to differentiate it from the ICL criterion.

We are now going to study in detail the two quantities on the r.h.s. of the above equation. The first probability density is obtained by integrating out the parameter

In order to have an explicit formula for this term, we impose the following Gamma prior conjugated density over the tensor :

where the hyper-parameters of the Gamma prior distribution have been set constant to and for simplicity.222The model can easily be extended to the more general framework: By using the Bayes rule

we get:

which can be integrated with respect to to obtain

(19)

We now focus on the second density on the right hand side

A Dirichlet prior distribution can be attached to in order to get an explicit formula, in a similar fashion of what we did with :

The integrated density can be proven to reduce to

(20)

4.3 Model B

When considering the constrained framework described at the end of the previous section, the ICL is defined

and it is maximized to provide estimates of and . The first density on the right hand side is obtained by integrating out the hyper-parameter . This integration can be done explicitly by attaching to the following prior density function

The second integrated density on the right hand side can be read in (20) and the third is obtained by integrating out the parameter , whose prior density density function is assumed to be

The exact ICL is finally obtained by taking the logarithm of

(21)

4.4 Greedy search

By setting conjugated prior distributions over the model parameters, we obtained an ICL (equation (18)) in an explicit form. Nonetheless explicit formulas to maximize it, with respect to and , do not exist. We then rely on a greedy search algorithm, that has been used to maximize the exact ICL, in the context of a standard SBM, by Côme and Latouche (2015). This algorithm basically works as follows:

  1. An initial configuration for both and is set (standard clustering algorithms like k-means or hierarchical clustering can be used).

  2. Labels switches leading to the highest increase in the exact ICL are repeatedly made. A label switch consists in a merge of two clusters or in a node switch from one cluster to another.

Remark 3.

The greedy algorithm described in this section, makes the best choice locally. A convergence toward the global optimum in not guaranteed and often this optimum can only be approximated by a local optimum reached by the algorithm.

Remark 4.

The exact ICL (as well as the ICL criterion) penalizes the number of parameters. Since the tensor has dimension , when , which is fixed, is very hight, the ICL will take its maximum for . In other words the only way the ICL has to make the model more parsimonious is to reduce up to one. By doing so, any community (or other) structure will not be detected. This over-fitting problem has nothing to see with the possible limitations of the greedy search algorithm and it can be solved by switching to model B.

Once has been fixed, together with an initial value of , a shuffled sequence of all the nodes in the graph is created. Each node in the sequence is moved to the cluster leading to the highest increase in the ICL, if any. This procedure is repeated until no further increase in the ICL is still possible. Henceforth, we refer to this step as to Greedy-Exchange (GE). When maximizing the modularity score to detect communities, the GE usually is a final refinement step to be adopted after repeatedly merging clusters of nodes. In that context, moreover, the number of clusters is initialized to and each node is alone in its own cluster. See for example Noack and Rotta (2008). Here, we follow a different approach, proposed by Côme and Latouche (2015) and Blondel et al. (2008): after running the GE , we try to merge the remaining clusters of nodes in the attempt to increase the ICL. In this final step (henceforth GM), all the possible merges are tested and the best one is retained.

The ICL does not have to be computed before and after each swap/merge: possible increases can be assessed directly. When switching one node (say ) from cluster to , with , the change in the ICL is given by333Hereafter, the “*” notation refers to the statistics after switching/merging.

The only statistics not simplifying, are those involving and , hence the equation above can be read as follows

(22)

where is the term inside the product on the right hand side of equation (19) and and refer to new configuration where the node in in .

When merging clusters and into the cluster , the change in the ICL can be expressed as follows:

(23)

When working with model B, we need to initialize and . Then a shuffled sequence of time intervals is considered and each interval is swapped to the time cluster leading to the highest increase in the ICL (GE for time intervals). When no further increase in the ICL is possible, we look for possible merges between time clusters in the attempt to increase the ICL (GM for time intervals). Formulas to directly assess the increase in the ICL can be obtained, similar to those for nodes swaps and merges. In case of model B, different strategies are possible to optimize the ICL:

  1. GE + GM for nodes at first and then for times (we will call this strategy TN, henceforth).

  2. GE + GM for time intervals at first and then for nodes (NT strategy).

  3. An hybrid strategy, involving alternate switching of nodes and time intervals (M strategy).

We will provide details about the chosen strategy case by case in the following.

5 Experiments

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 Section 4.4 being implemented in C++. A Euclidean hierarchical clustering algorithm was used to initialize the labels and was set to .

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 of labels in situations where the standard SBM fails. Thus, we simulate interactions between 50 nodes, grouped in two hidden clusters and , over 100 time intervals of unitary length. The generative model considered for the simulations depends on two time clusters and containing a certain number of time intervals . If is in then is drawn from a Poisson distribution . Otherwise, is drawn from a Poisson distribution . The matrices and are given by

where is a free parameter in . 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 , the stronger the contrast between the interactions pattern inside and outside the cluster. In this paragraph, is set equal to 2 and the proportions of the clusters are set equal (). The number of time intervals assigned to each time cluster is assumed to be equal to . In the following, we consider

This generative model defines two integrated intensity functions (IIFs), say and . 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 can be observed in Figure 1.

\thesubsubfigure
\thesubsubfigure
Figure 1: Real 1 and estimated 1 integrated intensity functions (IIFs) according to the considered generative model (). In blue we have , for , in red .

A tensor , with dimensions , is drawn. Its component is the sampled number of interactions from node to node over the time interval . 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 (providing estimates in 11.86 seconds on average) the latter on the corresponding adjacency matrix. This experiment was repeated 50 times and estimates of random vector were provided at each iteration. Each estimate is compared with the true 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 , the two integrated intensity functions can be estimated through the estimator in equation (15). Results can be observed in Figure 1, where the estimated functions (coloured dots) overlap the real functions 1.

Over fitting

We now illustrate how the model discussed so far fails in recovering the true vector 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 :

Despite the lower contrast (from to in and ), with and time sub-intervals of unitary length, the TSBM model still always recovers the true vector . Now we consider a finer partition of by setting and as well as scaling the intensity matrices as follows

Moreover, we set

and is the complement of , as previously. Finally, we sampled 50 dynamic graphs over the interval from the corresponding generative model. Thus, each graph is characterized by a sampled tensor .

Unfortunately, the model is not robust to such changes. Indeed, when running the greedy ICL algorithm on each sampled tensor , 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 paragraph 4.4, the ICL penalizes the number of parameters and since the tensor has dimension , for a fixed , when moving from the larger decomposition () to the finer one (), the number of free parameters in the model is approximatively444The dimension of the vector does not change. multiplied by 10. The increase we observe in the likelihood, when increasing the number of clusters of nodes from to , 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 and a single cluster is detected.

Model B allows to tackle this issue. When allowing the integrated intensity functions and to grow at the same rate on each interval belonging to the same time cluster , we basically reduce the third dimension of the tensor from to .

The greedy ICL algorithm for Model B was run on each sampled tensor , providing estimates of and in minutes, on average. A hierarchical clustering algorithm was used to initialize the time labels , and the initial number of time clusters was set to . 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 Figure 2.

Figure 2: Box plots for both clusterings of nodes and time intervals: 50 dynamic graphs were sampled according to the considered generative model, estimates of and are provided by the greedy ICL (model B).

Note that the results were obtained through the optimization strategy TN. The other two strategies described in section 4.4, namely the NT strategy and the M strategy, led to similar results in terms of final ICL and ARIs.

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 . The previous scenario shows that the TSBM can recover the true vector 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 nodes and 50 time intervals

These time intervals are grouped in two time clusters and , the former containing the first 25 time intervals, the latter the last 25 time intervals. If is in then is drawn from a Poisson distribution . Otherwise, is drawn from a Poisson distribution . The matrix is given by

and is a free parameter in . Hence, we have two different integrated intensity functions, say and with the same roles as in the previous section. These two functions are plotted in Figure 3, for a value of .

\thesubsubfigure
\thesubsubfigure
Figure 3: Real 3 and estimated 3 integrated intensity functions (IIFs) according to the considered generative model. In blue we have , for , in red .

We investigated six values for the parameter

For each value of , we sampled 50 tensors , of dimension , according to the generative model considered. Interactions are aggregated over the time interval 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 were obtained in a mean running time of 5.52 seconds. 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 Figure 4.

\thesubsubfigure ARIs obtained by greedy ICL.
\thesubsubfigure ARIs obtained with the Gibbs sampling procedure for SBM.
Figure 4: Box plots of ARIs for different levels of contrast (). We compare the proposed model with a standard SBM.

Although the SBM leads to slightly better clustering results for small values of (2.2, 2.3) and the TSBM for higher values of (2.5, 2.6), we observe that the two models have quite similar performances (in terms of accuracy) in this scenario.

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

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 takes values in the set and 50 dynamic graphs were sampled, according to the described settings, for each value of . We ran the TSBM on each dynamic graph obtaining 50 estimates of the labels vector (one for each ) and box and whiskers plots for each group of ARIs can be seen in Figure 5.

Figure 5: Box plots of ARIs for different levels of contrast (). Data have been sampled by non-homogeneous Poisson processes counting interactions in a dynamic graph whose nodes are grouped in three clusters and interactivity patterns vary across two time clusters.

By comparing this figure with Figure 4, 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 we are using is related to the ML estimator, we can imagine a slower convergence to the true value of .

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 dataset used in this section was collected during the ACM Hypertext conference held in Turin, June 29th - July 1st 2009. We focus on the first conference day (24 hours) and consider a dynamic network with 113 nodes (conference attendees) and 96 time intervals (the consecutive quarter-hours in the period: 8am of June 29th - 7.59am 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 meters for a time period of at least 20 seconds555More informations about the way the data were collected can be found in Isella et al. (2011) or visiting the website http://www.sociopatterns.org/datasets/hypertext-2009-dynamic-contact-network/. . The data set we considered consists of several lines similar to the following one

ID ID Time Interval () Number of interactions
52 26 5 16

It means that conference attendees 52 and 26, between 9am and 9.15am, have spoken for .

We set and the vector 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 and were provided in 13.81 seconds, on average. The final values of the ICL can be observed as box plots in Figure 6 .

Figure 6: Box plot of the ten final values of the ICL produced by the greedy ICL algorithm for different initializations.

The estimates associated to the highest ICL correspond to 5 node clusters. In Figure 7, we focus on the cluster , containing 48 nodes. In Figure 6(a) we plotted the time cumulated interactions inside the cluster. As it can be seen the connectivity pattern for this cluster is very representative of the entire graph: between 13pm and 14pm and 18pm and 19.30pm there are significant increases in the interactions intensity.

(a) Cumulated aggregated connections inside cluster .
(b) Estimated IIF for interactions inside cluster .
Figure 7: in Figure 6(a), cumulated aggregated connections for each time interval for cluster . In Figure 6(b) the estimated IIF for interactions inside cluster . Vertical red lines delimit the lunch break and the wine and cheese reception.

The estimated integrated intensity function (IIF) for interactions inside this cluster can be observed in Figure 6(b). 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 gathering666More informations at http://www.ht2009.org/program.php.:

  • 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 dataset by using the optimization strategy M described at the end of Section 4.4 (other strategies lead in this case to similar results) and was set equal to 20. The time clustering provided by the greedy ICL algorithm can be observed in Figure 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/color. Two important things should be noticed:

  • 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 to a certain intensity level, so time intervals in the same cluster, not necessarily adjacent, share the same global interactivity pattern.

  • 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.

(a) Aggregated connections.
(b) Clustered time intervals.
Figure 8: in Figure 7(a), aggregated connections for each time interval for the whole network. In Figure 7(b) interactions of the same form/color take place on time intervals assigned to the same cluster (model B).

6 Conclusion

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 non parametric 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.

Appendix A Computational complexity

In this section we provide details about the computational complexity of the main model presented in this paper, namely the model A. Assuming that the gamma function can be computed in constant time (see Press et al., 2007), we focus on the three statistics appearing in equation (3.3), namely

  1. ,

  2. ,

  3. .

The whole computation task consists in evaluating the increase in ICL induced by nodes exchanges and merges. Those computations involves the tree quantities listed above. The tensor is stored in a three dimensional array, never resized, occupying a memory space. Hence, at any time during the algorithm its elements can be accessed and modified in constant time. The tensor is handled similarly and clusters sizes (we recall that corresponds to the size of cluster ) are also stored in arrays. In order to evaluate the ICL changes, induced by an operation, we need to maintain aggregated interaction counts for each node: for a node we have, e.g.

the number of interactions from node to cluster inside the time interval . Similarly

denotes the number of interactions from cluster to node inside the time interval . Other related quantities are considered. These structures occupy a memory space of .

Exchanges

In order to evaluate the ICL increase induced by the switch of a node (say ) from cluster to cluster , we perform the following operations:

  • (respectively ) is reduced by () and () is increased by the same amount;

  • (respectively ) is reduced by () and () is increased by the same amount;

  • () is reduced (increased) by one.

Although these operations are in constant time, they are involved in a sum with elements (this can be seen in equation (22)), so that the total cost of the test is . Since node can be switched to remaining clusters and the graph has nodes, the cost of a full exchange routine is .

Remark 6.

When a node is actually switched from its cluster to another one, all data structures are updated but the update cost is dominated by the cost of the testing phase described above.

Notice that we have evaluated the total cost of one full exchange routine, i.e., in the case where all nodes are considered once. Reductions in the number of clusters (very likely to be induced by exchanges in case is high) are not taken into account.

Merges

The entire merge routine, consisting in a test phase and an actual merge, has a computational cost that is dominated by the cost of exchanges. Consider a cluster . We first look for the cluster (say ) leading to the best merge (highest increase in the ICL) with . This operation has a cost of : for each the evaluation of the increase in ICL has a cost of (see equation (23)) and