Edge Conditional Node Update Graph Neural Network for Multi-variate Time Series Anomaly Detection

Recently, graph neural networks (GNN) [30] have become a successful method for modeling intricate patterns in data containing graph structures, such as social networks and medical data. A graph comprises a collection of nodes linked by edges that signify the connections between these nodes. GNNs are particularly adept at modeling complex and malleable data structures, surpassing the capabilities of conventional deep neural networks.

Typically, GNNs assume that node representations are constructed based on the influences of neighboring nodes in a graph. Graph convolutional networks (GCN) [31] aggregate the feature representations of nodes by directly considering the feature representations of neighbors.

The attention mechanism, which is a crucial component in numerous deep learning applications, is widely used in sequential models to efficiently capture relevant information. When applied to graphs, this concept gives rise to graph attention networks (GAT) [32] that introduce a weighted aggregation of node representations during the aggregation step in GCN. In GAT, nodes are assigned varying weights through training, prioritizing more important nodes with higher weights and assigning lower weights to less significant ones.

Time-series anomaly detection has been the subject of extensive research for slightly some time, resulting in several proposed methodologies. Classical approaches are characterized by their simplicity and are ill-equipped to handle increasing system complexity. As deep learning has demonstrated exceptional performance in diverse tasks and various domains, numerous methods have emerged to address the challenges of time-series anomaly detection using deep learning techniques.

A fundamental structure for anomaly detection using deep learning is the autoencoder, which trains a network to reconstruct the input data of normal behavior and detects anomalies based on the difference between the input and the reconstructions. USAD [13] is an autoencoder-based model equipped with two decoders. These two decoders are trained adversarially, with one decoder focused on producing accurate reconstructions and the other focused on producing reconstructions that differ significantly from the input data. This dual-decoder setup is designed to effectively detect anomalies that exhibited small deviations from normal data. ACLAE-DT [14] employs an attention-augmented ConvLSTM and a dynamic thresholding mechanism in its reconstruction-based framework. In particular, it constructs as input an inter-correlation feature map that represents the correlations between all pairs of input data and entity embedding vectors, thereby capturing the temporal and contextual dependencies within the data. CAE-AD [15] introduces two contrasting methods: contextual contrasting and instance contrasting. The contextual contrasting method enables the attention-based projection layer to embed shared information between adjacent time steps. The instance contrasting method allows the encoder to capture shared information within a specific time step, utilizing data augmentation in both the time and frequency domains. PAMFE [16] reconstructs clean inputs from noisy ones by using multiple parallel dilated convolutions with different dilation factors. This approach, combined with a feature fusion module, allows the model to capture features at different scales and increases its robustness to noise.

Generative Adversarial Networks (GAN)[33] consist of a generator that generates data and discriminator that distinguishes between real and generated data. Methods using GANs typically train a generator and discriminator on normal time-series data, and compute anomaly scores based on the outputs of the generator or discriminator. MAD-GAN [17] introduces inverse mapping to traditional GANs with LSTM-based generators and discriminators. It uses inverse transformations to compute anomaly scores using both the reconstruction loss and discriminator results. MAD-GAN aims to capture the hidden relationships and interactions that exist between sensors by considering the entire set of sensors simultaneously. FGANomaly [18] generates pseudo-labels indicating whether the input is normal, based on the deviation between the original input and its reconstruction. This model is adversarially trained using data labeled ’normal’ in the pseudo-labels, which increases its robustness to contaminated data. STAD-GAN [19] is designed within a teacher-student framework, combining a GAN structure with a classifier to overcome noise vulnerability and poor performance on anomaly-containing training data. Initially, the teacher model generates a refined dataset with pseudo-labels through unsupervised training. The student model then trains with this dataset, enhancing pseudo-label accuracy. This iterative process improves the classifier’s anomaly detection capabilities. TGAN-AD [20] is a GAN model that includes a transformer-based generator and a discriminator designed to capture contextual information. The generator is trained to identify the latent vector with the highest covariation between the input and generated data, ensuring that the generated data closely resembles the original.

