{'fraction_non_alphanumeric': 0.04577832432307014, 'fraction_numerical': 0.008190445253950207, 'mean_word_length': 4.832205683355887, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
H-VGRAE: A Hierarchical Stochastic Spatial-Temporal Embedding Method for Robust Anomaly Detection in Dynamic Networks

Detecting anomalous edges and nodes in dynamic networks is critical in various areas, such as social media, computer networks, and so on. Recent approaches leverage network embedding technique to learn how to generate node representations for normal training samples and detect anomalies deviated from normal patterns. However, most existing network embedding approaches learn deterministic node representations, which are sensitive to fluctuations of the topology and attributes due to the high flexibility and stochasticity of dynamic networks. In this paper, a stochastic neural network, named by Hierarchical Variational Graph Recurrent Autoencoder (H-VGRAE), is proposed to detect anomalies in dynamic networks by the learned robust node representations in the form of random variables. H-VGRAE is a semi-supervised model to capture normal patterns in training set by maximizing the likelihood of the adjacency matrix and node attributes via variational inference. Comparing with existing methods, H-VGRAE has three main advantages: 1) H-VGRAE learns robust node representations through stochasticity modeling and the extraction of multi-scale spatial-temporal features; 2) H-VGRAE can be extended to deep structure with the increase of the dynamic network scale; 3) the anomalous edge and node can be located and interpreted from the probabilistic perspective. Extensive experiments on four real-world datasets demonstrate the outperformance of H-VGRAE on anomaly detection in dynamic networks compared with state-of-the-art competitors.


