Latent Evolution Model for Change Point Detection in Time-varying Networks

In this section, we briefly introduce some related literatures of CPD, which can be used to discover outliers in network and sequences [18, 19, 20, 21]. [22] provided the first attempt on transforming the CPD problem to detect the times of fundamental changes. The proposed LetoChange method used Bayesian hypothesis testing to determine whether a model with parameter changes is acceptable. [3] proposed a Markov process based method, which treats the time-varying network as a first order Markov process. This work used joint edge probabilities as the ”feature vector”. [10] proposed a tensor-based CPD model. Tensor factorization is widely recognized as an effective method to capture a low-rank representation in the time dimension. In this method, after getting the latent attributes of each time point, the clustering and outlier detection methods were used to identify change points. Activity vector method was proposed by [9], it focused on finding anonymous networks which are significant deviation from previous normal ones. Activity vector leveraged the principal eigenvector of the adjacency matrix to present each network, it considered a short term window to capture the normal graph pattern. Currently, [7] proposed a Laplacian-based detection method for CPD, named LAD. It utilised the Laplacian matrix and SVD to obtain embedding vector of each network. In additional, two window sizes are involved to capture both short and long term properties. [15] proposed an efficient approach for network change point localization. However, this method only considered the changes in link connections that cannot capture the dynamics for weight graphs.

As aforementioned discussion, most of the existing relevant methods did not fully consider the dynamic trends in the time-varying networks when detecting change points.

Recent studies have indicated that the latent space learning methods provide excellent solutions for the dynamic network problems, e.g., flow prediction [17], community detection systems [23], social network evolution learning [24], sequential analysis [25, 26, 27, 28], and traffic speed prediction [16]. The benefits of a latent space model is to obtain the low-rank approximations from the original features, so that obtain a more compact matrix by removing the redundant information of the large scale matrix [29, 30]. [31] proposed a dynamic model for a static co-occurrence model CODE, where the coordinates can change with time to get new observations. There are several advanced investigations focusing on traffic analysis under the latent space assumption. For example, to address the issues of finding abnormal activities in crowded scenes, [32] devised a spatio-temporal Laplacian eigenmap strategy to detect anomalies. [33] leveraged a tensor structure to resolve the high-order time-varying prediction problem. The key characteristics of the original time sequences can be represented by the learned latent core tensor. [16] developed a real-time traffic speed prediction model that learns temporal and topological features from the road network. It then used an online latent space learning method to solve the proposed constrained optimization problem. [34] took the advantage of dynamic latent space learning to address the temporal analysis problem very efficiently on the large-scale streaming data. In other applications, [4, 17] used the online matrix factorization model to solve the crowd flow distribution problem.

In summary, latent space models can effectively capture dynamic evolving trends from real-world networks, as well as learning the latent embedding for such networks.

Considering inherent complexities of the time-varying network, we leverage the latent space model to represent dynamic graphs. Latent space models are widely applied into solving network-wide problems, e.g., community detection [23], social and heterogeneous networks construction[1], and network-wide traffic prediction[17, 35], etc. The benefit of low-rank approximation of latent space model is to get a more compact matrix by removing the redundant information of the large-scale matrix [36]. For a slice $\bm{G}$ $\in$ $\mathbb{R}^{n\times n}_{+}$ (the timestamp $t$ is omitted for brevity), we use tri-factorization that decomposes $\bm{G}$ into three latent representations, in the form $\bm{U}\bm{C}\bm{V}$, where $\bm{U}$ $\in$ $\mathbb{R}^{n\times k}_{+}$ and $\bm{V}$ $\in$ $\mathbb{R}^{k\times n}_{+}$ denote the latent attributes of start and end nodes respectively, $\bm{C}$ $\in$ $\mathbb{R}^{k\times k}_{+}$ denotes the attribute interaction patterns, $n$ is the number of nodes and $k$ is the number of dimension of latent spaces. Therefore, we approximate $\bm{G}$ by a non-negative matrix tri-factorization form $\bm{U}\bm{C}\bm{V}$:

Note that, some tri-factorization methods prefer to use the factorization form $\bm{G}\approx\bm{U}\bm{C}\bm{U}^{\top}$ because they usually focus on the undirected graph [23, 37]. However, in our problem, some real-world networks are formulated as the directed graph. In this way, $\bm{G}\approx\bm{U}\bm{C}\bm{V}$ provides a more flexible formulation to represent the graph $\bm{G}$ [38, 39]. Thus, our basic problem formulation differs from that of methods in [16, 4, 17].