Recent research has extended the success of transformer models [34] in natural language processing to multivariate anomaly detection by adapting them to the time-series domain. TransAnomaly [21] integrates the transformer architecture with a variational autoencoder to tackle the issue of future dependency in sequential data modeling. This approach leads to a model that not only captures future dependencies more effectively but also offers reduced computational complexity, improved parallel processing, and enhanced explainability. GTA [22] extracts graph structures through a learned network and encodes spatial features using GCN. Additionally, it employs transformer models to capture both global and local temporal features. This integrated approach facilitates anomaly detection by considering spatial and temporal dependencies effectively. TranAD [23] is a transformer-based model using two decoders. The first decoder performs reconstruction, whereas the second decoder uses the difference between the results of the first decoder and the input to perform additional reconstruction. Both results contributed to the adversarial training process. D³TN [24] predicts the next time step value using a GRU-based predictor that processes features comprising global spatial features generated by a disentangled convolution layer, local spatial features, and temporal features derived from a self-attention-based layer.

Graph neural networks offer the advantage of explicitly defining the relationships between nodes. Leveraging this capability, they utilize signals from each sensor as node features to learn relationships between sensors, thereby aiding in anomaly detection. MTAD-GAT [25] is a graph attention-based model that constructs a graph structure without prior knowledge. It uses a fully connected graph across both spatial and temporal dimensions. To capture the unique characteristics of these dimensions, MTAD-GAT employs two distinct GAT models: one for spatial features and the other for temporal features. GDN [28] introduces node-embedding vectors to obtain a graph structure based on the similarity of the embedding vectors. Unlike conventional GATs, which rely solely on node features, a GDN concatenates embedding vectors with node features to compute the attention weight. GTAD [26] extracts temporal features using a temporal convolution network and spatial features through a graph attention network. Utilizing these features, the model is trained for data reconstruction and prediction. Anomaly detection is then performed by analyzing both the reconstruction and prediction errors. HAD-MDGAT [27], similar to MTAD-GAT [25], utilizes two GAT models to extract both spatial and temporal features. Leveraging these features, the model employs a GAN for input reconstruction and an MLP for data prediction. MST-GAT [29] captures multimodal correlations by employing both inter-modal and intra-modal attention mechanisms, which leverage modality-specific information. Additionally, it utilizes multi-head attention that operates independently of modality information. PGRGAT [35] employs prior information of sensor relationships to learn graph structure at an edge predictor. This graph structure is then used to encode features via the gated current graph attention unit network, which is a graph-based GRU module.

The dataset consists of multivariate time-series data collected from $N$ sensors over a period of time $T$, $\boldsymbol{s}=[\boldsymbol{s}^{1},\boldsymbol{s}^{2},\dots,\boldsymbol{s}^{T}]\in\mathbb{R}^{N\times T}$. Each column $\boldsymbol{s}^{t}$ at time $t$ corresponds to a vector representing measurements from $N$ sensors, denoted as $\boldsymbol{s}^{t}=[{s}^{t}_{1},\dots,{s}^{t}_{N}]\in\mathbb{R}^{N}$.

Both the training and testing datasets have the same number of sensors, but they differ in their time durations, denoted as $T_{train}$ and $T_{test}$ respectively. The training dataset comprises normal data, aligning with conventional unsupervised anomaly detection methods.

The model proposed in this study utilizes historical data to predict the measured values of the next time step. The input data to predict sensor’s value at time $t$, $x^{(t)}\in\mathbb{R}^{N\times w}$, is constructed using sliding window with window size $w$ over previous time series data, as follows:

After the prediction, our model calculates an anomaly score for that specific time step. If the anomaly score exceeds a predefined threshold, the model outputs 1 to indicate an anomaly and 0 if the data is considered normal.