I. INTRODUCTION
Anomaly detection for networks (graphs) has been an active topic in various research areas, such as social media, computer security, and e-commerce [1]- [3]. Numerous studies have been conducted on static networks in the past two decades by capturing deviations of structural patterns and show great success in practical applications [2]. However, with the rapid growth of data and users, the update of network topology and node attributes becomes increasingly frequent. Typical examples include IP-IP network traffic logs, item buying and comments in e-commerce sites, system components logs in cyber-physical systems, and so on. The patterns of anomalies are thus time-varying and difficult to capture if merely applying methods from static networks [4]. Due to the complex spatial-temporal (ST) patterns, identifying anomalous edges and nodes in such dynamic networks remains a challenging problem.
A general framework to detect anomalies in dynamic networks is to extract the features that contribute to the anomaly detection, and then compute anomaly scores based on the features. The critical problem is how to extract anomaly-related features. Some studies predefine the anomalies by considering structural changes, such as sudden changes in connectivity, node-degree distribution, and so on [5], [6]. However, the anomalous patterns are much more than sudden structural changes and could be very complicated in practical scenarios. It is thus necessary to implicitly explore the techniques that can extract anomaly-related features implicitly in the datadriven mode.
Recently, a novel feature selection technique called network embedding has shown promising results on extracting features of networks, especially on preserving structural information in the low-dimensional representations [7], [8]. The learned features can be utilized in downstream tasks, such as node classification, link prediction, as well as anomaly detection, without having to define them explicitly. Due to the anomaly diversities and the lack of labels for training, the anomaly detection methods based on network embedding has to work in the semi-supervised mode [9]. The model needs to be trained to capture normal patterns from the training set and detects anomalies according to the deviations from the learned normal patterns.
In this paper, we also focus on detecting anomalies in dynamic networks based on network embedding, similar to NetWalk [10] and AddGraph [11]. However, since dynamic networks often interact with the environment and human operators to show complex evolving patterns, their behaviors can exhibit both stochasticity and strong spatial-temporal (ST) correlations. The network embedding technique used in Net-Walk [10] and AddGraph [11] can only generate deterministic latent representation, which cannot model the stochasticity well and could lead to lack of robustness in anomaly detection when facing rapid dynamics. Moreover, limited with the utilized neural network models in NetWalk [10] and AddGraph [11], the ST information is still relatively coarse and could also cause the models sensitive to tiny changes. Therefore, the goal of this paper is to learn robust node representations to capture normal patterns of dynamic networks by modeling stochasticity and complex ST dependencies. Ideally, with only a few anomalous edges, the node representations would be dominant by the normal pattern, and thus the anomalous edges and nodes would have low reconstruction probabilities. There are two major challenges to achieve our goal.
The first challenge is how to combine stochasticity and ST dependency together in the latent representation and enhance the representation capability. The combination of stochasticity and temporal dependence has been conducted in [12] by unifying Recurrent Neural Network (RNN) and variational autoencoder (VAE) [13], and be extended to graph domain in Variational Graph Recurrent Neural Network (VGRNN) [14] by additionally extracting spatial features using graph neural networks (GNN). Unfortunately, constrained by GNN and RNN, which only consider the one-hop neighbors in node representation and temporal correlations in adjacent timestamps, the ST features extracted by VGRNN may lead the detector sensitive to local and short-term changes. And the capability of the latent representations is limited due to the single layer of stochastic variables in VGRNN .
To strengthen the latent representations, we propose a Hierarchical Variational Graph Recurrent autoencoder (H-VGRAE), which is scalable to deep structure and considers multi-scale ST information. H-VGRAE encodes stochasticity, multi-scale ST and content features into the hierarchical stochastic latent representation. H-VGRAE has three major advantages in architecture: a) Through hierarchical ST blocks connected with multiple stochastic layers, H-VGRAE can capture multi-scale ST patterns, which has been shown to be useful to strengthen robustness in anomaly detection tasks [15]. b) By sharing information about the conditional stochastic latent variables between the encoder and the decoder, H-VGRAE can model more complex stochasticity by stacking multiple layers of stochastic variables, alleviating the difficulty to train deep conditional stochastic layers. c) H-VGRAE can model non-Gaussian distributions by employing Normalizing Flow [16] in multiple correlated stochastic layers. Benefited from these advantages, H-VGRAE is capable of modeling complex patterns of normal behaviors in the dynamic networks and provides rich information to be adjust to high variabilities.
Given the learned node representations, the second challenge is how to detect anomalies accurately and interpretably. Considering the semi-supervised mode, the H-VGRAE is essentially maximizing the evidence lower bound (ELBO) of the likelihood of the adjacency matrix and node attributes in the normal training set. In the online detection stage, the H-VGRAE captures the normal parts of an anomalous snapshot in the latent representations and reconstructs the normal edges and node attributes with high probability. The multi-scale features and stochasticity modeling enhance the robustness of the latent representations. The remaining ones with low reconstruction probabilities are detected as anomalies, which allows estimating the contribution of each dimension of node attributes.
Moreover, most of deterministic embedding based methods rely on single anomaly detection mechanism, such as clustering in latent representation [10] and specifically designed anomaly score [11]. In H-VGRAE, the anomaly detection combines reconstruction and prediction by designing the prior to predict the latent variables based on historical states. Benefited from prediction prior, H-VGRAE essentially detects anomalies in two hybrid perspectives: 1) the input that cannot be reconstructed well is an anomaly; 2) the input that deviates from the prediction is an anomaly.
To summarize, the main contributions of this paper are as follows: • We analyze the likelihood of adjacency matrices and node attributes in multiple snapshots and obtain the evidence lower bound of the likelihood. • We propose a novel semi-supervised network embedding method, H-VGRAE, which can capture normal patterns in dynamic networks by maximizing the evidence lower bound. H-VGRAE can jointly utilize stochasticity, spatial-temporal information, and content information to empower the embeddings. • We design multiple mechanisms to strengthen the robustness of the node representations learned by H-VGRAE, including the extraction of multi-scale ST features, the generalization of non-Gaussian random variables, the predictive prior, and the information-sharing mechanism between the inference of posterior and prior to train deep hierarchies of conditional stochastic layers. • The anomalies are detected based on the conditional reconstruction probabilities. The level of anomalies can be compared in different dynamic networks. • Through extensive experiments on four real-world dynamic network datasets, we show the effectiveness and superiority of H-VGRAE in terms of AUC score. In dynamic networks with different scales, H-VGRAE exhibits great robustness, with AUC higher than 0.75 even with 10% anomalous samples in the most complex and attributed Github network.

II. RELATED WORK
In this section, we briefly review network embedding techniques and existing anomaly detection methods for dynamic networks.

A. Network Embedding
Network embedding assigns nodes in a network to lowdimensional representations, which effectively preserves the network structure and support network inference [8]. For dynamic networks, methods have been developed using various techniques such as matrix factorization [17], [18], random walk [10], [19], neural networks [14], [20]- [22], and stochastic process [23], [24]. Most of these methods focus on sudden changes in either structural features or content features (from node attributes) to learn dynamic node representations [25], and only a few of them model both kinds of changes simultaneously [14], [24]. Moreover, due to the simplification of temporal patterns, such as smoothness [18] and single granularity [14], the capability of existing models is limited.
Recent developments [12] in speech sequence have shown that the stochasticity can be more precisely captured by stochastic variables with properly estimated probability distributions than deterministic variables. Inspired by the stochastic models in sequence processing, a recent study VGRNN [14] introduces stochasticity into dynamic network embedding to model the uncertainty of the latent representations and obtain state-of-the-art performance in link prediction task. However, constrained by the stochastic RNN in their models, the layers of stochastic variables in VGRNN cannot be deepened and thus cannot capture multi-scale ST features well.