The reasons that we utilise the non-negativity constraint are: 1) the reasonable assumptions of latent attributes and the better interpretability of the results [40, 41]; 2) all weights in our networks are non-negative, thus latent spaces $\bm{U}$, $\bm{V}$ and $\bm{C}$ should be non-negative as well.

To take the time-series pattern into consideration, we formulate this sequential problem as a time-dependent model. Given a certain time window $T$ (i.e., a set contains $T$ previous networks, and $T$ presents the current time as well), for each timestamp $t$ $\in$ [1, $T$], our method aims to learn the time-dependent latent representations $\bm{U}_{t}$, $\bm{V}_{t}$ and $\bm{C}$ from $\bm{G}_{t}$. Although the latent attribute matrix $\bm{U}_{t}$ is time-dependent, we assume that the latent attribute interaction matrix $\bm{C}$ is an inherent property, and thus $\bm{C}$ is fixed for all timestamps. This assumption is based on the fact that temporal network exists invariant properties [16]. After considering the temporal pattern, our model can be formulated as:

Since our model is a time-varying forecasting system, the network is continually changing over time. We not only want to learn the time-dependent latent attributes, but also aim to learn the evolution patterns of latent attributes $\bm{U}_{t}$ and $\bm{V}_{t}$ from previous timestamps to the next, i.e., how to learn the evolution forms such as $\bm{U}_{t}$ $\rightarrow$ $\bm{U}_{t+1}$ and $\bm{V}_{t}$ $\rightarrow$ $\bm{V}_{t+1}$.

To this end, we define two transition matrices $\bm{A}\in\mathbb{R}_{+}^{k\times k}$ and $\bm{B}\in\mathbb{R}_{+}^{n\times n}$ to represent the smooth trends of $\bm{U}_{t}$ and $\bm{V}_{t}$ in T previous conditions. The learned transition patterns $\bm{A}$ and $\bm{B}$ can be treated as a global evolution trend from previous timestamps. Thus, matrices $\bm{A}$ and $\bm{B}$ can be recognized as the evolving patterns of the dynamic network. For example, $\bm{A}$ approximates the changes of $\bm{U}$ between time t-1 to t, i.e., optimizes $\bm{U}_{t}=\bm{U}_{t-1}\bm{A}$.

Then, let us involve all timestamps in a window-size $T$, our latent transition learning can be expressed as:

When we consider the temporal topology problem jointly, we have:

where $\lambda_{1}$ is the regularization parameter.

To this stage, transition matrices $\bm{A}$ and $\bm{B}$ can learn the short-term knowledge from previous networks in a window size. Considering that the long-term knowledge is helpful for the CPD and network evolution problems [3, 17], we also take the long-term pattern into consideration to improve the prediction performance in the next section.

Dynamic networks, such as climate and travel networks, usually have a strong periodic property. Such property is considered as the long-term guidance which has more gradual trends. Because of this point, we can leverage a long-term window to achieve the large-scaled network evolution, and propose a novel long-term guidance to extract periodic and stable information. Given a long-term set of networks $\{\bm{G}_{t}\}_{t=1}^{4T}$ (the window size extends from $T$ to 4$T$ for example, then the long-term guidance will consider the states of 4$T$ previous networks), we learn two guide matrices, $\bm{U}^{\rm lt}$ and $\bm{V}^{\rm lt}$, to represent the common characteristic of $\{\bm{G}_{t}\}_{t=1}^{4T}$:

where $r_{t}$ is the adaptive weight calculated by function: $r_{t}=\frac{e^{s_{t}}}{\sum^{4T}_{t=1}e^{s_{t}}}$, and $s_{t}$ is the similarity between $\bm{G}_{t}$ and the current network $\bm{G}_{4T}$. We use the matrix L₂-norm to calculate the similarities among networks, and then normalize them satisfying $\sum^{4T}s_{t}=1$. The assumption behind this is that $\bm{G}_{t}$ can provide more useful information if it is more similar to $\bm{G}_{4T}$. $\odot$ is the Hadamard product operator.

The main idea of the above process is analogous to the multi-view learning which projects multiple views into shared common spaces [42, 21]. The adaptive weight $r_{t}$ controls how much information can be extracted from $\bm{G}_{t}$. Such technique has been successfully used in common space learning [43, 44]. The learning approaches of Eq. (5) are:

After having $\bm{U}^{\rm lt}$ and $\bm{V}^{\rm lt}$, the strategy of the long-term guidance is formulated as:

where $\bm{U}_{T}$ and $\bm{V}_{T}$ denote the latent representations at current timestamp $T$.

Taking all techniques jointly into consideration, the final optimization function is expressed as:

where $\lambda_{2}$ is the regularization factor.

Make a Prediction. The learned latent space matrices $\bm{U}_{T}$ and $\bm{V}_{T}$, and the evolving patterns $\bm{A}$, $\bm{B}$ and $\bm{C}$ can be used to make the prediction after solving Equation (9). Then we have the predicted network $\hat{\bm{G}}_{T+1}$ as:

Due to the non-convexity of Eq. (9), we use an effective gradient descent approach with the multiplicative update method [45] to find its local optimization.

By computing the prediction network $\hat{\bm{G}}_{T+1}$ for timestamp $T$+1, we can compare it with the actual network $\bm{G}_{T+1}$ to check whether $\bm{G}_{T+1}$ is a change point or not. Specifically, we use the Laplacian-spectrum-based method which has shown a state-of-the-art performance in graph comparisons [46]. The eigenvalues of the Laplacian matrix capture the structural properties of the corresponding graph. Note that real-world applications contain both undirected and directed graphs. For the undirected graphs, we construct a normalized Laplacian matrix $\bm{L}_{t}$, defined as $\bm{L}_{t}=\bm{I}-\bm{D}_{t}^{-1/2}\bm{W}_{t}\bm{D}_{t}^{-1/2}$, where $\bm{I}$ is the identity matrix; $\bm{W}_{t}$ is the weight matrix of $\bm{G}_{t}$ and $\bm{D}_{t}$ is a diagonal matrix $\bm{D}_{t,(i;i)}$ = $\sum_{j}(\bm{G}_{t,(i;i)})$. For the directed graphs, we leverage the method proposed in [47] to build the Laplacian matrix:

where $\bm{\phi}$ is the matrix with the Perron vector [48] of $\bm{P}$ in the diagonal; $\bm{P}$ is a transition probability matrix [47].

Given the Laplacian matrices for all timestamps, i.e., $\bm{L}_{1}$, $\cdots$, $\bm{L}_{T+1}$ and $\hat{\bm{L}}_{T+1}$, where $\bm{L}_{T+1}$ is built by actual network $\bm{G}_{T+1}$; $\hat{\bm{L}}_{T+1}$ is built by the prediction $\hat{\bm{G}}_{T+1}$. We follow a well-performed metric that utilises the singular values calculated of Laplacian matrices $\bm{L}_{t}$ as the network attributes [7]. Then, the vector made up of singular values is denoted as $\overrightarrow{\bm{\sigma}}_{t}$, $t$ $\in$ [1, $T$+1]. Let $\overrightarrow{\bm{\sigma}}_{a}$ ($\overrightarrow{\bm{\sigma}}_{T+1}$) denote the vector of actual network, $\overrightarrow{\bm{\sigma}}_{p}$ (singular values of $\hat{\bm{L}}_{T+1}$) denote the vector of prediction network. The first anomaly score $\bm{Z}_{1}$ based on the cosine distance is defined as:

It is evident that the more dissimilar between the prediction and the actual network, $\bm{Z}_{1}$ is closer to 1, and vice versa¹¹1The range of $\bm{Z}_{1}$ is [0,1] because of the non-negativity of $\overrightarrow{\bm{\sigma}}$.

Here we analyse the time complexity of LEM-CPD algorithm. After the prediction process, the time complexity of the change point detection is mainly affected by SVD. Even though randomized SVD has an $O(n^{2}k)$ complexity, it is not involved in the update loop of variables (such as $\bm{U}_{t}$ and $\bm{V}_{t}$). In the prediction process, the time complexity of Eq. (11) - Eq. (12) is dominated by the matrix multiplication operations in each iteration. Therefore, for each iteration, the computational complexity of LEM-CPD is $O(Tkn^{2})$. Compared with the latest CPD method LAD ($O(n^{2}k)$) [7], the time complexity of our method is $O(Tkn^{2})$, which is higher than LAD.

Setting: We chose a widely used data generator Stochastic Block Model (SBM) [49] to generate synthetic data. We followed two settings in the data generation process in [7]. 1) Pure datasets: four time-varying networks with 500 nodes and 151 time points are generated, which contain 3, 7, 10, 15 change points respectively. Therefore, we can use HR@K, where K is setting to 3, 7, 10, 15 to evaluate all competitive methods. Hybrid dataset: we generate four homogeneous networks of the Pure datasets, while contains both change points and event points according to the hybrid setting in [7]. Thus, HR@K (K = 3, 7, 10, 15) is also suitable in this dataset. In the experiments, our short and long term window sizes are set to 3 and 12 respectively; the change point detection threshold is set to 0.5; the trade-off factor $\alpha$ is set to 0.2.

