Fast Change Point Detection on Dynamic Social Networks

1 Introduction

Dynamic network analysis is increasingly used in complex application domains ranging from social networks (Facebook network evolution [Leskovec and Sosiˇc, 2016]) to biological networks (protein-protein interaction [Shih and Parthasarathy, 2012]), from political science (United Nations General Assembly voting network [Voeten, 2012]) to communication networks ([Klimt and Yang, 2004; Asur et al., 2007]). Such dynamic networks are often represented using the snapshot model. Under this model, every network snapshot (represented by a graph) is defined at a logical timestamp. Two questions of fundamental importance are - (i) how does a network evolve? (ii) when does a network change significantly so as to raise suspicion that something fundamentally different is happening?

Various generative models [Peixoto and Rosvall, 2015; Zhang et al., 2016] have been proposed to address question (i) - to explain the evolution of a network. They study network evolution under certain generative models [Erd˝os and R´enyi, 1960; Karrer and Newman, 2011]. In reality, the generative model itself might change, as addressed in question (ii) above. Existing work [Akoglu et al., 2014; Ranshous et al., 2015] use complex methods to detect such changes. A drawback of these methods is that they are timeconsuming, and hence not scalable (in terms of both network size and number of snapshots). We seek to find an efficient and effective solution which can scale up both with network size and with the number of snapshots.

In this paper, we present a simple and efficient algorithm based on likelihood maximization to detect change points in dynamic networks under the snapshot model. We demonstrate the utility of our algorithm on both synthetic and real world networks drawn from political science (congressional voting, UN voting), and show that it outperforms two recent approaches (DeltaCon[Koutra et al., 2016], and LetoChange[Peel and Clauset, 2015]) along the axes of both quality and efficiency. Our work has the following contributions:

1. Our approach is general purpose - it can accommodate various snapshot generative models (see Table 1).
2. We model network evolution as a first order Markov process and consequently our algorithm accounts for the temporal dependency while computing the dissimilarity between snapshots.
3. Our algorithm is effective and efficient on both synthetic and real world networks and has constant memory overhead that can be tuned by a user controlled parameter.

2 RelatedWork

Ranshous et al. [Ranshous et al., 2015], and Akoglu et al. [Akoglu et al., 2014] recently survey network anomaly detection methods which share related goals. Closest to our work are those ideas that fall under the theme of "Event and Change Detection", defined as: given a network sequence, a dissimilarity scoring function, and a threshold, a change is flagged if the dissimilarity of two consecutive snapshots is above the threshold[Ranshous et al., 2015]. A key differentiator is that the work surveyed there is somewhat ad-hoc while we take a more principled approach and assume there is a latent generation model governing the network dynamics and seek to detect the change in the latent space. Moreover, we consider the temporal dependency across the snapshots while none of the works surveyed explicitly account for temporal dependency or the bias introduced therein. We also note that our work is orthogonal to methods that seek to identify the temporal scale resolution for networks[Caceres and Berger-Wolf, 2013].

DeltaCon [Koutra et al., 2016] uses a graph similaritybased [Berlingerio et al., 2012] approach to detect change points in dynamic networks. It derives the features of a snapshot based on sociological theories. And the feature similarity of each consecutive snapshot pair is calculated. That work is model agnostic (has no assumption on the generation model of networks), and is the state-of-the-art in terms of efficiency. We compare our algorithm against this approach.

Moreno and Neville [Moreno and Neville, 2013], Bridges et al. [Bridges et al., 2015] and Peel and Clauset [Peel and Clauset, 2015] develop network hypothesis testing based approaches. The advantage is that one can get a p-value of the test, which quantifies the confidence of the conclusion. However, these approaches have two shortcomings: firstly, they need to assume a specific generation model of the networks (mKPGM, GBTER and GHRG respectively); secondly, they are extremely slow, mostly due to the bootstrapping for pvalue calculation. La Fond et al.'s work [La Fond et al., 2014] can also generate a p-value. It is tested against DeltaCon without reporting running time and efficiency concern is also mentioned in the paper. These algorithms will not work in our setting where the detection is done real time under bounded memory constraints. We compare our model agnostic algorithm against [Peel and Clauset, 2015].