B. Anomaly Detection in Dynamic Networks
According to the utilization of network embedding, the anomaly detection methods in dynamic networks can be classified into two categories (for other classification criterion, see [2] for survey): • Anomaly detection without using network embedding: Most of existing methods define anomalies as sudden changes of structural patterns. The structural patterns are captured through connectivity [4], sketch [26], dense subgraphs [5], [27], [28], and so on. Some works design approximation mechanisms of the structural patterns to achieve near real-time pattern updating [5], [6], [29]. Then, the anomaly scores are defined according to the structural patterns and the anomalies are detected by setting a anomaly score threshold. • Anomaly detection using network embedding: Recently, some works begun to leverage network embedding techniques to implicitly define anomalous patterns in dynamic networks. In [30], a commute time distance is used to localize anomalous changes in dynamic graphs, while the method mainly focuses on sudden structural changes and cannot capture anomalies with complex temporal patterns. NetWalk [10] utilizes an autoencoder to generate low-dimensional latent representations of nodes based on random walks in each timestamp, and detect anomalies by clustering these node representations dynamically. AddGraph [11] extends NetWalk by using a spatial-temporal neural network to capture more complex patterns in dynamic networks, and employing negative sampling to build an end-to-end semi-supervised learning model.

A. Problem Formulation
Dynamic network model. Let T max be the maximum timestamp. A network stream with length T max can be represented as G = {G(t)} Tmax t=1 , where G(t) = (V(t), E(t)) is the entire snapshot at timestamp t, and V(t) and E(t) are the sets of nodes and edges respectively. the maximum number of nodes in all the T max timestamps. The topology of the network G(t) can be represented by an adjacency matrix means there exists an edge from node i to node j at timestamp t. In practical scenarios, the dynamic may also have attributes attached on each node, which can be denoted as X(t) = (x 1 (t), . . . , x N (t)) ∈ R N ×D , where D denotes the dimension of node attributes that is constant across time.
Problem statement. For dynamic network anomaly detection, the objective is to detect anomalous edges and nodes at any given timestamp t, i.e., in real time as G(t) occurs. As in time series modeling, historical observations can be beneficial for understanding current data. Therefore, the streams are used to detect anomalies, instead of just using A(t) and X(t).

B. Overall Structure
As shown in Fig.1, the overall structure of using H-VGRAE to detect anomalies is composed of two stages: offline model training and online anomaly detection. In the offline stage, the proposed H-VGRAE is trained on normal training data A(≤T train ) and X(≤T train ) (T train <T max ) to capture normal behaviors of the dynamic network by maximizing the evidence lower bound of the likelihood. The reconstruction probabilities are used as anomaly scores. Then the thresholds of anomaly scores to detect anomalous edges and nodes are set based on the anomaly detection performance on a validation dataset. In the online stage, the trained H-VGRAE is used to calculate the anomaly scores of a new network observation G(t )(t <T max ) at timestamp t , based on the all the current and previous online observations. If the anomaly score of some edge/node in G(t ) is higher than the threshold, it will be detected as an anomalous edge/node, otherwise, it is normal. For the detected anomalous edges/nodes, the interpretations are made by the reconstruction probabilities of these edges/node attributes.