Our model extracts the graph structure from node embedding vectors, predicts sensor values for the next time step, and detects anomalies. To achieve this, our model is composed of the following four main components:

1. 1.

   Node Condition Embedding uses embedding vectors to be utilized as conditions to transform source node representations and readout final node representations, and to construct the graphs structure.

2. 2.

   Graph Structure Extraction generates graph structure based on the similarity between all pairs of the node embedding vectors.

3. 3.

   Condition-based prediction forecasts each sensor value for the next step by employing node embedding vectors as conditions in the node conditional readout module and edge conditional node update module.

4. 4.

   Anomaly Detection calculates anomaly scores, if the score is higher than a threshold, the model determines the time step as an anomaly.

Fig. 1 shows an overview of our model.

We introduce node embedding vectors for getting relationships between sensors, as $d_{e}$ dimensional vectors:

where $v_{i}$ represents the $i$th node embedding vector.

In [28], the sensor-embedding vectors are trained to reflect the similar behavior of the sensors in graph attention. However, graph attention can cause node representations to become similar, which may contribute to the model’s limited capability in modeling sensor behavior.

To avoid this problem, our model trains the node embedding vector without attention, instead, employing the edge conditional node update module and node conditional readout module, as described in the following section.

To discover the relationships between sensors in systems where these relationships are not predefined, we employ the method of assessing similarities among all node embedding vectors. This similarity measurement is crucial for constructing the graph structure of the system. Specifically, for each node, we determine the Top-$k$ most similar nodes, based on the similarity score $e_{i,j}$, to form an adjacency matrix:

where TopK is the operator used to select the indices of Top-$k$ values among the set of inputs. The value $k$ of TopK is given by the user to adjust the sparsity level of the graph structure.

Condition based prediction process comprises two main modules: the edge conditional node update module (ECNUM) and the node conditional readout module (NCRM). ECNUM uses a graph-based method to generate a node representation, which is then utilized by NCRM to predict the final result. In this paper, we use MLP as ECNUM and NCRM.

A simple encoder network such as a linear or 1D CNN could be adopted to encode the input to $d_{f}$ dimensional vector, $z^{(t)}\in\mathbb{R}^{N\times d_{f}}$, from a given input $x$, as formulated in Eq. 5.

Subsequently, each row of $z^{(t)}$ is adopted as the node representation. In this paper, we utilize a linear layer as an encoder.

Given the values predicted from the learned graph structure, we detect anomalies using graph deviation scoring (GDS) from [28]. GDS calculates the anomaly score by comparing each sensor’s scores, using the absolute error as the basis for normalization.

The first step in calculating the GDS is to calculate the absolute error between the predicted value and the actual value at time $t$ and sensor $i$:

The number of absolute errors at time $t$ is $N$ number of sensors. It is not appropriate to compare these absolute errors directly since the absolute error of each sensor exhibits a different behavior. Therefore, we normalize each absolute error as follows:

where $\tilde{\mu}_{i}$ and $\tilde{\sigma}_{i}$ are the median and interquartile range (IQR) of the absolute error at sensor $i$. IQR is defined as the difference between the 1st and 3rd quartiles of set of values. This normalization method is known to be more robust against anomalies than normalization using mean and standard variation.

To compute the anomaly score at given time $t$, we use the max function to select the most anomalous sensor and score as follow:

The model cannot perfectly predict the sensor value; as a result, the anomaly score can exhibit a sharp change, even when the input is normal. To mitigate these changes, we use a simple moving average (SMA).

Finally, if the smoothed anomaly score is above a fixed threshold, the model labels the input as abnormal for a given time, $t$.

Many previous methods use automatic thresholding methods, such as extreme value theory [37], to determine the threshold value. However, these methods produce different results depending on the initial set of parameters. To evaluate the anomaly detection capability of the model, we opt for a grid search method to identify an optimal threshold that maximizes the F1 score.

The pseudo-code of the edge conditional node update GNN is presented in Algorithm 1