Discussion: Table 1 presents the HR@K accuracy of all test methods on eight datasets. It is apparent that there are several methods which can effectively handle the CPD problem in the synthetic datasets, e.g., LEM-CPD (ours), LAD, and EdgeMonitoring. It is mainly because the normal evaluation pattern in the synthetic datasets is very stable, so that these method can easily detect anomalies if some dramatically changes appeared. TensorSplat method cannot capture all change points in the sparse datasets. Unlike LADdoc, LAD and LEM-CPD, Activity vector chooses the principal eigenvectors of the adjacency matrix to represent the graph embedding directly, the latent space may not be learned well compared with the Laplacian spectrum method. Furthermore, as the results shown in Table 1, we can see that without the long-term guidance, LEM-CPD_-lt performs worse than our final model, which demonstrates the effectiveness of the long-term guidance. For methods that use the prediction process, such as LEM-CPD_-lt and LEM-CPD, can achieve the best detection accuracy as shown in Table 1.

Figure 3 provides the detection results and their $\bm{Z}$ scores in the Pure dataset with seven change points. The seven most anomalous points are correctly detected. We can clearly see that the $\bm{Z}$ scores of change points are markedly different from the normal ones. The scenario proves our point that the normal evolutionary pattern of the synthetic network is stable. Therefore, in our trade-off strategy, the score $\bm{Z}_{2}$ can greatly handle this situation, so that we use a small $\alpha$ in the experiments.

Setting: To evaluate the performance of our method in the real-world time-varying networks, we conduct experiments on three datasets. The selection of K of HR@K is based on the number of change points in the corresponding dataset. For example, there are two change points in the Senate dataset, then K will be chosen from 1 and 2.

1) Flow dataset: a large-scale, real-world dataset that is collected from the city transport. This data records the crowd flows among the entire city train network. In a weekday, there exist six change points during the rush and non-rush transition period. 2) Maintenance dataset: The Maintenance of metro lines will change the connections among stations. In this dataset, the one line of the city train network had been maintained on three separated days. Therefore, three change points are included in the dataset (21 days in total). 3) Senate dataset: a social connection network between legislators during the 97rd-108th Congress [50]. In this dataset, 100th and 104th congress networks are recognized as the change points in many references [3, 7]. The performances of different methods are summarized in Table 2.

Discussion: Our model, LEM-CPD significantly outperforms all other comparative methods as shown in Table 2, especially LEM-CPD has outperformed in the Flow dataset. In this test, we report the average accuracy of HR@3 and HR@6 in five weekdays. Compared with LAD and Activity vector, the prediction strategy can capture much more dynamic trends with time evolving. As the motivation illustrated in the Introduction, LADdoc, LAD and Activity vector may difficult to identify the natural gaps between the previous normal graph pattern and the actual network. TensorSplat and EdgeMonitoring achieved the second and third positions because they considered a part of evolution trend during their learning process, e.g., tensor-based methods are able to get the evolution information. In the real-world experiments, we set up a larger $\alpha$=0.6 to let the prediction score participate more in the final decision. There presents similar conclusions in the Maintenance dataset because the city train network has the strong dynamic characters. However, in the Senate dataset, LADdoc, LAD and EdgeMonitoring methods can also get 100% accuracy, it is because the 100th and 104th congress networks are significant deviation from others.

This section evaluates the performances of LEM-CPD by varying the critical parameters ($\alpha$, $\lambda_{1}$ and $\lambda_{2}$). We here show the experimental results for the Flow dataset. We discuss them separately but pick them up by the grid search method because four parameters have high dimensional correlations that are hard to visualize. Our illustration approach that discusses parameters separately has been widely used in many other research papers [16, 4].

Fig. 5 (a) shows the different performances with a varying setting for $\alpha$. The $\alpha$ indicates the balance factor to control the importance of prediction results. When we increase $\alpha$ from 0.1 to 0.6, the results improve significantly. However, the performance tends to stay stable at 0.5 $\leq$ $\alpha$ $\leq$ 0.9. In particular, LEM-CPD achieves the best result when $\alpha=0.6$, while it can get good performance if the $\alpha$ is set between 0.5 and 0.9. It indicates that the prediction result is important and useful to the change point detection.

Fig. 5 (b) and (c) reveal the effect of varying $\lambda_{1}$ and $\lambda_{2}$. These two parameters determine the strength of the transition matrices and the long-term guidance, respectively. $\lambda_{1}$=$2^{-1}$ and $\lambda_{2}$=$2^{3}$ yield the best prediction results for LEM-CPD.