C. Basics of Graph Autoencoder and Stochastic RNN
In network embedding, Graph Autoencoders (GAEs) [31] are commonly used for their ability of modeling structure information. Inheriting the encoder-decoder framework from the autoencoders, GAEs can learn low-dimensional latent representations of node features in the unsupervised mode. The encoder is implemented with Graph Neural Network (GNN), which is able to fuse the content and structural information of a node with its neighbors. Typical GNNs include GCN [32], GAT [33] and GraphSAGE [34]. Then, the structure information is reconstructed by a link prediction decoder, like bilinear decoder.
Stochastic RNN is a variant of Bayesian model [12], which contains a Variational Autoencoder (VAE) [13] at every timestamp. The latent random variable z(t) in Stochastic RNN are conditioned on the state variable h(t − 1). This condition helps the Stochastic RNN modeling highly variable temporal patterns of time series. With input x(≤t), the Stochastic RNN have the inference stage q φ (z(t)|x(≤t, z(<t)) and the generation stage p θ (x(t)|x(<t), z(≤t)) with a prior p prior (z(t)|z(<t)). For real-value inputs, all the three probability distributions can be set to diagonal Gaussian N (µ(t), σ 2 (t)I). To model temporal dependency, there is a recurrence stage that updates the hidden state where f (·) is a deterministic non-linear transition function such as LSTM [35] or GRU [36]; ϕ x (·) and ϕ z (·) are neural networks to extract features. The parameters of the Stochastic RNN can be learned in a end-to-end mode, by maximizing the ELBO of likelihood: Due to the simple assumption of diagonal Gaussian, the KL divergence may vanish during training, leading to the generation stage trapped in sub-optimal areas. Furthermore, as shown in [37], the Stochastic RNN inherits the limitation of VAEs, which are difficult to optimize for deep hierarchies due to multiple layers of conditional stochastic layers. The KL vanishing problem and the shallow structure can degrade the capacity of the Stochastic RNN, which is crucial to model highly variable temporal patterns.

IV. PROPOSED METHOD
In this section, we first present the network architecture of our proposed H-VGRAE. Then, we describe the offline and online stages in H-VGRAE, respectively. H-VGRAE is in the general framework of stochastic RNN, which is composed of an encoder to encodes the input into low-dimensional latent representations and a decoder to reconstruct the input based on the latent representations, as shown in Fig.2. The basic idea of H-VGRAE is the following. First, it combines the GAE and RNN to capture ST patterns and content information from the dynamic of network topology and node attributes. Second, the current ST features and historical state are encoded into low-dimensional random variables, which can model high variability of the data patterns. Third, inspired by the utilization of multiscale feature extractors in the speech processing literature [15], we design a hierarchical structure to extract multiscale spatial-temporal features and encode them to multilayer stochastic latent variables. We then design a information sharing mechanism between random variables in the encoder and the decoder to alleviate the difficulties of optimizing multiple stochastic layers. Forth, we employ the Normalizing Flow [16] on the random variables to model non-Gaussian distributions.

A. H-VGRAE
The model of H-VGRAE is based on the likelihood maximizing like Stochastic RNN, but have an additional input of X(t). For the T -length snapshots of A(≤T ) and X(≤T ), the likelihood with latent variable Z(≤T ) is where, It is to be noticed that the ELBO in Eq. 3 is different from the loss function of VGRNN [14]. The loss function of VGRNN does not have the term of p(X(t)|A(<t), X(<t), Z(≤t)), which means it cannot guarantee to maximize the essential ELBO.
The task is to maximize the ELBO by make estimations of the posterior, the prior and the reconstruction probabilities. In H-VGRAE, the inference of the three probability distributions is implemented by an Encoder-Decoder framework. A possible implementation of the inference network is composed of GAE and Stochastic RNN to extract spatial and temporal features respectively. However, limited by the local approximation of GNN and the memory capability of RNN, the inference of the three distributions can only focus on local spatial features and short-term temporal features. The limited reception field could cause inaccurate inference and make the anomaly detection sensitive to local and short-term changes.
To tackle this problem, we design a hierarchical spatialtemporal inference model to extract multi-scale ST features and encodes them to hierarchical stochastic latent variables. The network architecture of the encoder and the decoder in  H-VGRAE is shown in Fig. 3 and 4 respectively. In the encoder, it stacks multiple GNN layers and dilated RNN (DRNN) layers [38] to extract multi-scale ST and content features from A(≤t) and X(≤t), and then encodes current feature information and historical state into stochastic latent representation , d is the dimension of the latent variables. In the decoder, it uses Z(t) to reconstruct the node attributes X(≤t) and the adjacency matrices A(≤t). Similar with the encoder, the decoder first decodes the multi-scale content and ST features of Z(t), then reconstruct A(≤t) and X(≤t)) by a Bernoulli MLP and a Gaussian MLP respectively.