Finally, Li et al. [Li et al., 2016] propose an online algorithm, and consider temporal dependency. The problem they study is different from ours in that they study information diffusion on network with fixed structure and use continuous time. A recent work by Zhang et al. [Zhang et al., 2016] also studies the dynamic network in a Markov chain setting. They focus on community detection while we focus on change point detection.

3 ProblemFormulation

This paper studies how to detect the times at which the fundamental evolution mechanism of a dynamic network changes. We assume that there is some unknown underlying model that governs the generative process. Our change point detection algorithm is agnostic to this model. We assume that the observed network snapshots are samples that depend on some generative model and the previous snapshot. Networks have fluctuation across snapshots even when the generative model stays unchanged. Only when the generative model changes do we consider it a fundamental change. We represent the evolutionary process as a Markov Network (Figure 1).

In Figure 1, Mt is the network generation model at time t. It is a triad Mt = 〈Typet ,Θt, αt〉, where αt is the continuity parameter at time t, Typet specifies the model, while Θt represents the model parameters (Table 1 consists of some generative models we experiment on). Gt is the network (graph) observable at time t. We assume the number of vertices in Gt is fixed to be N for all snapshots (the union of all nodes is used when there is node addition/deletion, as in [Peel and Clauset, 2015]), so each Gt has 2(N2) possible configurations, and T is the total number of snapshots we observe. As per Figure 1 the configuration of the network at time t, Gt, depends on the generation model at time t, Mt (unobserved) and the network configuration at time t - 1, Gt-1 (observed). Hence the networks in the observed sequence are samples from a conditional distribution (samples are not independent). The continuity rate parameter αt controls the fraction of edges and non-edges that are retained from the previous snapshot, Gt-1. The network at time t is assumed to be generated in the following way: for each dyad, keep the connection status from time t-1 with probability 1-αt, and with probability αt, resample the connection according to the generation model at time t. Consequently, the smaller αt is, the more overlap between two consecutive snapshots there is. Note that two consecutive network configurations may differ substantially if αt > 0, even though the underlying generation model may be the same. Moreover, the changes of the generation model are assumed to be rare across the time span (Mt ̸=Mt+1 is a rare event).

Problem Definition Our goal is to efficiently find a set S ⊂ {2, . . . , T} such that t ∈ S ⇐⇒ Mt ̸= Mt-1, that is, to efficiently find all the time points at which the network generation model is different from the previous time point.

4 Methodology

Given the graphical formulation of the problem, exact inference is impossible since we do not know the underlying generative model, and our observations are stochastic. However, even without prior knowledge of the generative model, we can still design an approximate inference technique based on MCMC sampling theory.

The framework is straightforward, as mentioned in Section 2: we first extract a "feature vector" from each snapshot, then quantify the dissimilarity between consecutive snapshots, and flag out a change point when the dissimilarity score is above a threshold. We use the joint edge probability as the "feature vector" (Section 4.1), exploit Kolmogorov-Smirnov statistic, Kullback-Leibler divergence and Euclidean distance for dissimilarity measure (Section 4.2), and use a permutation test like approach to determine the threshold (Section 4.3).

4.1 Edge Probability Estimation

In this subsection, we discuss how to (approximately) estimate the joint distribution of the dyads 1 . We track the presence or absence of a small fixed number of dyads throughout the entire observed sequence of network snapshots. We break down the observation sequence into fixed-length windows, and for each window we infer the joint distribution of the dyads in our sample. We model each dyad to be a conditionally independent two-state Markov chain (Figure 2, we use α instead of αt in this section for brevity) given the sequence of generative models (Mt)t≥0 (this conditional independence assumption is satisfied for the choices of models in Table 1). Note that even for generative processes that may

1 we refer node pairs, which may or not be linked, as dyads

Figure 2: Two state Markov Chain

result in greater dependence among dyads (such as the configuration model), in many cases such dependence will be local, and if the number of dyads sampled is small, then these dyads will be spread out enough to be considered independent. Moreover, the conditional independence assumption significantly improves computational efficiency ([Hunter et al., 2012]). The marginal probabilities of these dyads can then be estimated using the observed samples within each time window.