Our model is evaluated using three real-world datasets. Two of these datasets, SWaT [38] and WADI [39], are sourced from water treatment system-based test beds. These datasets are generated by simulating operator-defined attack scenarios, with evidence labels provided for each scenario. The third dataset, PSM [40], is collected from real-world server machines and is manually labeled by engineers and application experts. The PSM dataset includes both planned and unplanned anomalies

The SWaT dataset is constructed using a water treatment test bed that emulates a small-scale system. It contains 51 features, 47,515 training samples, and 44,986 test samples, with an anomaly rate of approximately 11.97%.

The WADI dataset is built based on a water distribution test bed, which models a larger-scale system than SWaT. This dataset containes 127 features, 118,795 training samples, and 17,275 test samples, with an anomaly rate of approximately 5.99%.

The PSM dataset is collected internally from multiple application server nodes at eBay. It contains 25 features, 132,481 training samples, and 87,841 test samples, with an anomaly rate of approximately 27.75%

The SWaT and WADI datasets are downsampled following the methods outlined in GDN [28]. Data are downsampled every 10s using median values to speed up the training and the labels for the test dataset are determined based on the predominant label among the down-sampled labels. The missing data is replaced by mean value. As noted by [41], data in the first 5–6 hours after system start-up is considered unstable; therefore, we exclude the first 2,160 data points.

We implemented our model using PyTorch version 1.13.0, CUDA 11.6, PyTorch Geometric library 2.2.0, and PyTorch Lightning 1.9.0. A server equipped with an Intel(R) Xeon(R) Gold 6148 2.4 Hz CPU and an NVIDIA RTX A6000 GPU is used to train the model. The model is trained using Adam optimizer with a learning rate of 0.001 and beta values of 0.9 and 0.99. We employ an early stopping strategy with a maximum of 100 epochs and patience of 5 and select the final model with the lowest validation loss.

To identify suitable hyperparameters for the model, we employ the grid search method. The final hyperparameters are presented in Table 1.

We conduct a comprehensive comparison of our model with several baseline models. The comparison includes MAD-GAN [17] from GAN-based models, USAD [13] representing autoencoder-based models, MTAD-GAT [25] and GDN [28] from GNN-based models, as well as TranAD [23] from transformer-based models, and GTA [22] from a hybrid model between transformer and GNN. Since the point adjustment strategy proposed by [42] has the great possibility of overestimating the model’s performance, as noted in [43], we omit this strategy when calculating the performance.

As evaluation metrics, we use metrics commonly used in previous research, including F1 score (F1), recall (Rec), and precision (Prec).

The precision measures the ratio of correctly detected anomalies to the total number of detected anomalies. This is calculated by dividing the number of correctly detected anomalies by the total number of detected anomalies.

The recall measures the ratio of correctly detected anomalies to the total number of actual anomalies in the dataset. This is calculated by dividing the number of correctly detected anomalies by the total number of actual anomalies.

The F1 score is the harmonic mean of the precision and recall, providing a balanced measure that considers both in a single metric. The precision, recall, and F1 scores can be formulated as follows:

where $TP$, $FP$, and $FN$ denote the numbers of true positives, false positives, and false negatives, respectively.

Table 2 presents a comprehensive overview of the average performance and standard deviation of both our models and the baseline models, derived from a total of 10 experiments. Our model shows performance with an average F1 score of 0.8089 on the SWaT dataset, 0.5730 on the WADI dataset, and 0.6027 on the PSM dataset. Notably, these scores represent a significant improvement over the best-performing baseline models, with increases of 4.0% for the SWaT dataset, 36.2% for the WADI dataset, and 3.9% for the PSM dataset, respectively.

Table 3 shows the performance of the models with the best F1 score. Our model shows performances with an F1 score of 0.8541 on the SWaT dataset, 0.6405 on the WADI dataset, and 0.6179 on the PSM dataset. When comparing our model to the best F1 baseline models, we observed a significant performance improvement. Specifically, our model achieved remarkable increases of 5.4%, 12.4%, and 6.0% in F1 scores for the SWaT, WADI, and PSM datasets, respectively.