1) Encoder:
The structure of the encoder is shown in Fig.3. The encoder is composed of a content feature extraction block, a multiscale ST feature extraction block and a hierarchical stochastic block. The content feature extraction block ϕ x (·) is be implemented with MLP. In the M -scale ST feature extraction block, there are M stacked GNN layers, and each is followed by a DRNN layer. Following each DRNN layer, the latent variable where the activation functions of GNN, GNN Linear and GNN Softplus are ReLU, none and Softplus.
The procedure of the encoder at time t is as follows.
Step 1: The node attributes X(t) are inputted to ϕ X (·) to extract content features C X (t) = ϕ X (X(t)). ϕ X (·) are deep neural networks which operate on each node independently, and are crucial for learning complex content features.
Step 2: C X (t) are inputted into a sequence of GNNs (M in total), and i-th scale of spatial features is The GNN can be implemented with GCN [32], whose calculation process in the i-th layer can be formulated as follows: where σ(·) is a activation function with non-linearity. In our model, we choose the ReLU(x) = max(0, x) to be the activation function.
is the regularized adjacency matrix with self loops, whereÃ(t) = A(t)+I N ×N is the adjacency matrix with self loops andD is the diagonal degree matrix, whereD ii (t) = jÃ ij (t) denotes the degree of node i andD ik (t) = 0 for k = i. A single layer of GNN can represent each node by fusing the data of this node and its 1-hop neighbors. As proved in [32], stacking i layers of GNN is equivalent to represent each node by considering its 1-hop to i-hop neighbors' data, which is corresponding with the i-th spatial scale in our model.
Step 3: Each S The correspondence of spatial scales and temporal scales is based on the nature of the dynamic networks, that the information propagation process is controlled by the distance between nodes. Supposing that the messages transmitted by some node include both its own messages and its neighbors' messages, and the time consumed during transmission is only proportional to the hops between two nodes. The more distant nodes which are considered in upper GNNs should have less transmission frequencies, which means those messages should correspond to the long-term trend.
To model these different frequencies, DRNN (i) (·) [38] are used in the H-VGRAE. The structure of a three-layer DRN N with exponentially increased dilation factors 1, 2, and 4 is shown in Fig.5. By stacking multiple DRNN layers, the reception fields are extended in different scales. The hidden units at different timestamps are connected based on the layer level, in which the hidden unit h (i) (t − 2 i−1 ) is connected with h (i) (t). By such skip connections, the stacked DRNNs can extract multi-scale temporal features and alleviate the vanishing gradient problem of standard RNN. DRNN (i) (·) can be implemented with two popular variants of RNNs, which are Long Short-Term Memory (LSTM) [35] and Gated Recurrent Units (GRU) [36]. In general, the performance of GRU is as good as LSTM and GRU has fewer parameters and simpler structure. In the case of using GRU, the DRNN (i) (·) for node k is formulated as follows: Step 4: can also leverage the information of the (i + 1)-th latent variable Z (i+1) (t). The up-down structure can pass the high-level information to the low-level feature extraction, to assist the inference local and short-term latent variable with global and long-term information. Then the following inference of Z (i) g (t) is the same as Eq.8.
Step 5: Gaussian Z  According to the normalizing flow technique [16], if Z (i) g,k (t)) for each node k and f (·) is invertible and smooth, the distribution of Z There are various choices of f (·). A popular and simple one is planar NF, in which (10) To make the non-Gaussian Z (i) (t) more variable, the planar NF can be treated as a layer and be stacked for L times g,k (t)))). This extension enhances the representation ability of the latent representation.
2) Decoder: The structure of the decoder is nearly symmetric of that of the encoder, as shown in Fig.4. The inputZ (i) (t) to the decoder is sampled from q(Z (i) (t)|Z (i+1) (t), A(≤t), X(≤t), Z(<t)) for i<M and q(Z (i) (t)|A(≤t), X(≤t), Z(<t)) for i = M . The decoder is similar with that of the encoder is also composed of a content feature extraction block, multi-scale spatial-temporal blocks and density estimators on the probability of A(t) and X(t). Considering the data structure, we choose Bernoulli MLP to reconstruct A(t) and Gaussian MLP to reconstruct X(t).
The procedure of the decoder at time t is as follows.
It is difficult to sample Z (i) (t) directly, due to its complex and implicit distribution. To tackle this problem, we can first sample a valueZ Enc,k (t)I) easily, using the reparametrization trick. Let Z Enc,k (t), where ε ∼ N (0, I). Then, the sampled valueZ g,k (t)).
Step 3: Input C  Dec (t). The computation process is the similar with the decoder in Eq.6.
Step 4: Each S     Step 5: Input the decoded ST features to a FC layer to fuse them, and use the fused features F(t) to reconstruct adjacency matrix A(t). Considering the data structure, we use a Bernoulli MLP to obtain p(A(t)|A(<t), X(≤t), Z(≤t)) and a Gaussian MLP to obtain p(X(t)|A(<t), X(≤t), Z(≤t)).
The Bernoulli MLP utilized is in the simple form: and the Gaussian MLP is:

3) Prior:
As shown in Eq.2, the decoder need to generate a prior distribution of the latent variable. In static scenarios, the prior is commonly set to be N (0, I) [13]. However, as shown in Stochastic RNNs, N (0, I) would be too simple to capture temporal information in the scenario of dynamic networks [12]. Moreover, utilizing complex decoders like RNNs may lead to the KL term in Eq.2 vanishing [39]. This means the latent variable and the input are independent. The encoder q(Z(t)|A(≤t), X(≤t), Z(<t)) = q(Z(t)) and the decoder p(A(t)|A(<t), X(≤t), Z(≤t)) = p(A(t)|X(t)) would just fit the dataset. This makes the training ends to a two-step learning procedure and leads to sub-optimal model.
To avoid such problem, we let the prior p prior (Z(t)|A(<t), X(<t), Z(<t)) come from the prediction based on previous timestamps. Similarly, we first generate the Gaussian form,  prior,g (t) is also extended to non-Gaussian prior Z (i) (t) by using the same planar NF in the encoder. Benefit from the flexible settings of the posterior and prior, the KL vanishing problem would be alleviated.
Similar with the posterior, the prior is also be factorized as In the view of unsupervised representation learning, the predictive prior is an auxiliary objective for the latent variables. Furthermore, due to the different scales of the priors, the prediction process is also hierarchical. Ideally, with the help of the predictive prior, the encoder can infer not only the best latent representation to reconstruct current input, but also the best to forecast future snapshots.
Unfortunately, as illustrated in [37], highly flexible models with deep hierarchies would be difficult to optimize due to multiple conditional stochastic layers. To tackle such problem in H-VGRAE, we aim to find more relations between the encoder and the decoder to help the optimization. Following the information flow direction, Z (i) (t) should be conditioned on Z (i+1) (t) (i<M ). As shown in Fig.6, both the encoder and the prior should infer the distributions of Z (i) (t)|Z (i+1) (t), but obtain q(Z (i) (t)|Z (i+1) (t), A(≤t), X(≤t), Z(<t)) and p prior (Z (i) (t)|Z (i+1) (t), A(<t), X(<t), Z(<t)) respectively. If the encoder and the decoder could share information about the two distributions, the training can be easier to converge.
However, the posterior and prior cannot be explicitly modeled, due to non-Gaussian property by the NF. To achieve the goal of information sharing, we have to trace back to the Gaussian form. Considering the NF, the posterior can be rewritten as Similarly, For the inferred Gaussian distributions Enc,k (t)I) and N (µ prior,k (t)I) We then can employ information sharing mechanism due to their explicit distribution representations.
Inspired by [37], we use the precision-weighted mechanism to introduce prior information into the encoder. The updated mean and variance of q(Z g,k (t), A(≤t), X(≤t), Z(<t)) is then, where k means the k-th node. In such a way, the optimization process would be easier.

C. Online Detection
For an observation of the dynamic network G(t ) at timestamp t , by inputting A(t ) and X(t ) to the trained H-VGRAE, the reconstruction probabilities of edges and node attributes can be obtained. If log p(A ij (t )|A(<t ), X(<t ), Z(≤t ))<α A , the edge for node i to node j is determined as anomalous; otherwise, it is normal.
For detecting anomalous nodes, there are two reasons causing a node to be anomalous, which are anomalous attributes and the edges from or to these nodes. The anomalous nodes caused by the anomalous edges can be detected as above. The anomalous nodes caused by anomalous attributes can be detected by comparing p(x k (t)|A(<t), X(≤t), Z(≤t)) with α X .

D. Anomaly Interpretation
The goal of anomaly interpretation is to interpret what caused an edge or a node to be an anomaly in terms of a comparable and illustrative metirc. Methods in [10] and [11] give illustration of the anomalous edges with heuristically defined anomaly scores. However, these anomaly scores are not normalized, which makes the interpretation and comparison difficult.
In our model, the anomaly score is the reconstruction probability, which is in the range of ([0, 1]). Because each edge and node are conditional independent in the reconstruction, we can rank all the anomaly scores of the edges and nodes, and locate them with low reconstruction probabilities.
Furthermore, to evaluate the anomalous levels of different dynamic networks, we can calculate the anomaly score in different granularity. In the coarse granularity, we can use the proportion of the anomalous edges γ e (t ) = |E a (t )|/|E(t )| and the anomalous nodes γ n (t ) = |V a (t )|/|V (t )| as the anomaly coefficient, where |E a (t )| and |E(t )| are the number of detected anomalous edges and total edges, |V a (t )| and |V (t )| are the number of detected anomalous nodes and total nodes. In the fined granularity, we can calculate the difference between p(A(t )|A(<t ), X(<t ), Z(≤t )) and the distribution p(A thres (t )) with each element p(A thres ij (t )) = α. The anomaly coefficient of edges can be calculated by using KL divergence η e (t ) = KL(p(A(t )|A(<t ), X(<t ), Z(≤t ))||p(A thres (t ))) Similarly, the anomaly coefficient of nodes can be calculated by η n (t ) = KL(p(X(t )|A(<t ), X(<t ), Z(≤t ))||p(X thres (t ))) Then the anomalous levels of different dynamic networks can be determined by η e (t ) and η n (t ).