We formalize the estimation procedure below. Given a network sequence G1, . . . , GT , we group the networks into sliding windows. We define Wt to be a subsequence of s consecutive observed networks ending at network Gt, so Wt ≡ (Gt-s+1,Gt-s+2, . . . ,Gt). We use equal sized sliding windows with a step size η, and we obtain a window sequence Ws,Ws+η, . . . ,Ws+(i-1)η, . . . ,WT . Non-overlapping window setting uses η = s. In each window i, we can estimate the joint edge distribution (for the selected dyads) P(e1, e2, . . . , ek | Ms+(i-1)η), where ej = 1 indicates an undirected edge between the j-th dyad, and k << (N2 ) is the number of dyads tracked. For each of the models in Table 1, the joint distribution can be factorized into P(e1, e2, . . . , ek | Ms+(i-1)η) = ΠjP(ej | Ms+(i-1)η). (conditional independence, see method description above)

We can view a dyad across time as a two state Markov chain, and the chain length is the window size. We call a dyad across time a chain in the following text for brevity. Let pj ≡ P(ej = 1 | Ms+(i-1)η), qj ≡ 1 - pj , and suppose we are interested in k chains.

1) Maximum Likelihood Estimator (MLE)

The joint probability of the chains is (Figure 2)

where N (j) 01 is the number of transitions from 0 to 1 for a chain (non-edge to edge for the dyad j within the window), ck stands for the combinatorial coefficients independent of α, ⃗ p. and ∑p,q∈{0,1} N (j) pq = s - 1 for all js. And hence the log-likelihood (omitting the coefficient ck) is:

MLE for a single chain First consider there is only one chain. Solving the zero-derivative Equations (4) and (5), leads to estimators of α, p. And the estimators indeed lead to a negative definite Hessian, and therefore is the MLE. Hence we have

MLE for multiple chains The MLE for multiple chains essentially involves solving a high degree polynomial, which in general does not have a closed form solution.

where 1-α is the continuity rate. If α = 0 then all snapshots are identical, which is uninteresting, so we have α ̸= 0 in Equation (5).

A common way to solve such maximization problems is to employ numerical methods such as gradient descent. The drawback of such an approach is that it can be computationally expensive with hundreds of dyads and windows. Therefore, we settle for an approximation of the MLE that empirically approximates numerical values well 2 . Intuitively, the estimator for α should depend on all the chains, but chains that spend more time in both states 0 and 1 provide more information about α than chains that spend most time in one state (the latter may be due to small α or to a value of pj far from 1/2). Since we can easily compute the MLE for α for a single chain, we estimate α with a weighted average of the MLEs from the individual chains, with chains that spend more time in both states being weighted more heavily. We then estimate each pj by the MLE for the jth chain, since the chains are conditionally independent given α. This results in the following estimators.

2 At significant level 0.05, two sample t-test shows the approximated values equal to the numerical values.

Empirically we find the exponent p = ∞ works best, which means ˆαapprox is a simple average of the αjMLE corresponding to the chains with maximal value of N (j) 0∗ N (j) 1∗ . The continuity rate describes the temporal dependency among networks, and can help us determine a proper window size.

Drawbacks of MLE Though MLEs are consistent in general, there is no guarantee of unbiasedness for these particular estimators with limited samples. Moreover, they have three random quantities (N01,N10,N1∗,N0∗ in (6) have 3 degrees of freedom for fixed s) and hence require more samples to estimate, making it prohibitive in practice.

2) Simplified Estimator

To overcome the drawbacks of MLE, we propose a simple estimator for the edge probability which is consistent and unbiased, has only one random quantity and therefore requires fewer samples. The simple estimator essentially estimates the edge frequency in each window. If we know changes happen rarely, and the process stays in equilibrium in most of the time, we can show the following estimator to be consistent and unbiased in equilibrium:

which is the proportion of snapshots in which the dyad j being an edge within the window.

Proposition 1 In equilibrium, ˆ pjeq is consistent as chain length (window size) increases.

By ergodic theorem [Givens and Hoeting, 2012], π(j) 1 almost surely =========== lims→∞N (j) 1 /s = lims→∞ ˆ pjeq, where the first equation means almost sure convergence, and implies convergence in probability (estimator being consistent).

Proposition 2 In equilibrium, ˆ pjeq is unbiased.