These results consistently demonstrate a notably high performance improvement in the complex WADI dataset, with a 36.2% increase in average F1 score and a 12.4% increase in the best F1 score compared to baseline models. Moreover, the comparison with the graph-based baseline GDN model, which is similar to our approach, further validates the benefits of our model’s strategy to use distinct source node representations for each target node.

To evaluate the effectiveness of our model in detecting and localizing anomalies, we analyzed the second anomaly in the WADI dataset. This anomaly is triggered by an attack that manipulates the flow sensor 1_FIT_001_PV, causing it to report values different from the actual measurements. During the attack, the values remain within the normal sensor measurement range. As shown in Fig. 6, it is almost infeasible to detect anomalies using the sensor independently.

However, during this attack period, the proposed model correctly identified 1_MV_001_STATUS as a sensor with a high anomaly score, as shown in Fig. 7. Traditional models typically rely on attention weights to identify sensors associated with 1_MV_001_STATUS. In contrast, our model uses layer-wise relevance propagation (LRP)[44] instead of attention mechanisms to provide information about relevant sensors.

If we apply the standard LRP to our model, the edge conditional node update and the node conditional readout module do not propagate the relevance scores corresponding to the node-embedding vectors. To address this problem, we reassigned the relevance scores of the node-embedding vectors to the feature relevance scores. For the node conditional readout module, we reassigned the relevance score of the one-node embedding vector to the feature relevance score, as follows:

where, $R^{x_{j}}_{i}$, $R^{v_{j}}_{i}$ are $i$th element of relevance score of feature $x_{j}$, and node embedding of node $v_{j}$, respectively.

The edge conditional node update module utilizes the summation of the relevance scores of the two node embedding vectors. The reassignment considering redundancy is conducted as follows.

where $|\mathcal{N}(i)\cup\{i\}|$ is the number of edges connected to node $i$, $R^{S_{i,j}}$ is the relevance score vector for source node embedding at edge from node $i$ to node $j$, $R^{T_{j,i}}$ is the relevance score vector for target node embedding at edge from node $j$ to node $i$.

Due to the nature of the data, the inputs of some sensors may become zero. In this case, the relevance scores could also become zero, making it difficult to make accurate comparisons. To address this issue, we summed the relevance scores for each node from the relevance scores that correspond to the node representation immediately after input transformation. Since every input value is transformed independently into a node representation, the total relevance scores for each node remain consistent.

The graph structure learned by our model and LRP results centered on the 1_MV_001_STATUS sensor are shown in Fig. 8. In the figure, the red square represents 1_MV_001_STATUS and the red triangle represents 1_FIT_001_PV. The red and blue edges denote the positive and negative relevance scores, respectively, and the edge thickness represents their relative magnitudes. This result confirms the high relevance of 1_MV_001_STATUS and 1_FIT_001_PV.

We assess the sensitivity of our model with respect to the window length and Top-$k$ factor. We report the average F1 score and standard deviation based on 10 repetitions of the experiments. For window length, experiments are conducted with window lengths of 3, 5, 7, 10, 15, 20, and 30. Regarding Top-$k$, experiments are carried out with values of 5, 10, 15, 20, 25, 30, 35, and 40. However, for the PSM dataset, experiments are conducted only up to 25, as the dataset has a number of features lower than 30.

Figure 9 illustrates the results of the sensitivity analysis to the sliding window length. The figure shows that for the SWaT and WADI datasets, optimal performance is achieved with a window length of 5. In contrast, the PSM dataset exhibits stable performance across varying window lengths.

Figure 10 shows the results of the sensitivity analysis for the Top-$k$ factor. The figure shows that the performance of the models initially increases, reaches an optimal point, and then decreases. It is noteworthy that for the PSM dataset, there is no observed drop in performance, as the optimal point is identified at 25.