V. EXPERIMENTS
In this section, we evaluate the performance of the proposed H-VGRAE and compare it to state-of-the-art anomaly detection methods on dynamic networks.
UCI message is a directed network based on an online community of students at the University of California, Irvine. Each node represents a user and each directed edge represents a message interaction between two users. We aggregate all the edges in a day to form a snapshot and then discard the snapshots with few edges. This network does not have any node attribute.
HEP-TH is a dataset created on the all the citations of the papers in High Energy Physics Theory conference from January 1993 to April 2003 (124 months). If a paper i cites paper j, the graph contains a directed edge from i to j. For the first 40 month, we create a citation graph using all the papers published up to that month to be a snapshot. No node attribute is provide for this network.
Social evolution is collected from Jan 2008 to June 30, 2009 and released by MIT Human Dynamics Lab [44]. Each node represents a user and each edge is created by Close Friendship records. We consider the collected information from Jan 2008 until Sep 10, 2008 (i.e. survey date) to form the initial network and use cumulative data for 5 days periods of to form a snapshot of dynamic network (54 snapshots) and the remove nodes with few edges. For this dataset, we consider Calls and SMS records between users as node attributes.
Github is a collected from Jan 2013 to Dec 2013. Each node represents a Github user and each edge is set based on the Follow records. We aggregate 15 days of records to be a snapshot and form 24 snapshots. For this dataset, the Star and Watch records are treated as node attributes.
For UCI message and HEP-TH, which have no node attributes, we consider the N -dimensional identity matrix as node attributes at time t.
2) Baselines: The competing methods compared with H-VGRAE in this paper are summarized as follows.
• GOutlier [4]: GOutlier detects anomalies based on a structural connectivity model. It builds a generative model for edges in a node cluster, and the model can also be used to produce anomalous score for a given edge. GOutlier cannot use node attributes. • CM-Sketch [26]: CM-Sketch defines the local and global structural feature and historical behavior near an edge to measure whether the edge is anomalous or not. It utilizes Count-Min sketch for approximating these properties. CM-Sketch cannot utilize node attributes. • NetWalk [10]: NetWalk first uses a standard autoencoder to generate node embeddings from random walks, and then detects anomalies by clustering the node/edge embeddings. NetWalk cannot handle dynamic networks with node attributes. • AddGraph [11]: AddGraph is also a network embedding based method that designs a graph-RNN to capture ST patterns, and employ a selective negative sampling mechanism and margin loss in training of AddGraph in a semi-supervised fashion. AddGraph can handle dynamic networks with node attributes. • VGRNN [14]: VGRNN is a stochastic network embedding model with single scale ST features and share weights between the encoder and the decoder. VGRNN can utilize node attributes but only maximize the likelihood of edges.
For the models that cannot handle node attributes, only adjacency matrices are used. Moreover, we design ablation experiments to show the effectiveness of the components of H-VGRAE.   The dropout rate is 0.2 and the weight decay parameter for regularization is 1e-5.
For performance comparison, we evaluate different methods based on their ability to correctly detect true and false anomalous edges. To evaluate the robustness of the anomaly detectors, the metric is set to be area under the ROC curve (AUC) score. In all of our experiments, we test the models on the last 10 snapshots of dynamic networks while learning the parameters of the models based on the rest of the snapshots. In the training stage, we use 50% edges of the training snapshots to train the models. Due to the difficulties to obtain ground truth anomalies in the test phase, we inject 1%, 5% and 10% of anomalous edges following the method in [10]. All the experiments are conducted on the deep learning platform, PyTorch [45], on a Linux server.

B. Results and Discussions
Table II summarizes the results for anomalous link detection on different datasets. It is evident that the network embedding based methods (NetWalk, AddGraph, VGRNN, and H-VGRAE) outperform traditional models (GOutlier and CM-Sketch) on all the four datasets. This illustrates the effectiveness of complex pattern captured by network embedding, compared with the artificial features in GOutlier and CM-Sketch. In all the datasets, our H-VGRAE consistently outperforms other competitors by a large margin.
The effectiveness of multi-scale ST features and stochasticity modeling is also shown in Table II. The ST features are crucial for anomaly detection in dynamic networks, which can be shown by comparing the method without extraction of ST features (NetWalk) and the methods with ability of ST modeling (AddGraph, VGRNN and our H-VGRAE). We can see that the ST feature modeling methods have higher AUC scores. Moreover, the robustness achieved from stochastic latent representation is also validated by comparing the deterministic methods (NetWalk and AddGraph) and the stochastic methods (VGRNN and our H-VGRAE). This validates that random variables can carry more information than deterministic hidden states specially for capture complex temporal patterns in dynamic networks. From the perspective of ensemble learning, the multiple sampls of latent variables essentially achieves the ensemble learning in a single model.
To show the effectiveness of multi-scale ST feature extraction, we can compare models in a two steps. Firstly, the single scale version S-VGRAE-I outperforms VGRNN for only near 0.001 in average. The main difference between S-VGRAE-I and VGRNN is that VGRNN shares weights between decoder and the encoder while S-VGRAE-I does not. With more parameters in S-VGRAE-I, the outperformance is not evident, which means that it is not the number of parameters that limits the performance of VGRNN. Secondly, the multi-scale