The above propositions imply that the larger the window size the better the estimation, and that in equilibrium, the temporal dependency (continuity rate) has no impact on estimating the onset probability of a Markov chain, and hence no impact on estimating the edge probability of a snapshot.

Although MLE is close to the true value when the chain is long enough (100 or longer), we do not use so large a window size in practice (20 usually, no larger than 50). Experiments (Figure 3) show that the simplified estimator is much better than MLE for change point detection in practice.

4.2 Distance Measure

Now, we need to compare the probability distributions of edges across consecutive windows. Kolmogorov-Smirnov (KS) statistic and Kullback-Leibler (KL) divergence are two common measures for comparing distribution. Their calculations require the enumeration of the whole state space and hence exponential to the number of variables for joint distributions. Although KS statistic is designed for univariate distribution, we can map the joint distribution, which has multivariate binary variables, to one dimension by decoding the binary vectors as an integer and use KS statistic. We bootstrap from empirical distributions of two consecutive windows respectively and use two sample KS test to quantify the difference of two distributions. We can use divide-and-conquer to alleviate the exponential complexity: partition the dyads into g groups, compute KL/KS dissimilarity within each small group, and record the median among all the g groups as the final dissimilarity.

Both of the above measures have good quality in terms of change point detection (Figure 4), but KS statistic is extremely slow (see extended version [Wang et al., 2017]), mostly due to the large sample bootstrap from each window. Euclidean distance, though lacks probabilistic interpretation, has linear complexity and has reasonable quality in practice.

4.3 Threshold Determination

Suppose we have w windows, then we compare w-1 pairs of distributions and get w - 1 difference/distance scores. How do we choose a threshold to determine at which window the network changes? We use a permutation test [Pitman, 1937] based approach to determine the threshold. For a desired significance level αs, we bootstrap from the w - 1 distance scores, and use the upper 100αs% quantile as the threshold.

4.4 Complexity Analysis

The algorithm is linear to the number of windows and constant to the network size for moderately large network. Only a small fraction of dyads in the network is sampled and tracked. The sampling of the dyads is only performed once at the beginning, and hence irrelevant to the number of snapshots. For each snapshot, selecting a specific set of dyads has linear cost to the number of edges. Each window is only scanned once and therefore the time cost is linear in the number of windows. Moreover, since the number of windows is linear to the total number of snapshots even in the worst case (windows are overlapping, and window step is one), the algorithm is linear to the number of snapshots. Therefore, the time complexity is O( ¯MT), where ¯M is the averaged number of edges in each snapshot, and T is the number of snapshots.

The memory cost is low, and can be viewed as constant: for each snapshot, only the information of the tracked dyads is stored; information of dyads within the same window is aggregated; dyads information in the old window is overwritten once it is compared against the new window. And the space complexity is O(c), where c is a prescribed sample size. Theoretically the sample size should be proportional to the network size for good estimation. Our experiments show that a fixed sample size (to track 250 out of 1.2G dyads ≈ (50k 2 )) works well on a moderately large network.

5 Experiments and Results

We did thorough evaluation of our edge probability estimation based change point detection algorithm (called EdgeMonitoring for simplicity) on synthetic and real world datasets. For the synthetic datasets, the generative process is known, and we can compute the ground truth in the form of likelihood, which is naturally a baseline choice. We also use the state-of-the-art DeltaCon [Koutra et al., 2016] and LetoChange [Peel and Clauset, 2015] as two baselines.

5.1 Synthetic Data

Data generation 3 We generate a sequence of networks from a fixed generation model. The snapshots are not independent, each snapshot depends on the preceding one through the continuity parameter α(αt ≡ α). For each snapshot, each edge Figure 4: SBM, ground truth changes explained in Table 2. 1-α = 0.51 and window size = 20. DeltaCon (Fig b-top) and EM-KL (Fig b-bottom) have the smallest variance, but DeltaCon has two false negatives at 4©and 5©.

3 Generated using SNAP[Leskovec and Sosiˇc, 2016]