C. Latent Representation Visualization
To explain how H-VGRAE works for anomaly detection, we visualize the Z-space representations. Because H-VGRAE is a reconstruction based model, the latent Z-space representations directly influence the reconstruction probabilities. During offline model training, H-VGRAE learns the latent representations of normal behaviors of training data. In the online detection stage, H-VGRAE also encodes the input observations according to the normal pattern. If an input snapshot during online detection contains anomalous edges, its Z-space representation should still concerntrate on normal parts, so the reconstruction probabilities of these anomalous edges would be low. Due to the constraint of visualization, we set H-VGRAE(M = 3) with d = 1 and visualize  (Z (1) (t), Z (2) (t), Z (3) (t)) at each node on UCI Message dataset. Figure 9 shows the 3-D Z-space variables learned from the UCI message dataset by H-VGRAE(M = 2). (Z (1) (t), Z (2) (t), Z (3) (t)) are set to their corresponding mean value. We find that the Z-space representation of anomalous nodes highly overlap those normal nodes, indicating that their Z-space representations are quite similar. This phenomenon demonstrates that H-VGRAE learns the normal behaviors well and tries to encodes the anomalous inputs in a normal way. In such a way, the reconstruction is also based on the normal behaviors. Then, the anomalous edges are the inputs deviated from the normal behaviors, which means these edges should not exist according to the normal pattern and the reconstruction probabilities of the anomalous edges are low.

D. Parameter Sensitivity
We examine the influence of hyper-parameters on H-VGRAE, including the number of scales M , the size of dimensions of latent variables d and the training ratio of edges.
To evaluate the influence of M and d, we set the range of M to {1, 2, 3, 4, 5, 6} and the range of d to {4, 8, 16, 32, 64}. Other parameters are set to default. The performance of H-VGRAE is examined on the anomaly detection tasks on the datasets of UCI Message and Github with 5% anomalies. The results are summarized in Figure 10. The AUC score of H-VGRAE increases significantly from M = 1 to M = 2. Further improvement is not evident from M = 2 to M = 4, and the performance degrades after M = 5. Due to the exponentially increase of dilation factor, M = 5 means the skip connection would be 16 timestamps, which is larger than the test snapshots 10. Thus the upper-scale DRNN(·) in M = 5, 6 would degrade to FC layers and loss the ability of modeling temporal dependency. Moreover, with the increase of the number GNN layers, the node representations would be smooth and ends to the same if stacking too many GNN layers [46]. So the upper scales would be noise to interrupt the inference.
With the increase of d, the AUC score of H-VGRAE increase gradually from d = 4 to d = 16, and then drops after d = 32. The reason is that, when d is small, useful information may be missing; while when d is too high, there can be some noise captured by the latent representation. The dimension of the previous layer is 32. When d<32, the overall structure of H-VGRAE would be of the sandwich structure to help the model sufficiently apply bottleneck strategy [47]. The structure can perform the mechanism similar with dimension reduction to preserve critical information.
Then, we evaluete the influence of the training ratio of the edges. The range of the training ratio is set to {10%, 20%, 30%, 40%, 50%, 60%} and other parameters are set to optimum. We use the UCI message dataset and Github dataset with 5% anomalies in test data, and record the AUC scores of the 10 snapshots in the test stage. As shown in Figure  8, with the increase of the training ratio from 10% to 40%, the AUC scores increase sharply, and then the performance stays relatively stable.

VI. CONCLUSION
We have presented a robust semi-supervised model, H-VGRAE, to detect anomalies in dynamic networks by learning hierarchical stochastic network embeddings from multiscale spatial-temporal features. The H-VGRAE can handle scenarios with time-varying edges and node attributes. To break the limits of Gaussian assumption and avoid the KL vanishing problem, the normalizing flow technique and the  predictive prior mechanism is employed. The difficulty of optimizing deep hierarchies of stochastic conditional layers is alleviated by the designed information sharing mechanism between the conditional random variables in the encoder and the decoder. Based on the reconstruction probabilities, the anomalous edges and nodes can be detected with interpretable latent representations. Through extensive experiments, the proposed H-VGRAE outperforms state-of-the-art approaches on four real-world datasets of dynamic networks.
In the future, we will improve H-VGRAE by exploring feature extractors that can tackle the problem of feature degradation with the increase of scales. For the dynamic networks with time-varying number of nodes, new mechanisms of extracting temporal features should be designed to handle variable dimensions of hidden states.