is selected independently with probability α, and if selected, the edge is again sampled from the generative model (Table 1). We introduce the change points by changing the generative model in the middle of the sequence of snapshots. Note that this change may be simply a change of parameter values for a given model (Eg. ER0.4 to ER0.6), or a change in the model type (Eg. SBM to ER), as well. Since our algorithm makes no assumptions about model specifics, we are able to detect both kinds of changes. We only inject parameter change in the synthetic experiment since the latter change is easily detectable. Sample changes are displayed in Table 2. The likelihood of the snapshot sequence is also provided.

We ran experiments with network sizes ranging from 1k to 50k, window size to be from 10 to 100 and continuity rate 1 - α to be 0.51 and 0.9. We generated a total of 5000 snapshots and sampled 250 edges uniformly at random to track. Both overlapping window (s=2η) and non-overlapping window have similar results, yet the latter is faster simply due to fewer windows. Hence we display non-overlapping window results only. For KL and KS, edges are grouped into 25 equalsized groups. We use upper 5% quantile as the threshold.

Results Figure 4 shows the qualitative comparison and Table 3 reports the efficiency. Figure 4a shows the likelihood of the network drops dramatically after the generative model changes, and recovers to new equilibrium afterwards. Our EdgeMonitoring (EM-Eu, EM-KL) approach can successfully identify all change points with 5X speed up over DeltaCon. The changes are explained in Table 2. We can see that EM-KL has the best performance: little fluctuation and perfect precision and recall. DeltaCon, though has smallest fluctuation, misses two change points. Both EM-KS and EM-Eu have large fluctuation. The quality of EM-KS heavily relies on the joint probability estimation, and we do see smaller fluctuation and higher recall for larger window size. EM-Eu in general has large fluctuation. EM-KL has the best overall performance, in terms of both quality and time efficiency. We believe grouping together with median selection contribute to its superiority.

5.2 Real World Data

Senate cosponsorship network([Fowler, 2006]) We construct a co-sponsorship network from bills (co-)sponsored in US Senate during the 93rd-108th Congress. An edge is formed between two congresspersons if they cosponsored the same bill. Each bill corresponds to a snapshot, and forms a clique of co-sponsors. A window is set to include all bills in a single Congress (Biennially).

We randomly selected 250 dyads and tracked their fluctuations across the Congresses. We start from the 97th Congress since full amendments data is available only from 97th session onwards. Figure 5 compares EdgeMonitoring+KL, DeltaCon and LetoChange. All methods were able to detect the most significant change point at the 104th Congress. Fowler [Fowler, 2006] points out that there was a "Repub-

• 1 EM for EdgeMonitoring (running time includes both KL and Euclidean), DC for DeltaCon, LC for LetoChange. EM and DC are implemented in MATLAB while LC in Python. All run on a commercial desktop with 60hrs as time limit. Each running time averaged over 5 runs.
• 2 BTER dataset has 800 snapshots

lican Revolution" in the 104th Congress which "caused a dramatic change in the partisan and seniority compositions." The author also points out the significance of the 100th (highest clustering coefficient, significant collaboration) and 104th Congress (lowest clustering coefficient, low point in collaboration) as inflection points in the Senate political process. Both our EdgeMonitoring approach and LetoChange classify these two Congresses as change points, but the latter takes much more time. DeltaCon picks up on one (104th) and not the other (100th). This provides evidence that our algorithm is able to capture the changes in network evolution effectively while being significantly faster than the state-of-the-art.

6 Conclusion

We develop a change point detection algorithm for dynamic networks that is efficient and accurate. Our approach relies on sampling and comparing the estimated joint edge (dyad) distribution. We first develop a maximum likelihood estimator, and analyze its drawbacks for small window sizes (the typical case). We then develop a consistent and unbiased estimator that overcomes the drawbacks of the MLE, resulting in significant quality improvement over the MLE. We conduct a thorough evaluation of our change point detection algorithm against two state-of-the-art approaches on both synthetic and real-world data. Our results indicate that our method is up to 9X faster than the state-of-art while achieving better quality. In the future we plan to extend our work to track higher order structures of the network such as 3-profiles [Elenberg et al., 2015] or 4-profiles and see how they evolve over time.

Acknowledgments

This work is supported in part by National Science Foundation (NSF) grants DMS-1418265, IIS-1550302 and IIS1629548. Opinions and findings herein are of the author(s) and not necessarily the views of the NSF.