GRELEN: Multivariate Time Series Anomaly Detection from the Perspective of Graph Relational Learning

1 Introduction

Time series anomaly detection is an important research topic and has wide applications. For example, in industry, sensors are mounted in the system for production line monitoring. Generally, a system has different sensors to describe its global state and requires multivariate time series anomaly detection techniques to trigger system-level alarm [Su et al., 2019].

As sensor dimension increases, dependence relationship between sensors becomes more important for efficient anomaly detection. For example, Fig.1 plots time series signals of four sensors in a secure water treatment system (SWaT). It shows one anomaly happened on sensor LIT101. Then sensor FIT-101 and actuator MV-101 follow the anomaly and drop to low level values, while sensor P-101 remains uninfluenced. Yet from the detection perspective, it does not mean sensor P-101 is useless, since its dependence relationship between others is actually changed and will help detect anomaly. Similarly, the relationship between attacked sensor and other unattacked sensors will also be the indicator of anomaly. It motivates that monitoring dependence relationship would be more powerful in complex systems, compared with merely focusing on monitoring each individual sensor's behavior.

∗Corresponding author.

So far numerous works have been proposed for multivariate time series anomaly detection. Neural networks are widely used given high data dimension and generally have better performance than traditional statistical methods. Considering the anomalies are generally unpredictable with various patterns, unsupervised-learning-based methods are more attractive. The basic framework is to use no-anomaly data to establish a model describing the normal system pattern. Then testing data that cannot be well fitted by the pre-trained model would be regarded as anomalies.

However, most existing methods cannot model the anomaly from the perspective of multivariate time series dependence relationship. Though there are also some graph learning based methods considering dependence relationship in modeling, their intrinsic information loss will consequently undermine the detection power, which will be illustrated in Section 4.6.

In this work, we propose a novel anomaly detection framework, called Graph Relational Learning Network (GReLeN), to fill the gaps mentioned before. Our contributions could be summarized as follows:

• (a) GReLeN utilizes VAE structure to learn a probabilistic relation graph for multiple sensors, where the latent vari-

• able of VAE is used to capture the dependence relationship between sensors.
• (b) We propose to detect anomaly from the perspective of graph relational learning for the first time. The anomaly detection score is established based on the explicit dependence relationship learned from model.
• (c) The experiments on four real world datasets show that our proposed model could describe the latent dependence among sensors well and has good interpretability corresponding to domain knowledge. Our proposed model also shows superiority regarding detection accuracy and anomaly diagnosis.

2 RelatedWork

2.1 Multivariate Time Series Anomaly Detection

Considering the quick development of neural network, only deep-learning based models are reviewed in this subsection. The current works could be classified as prediction-based methods and reconstruction-based methods.

(a) Prediction-based models: Prediction-based models utilize advanced deep learning components to capture spatialtemporal dependence, and use the prediction error as anomaly score for detection. Well designed modules could achieve more accurate prediction, and consequently be helpful to detect abnormal ones. [Deng and Hooi, 2021] proposes a GNN based method to aggregate the information between sensors. [Zhao et al., 2020] combines feature-oriented Graph Attention Network (GAT) and time-oriented GAT to handle spatial dependence and temporal dependence while predicting.

(b) Reconstruction-based models: Reconstruction-based methods hope to find a latent representation for the entire time series for data reconstruction. Loss function of model is a common choice for anomaly score. [Li et al., 2019] uses Long Short-Term Memory (LSTM) as basic cells, and considers the entire variable set concurrently with a Generative Adversarial Network (GAN) framework. In [Su et al., 2019], the proposed OmniAnomaly uses stochastic Recurrent Neural Network (RNN) to find robust representations for multivariate time series. [Audibert et al., 2020] proposes an autoencoder architecture whose adversarial-style learning is inspired by GAN. Recent work [Abdulaal et al., 2021] utilizes spectral analysis on the latent representation and produces a synchronized representation for multivariate data. However, in these works, no dependence relationship between variables has been considered and modeled explicitly.

2.2 Graph Learning inMultivariate Time Series

GNN is commonly used in multivariate time series prediction, where they regard each single time series as a node and their dependence as edges in graph. A common assumption is that the graph structures should be pre-defined based on domain knowledge. It makes the model very sensitive to the choice of graph, and loses generality. Recent studies adopt a data-driven way (i.e. graph learning strategy) to learn the graph structure between different time series automatically. In [Shang et al., 2020; Wu et al., 2020], graph structure learning is used for improving multivariate time series prediction. In [Kipf et al., 2018], a Neural Relational Inference

(NRI) model is used to learn the system dynamics from observational data based on VAE. However, their node-to-edge and edge-to-node operations on graph result in heavy computation complexity, which hinders their application in highdimensional time series.

[Deng and Hooi, 2021] is the most similar work to ours. It uses node embedding to learn graph structure and do prediction-based detection. However, its learning strategy assumes the graph structure should be determined in advance instead of automatically learned, which may lead to low model robustness. Besides, it assumes the graph structure is sparse and unchanged over time, which leads to information loss. Furthermore, its learned graph structure only acts as an intermediate module to improve prediction, but is not used to indicate any anomaly explicitly. More comparisons will be explained in Section 4.6.

3 ProposedModel

3.1 ProblemFormulation

Consider a system with N sensors, our sequential data collected over T timestamps could be denoted as X ∈ RN×T . We denote the data collected at each sensor i and each timestamp t by using subscript and superscript respectively. That is, Xi = Xi,:,Xt = X:,t. In our unsupervised assumption, the system operates in normal condition during the first Ttrain timestamps, which could be regarded as training data. Our task here is to identify potential anomaly in the following timestamps t > Ttrain. All the data coming after Ttrain will be regarded as testing set.

A sliding window w is used for constructing the samples: St = Xt-w+1:t . The task of our anomaly detection model f is to provide a set of binary labels indicating whether there's any anomaly for certain timestamp of the testing set: yt = f(St ), yt ∈ {0, 1} , 1 ≤ t ≤ Ttest. The mathematical notations are shown in Table 1.

3.2 Anomaly Detection Process

As shown in Fig.2, our proposed anomaly detection structure is composed of an offline training module and an online testing module. The offline training phase uses normal training data to learn the normal sensor dependence relationship together with the data reconstruction mechanism. Normal pattern of time-varying multivariate system and the pretrained model could be obtained. In the online testing phase, the model is used for testing data stream without additional training. By monitoring the sensor dependence relationship of the testing data, anomaly detection and diagnosis could be achieved.

3.3 Model Overview

The overall architecture of our model GReLeN is shown in Fig.3. GReLeN follows the main structure of VAE. In particular, three sub-modules are contained in our proposed framework.

• (a) Encoder: In the encoder part, we use linear layers to extract temporal features for each time series input. Selfattention-based [Vaswani et al., 2017] mechanism would be applied for efficient high-level feature capturing and dependence relationship inference.
• (b) Latent variable with sampling: The output of encoder is used to parameterize the distribution of dependence relationship among different time series. Reparametrization trick would be used for enabling backpropagation.
• (c) Decoder: In the decoder part, our goal is to reconstruct the input series recurrently. A spatial-temporal cell would be applied here based on the relational graph we have sampled by latent variable.

3.4 Encoder

The input of encoder is a sample series St generated by sliding window w. The target of encoder is to learn the distribution qϕ that the latent dependence relationship Z t |St follows. In particular, we can denote that qϕ is parameterized by Θt = {θt i,j , i, j = 1, . . . , N}.

We first use a linear layer to extract temporal feature Ht from the raw input series, where Wl ∈ RT×c1 is the learnable parameters and c1 is the dimension of high-level extracted temporal features: Ht = StWl.

Self-attention-based mechanism would be applied to calculate Θt based on Ht = {hti , i = 1, . . . ,N}. More specific, high-level extracted feature will be projected into latent query, key subspace Qt ,Kt with learnable parameters Wq,Wk. Then the dependence relationship distribution parameters can be obtained by dot-product attention:

Further, it is natural to extend the multi-head attention strategy here to learn multiple types of dependence and improve the representation ability of latent variables. For the attention with h0 heads, we can get Θt as follows:

Θt = Softmax (Concat (head1, ..., headh0 )) = [θt i,j ]1≤i,j≤N Qth = HtWq,h, Kth = HtWk,h, headh = QthKth T , 1 ≤ h ≤ h0 where Wq,h and Wk,h are learnable parameters of the hth

head, Θt ∈ Rh0×N×N .

Here, each pair-wise parameter θt i,j is a vector with length h0, since it means we have h0 types of dependence in total. It's intuitive that we assume θt i,j is the parameter for categorical distribution (e.g. multinomial distribution). Therefore, we use the softmax function to guarantee the sum of the probability of all these types equals 1.

3.5 Sampling

The above encoder formulates the distribution qϕ to generate latent variable zt ij ∈ Rh0 , which is supposed to be a onehot vector indicating which type of dependence it belongs to. However, the discretization of zt ij hampers backpropagation and brings additional trouble in the training process.

A recent solution is to use Gumbel-Softmax categorical reparameterization trick [Jang et al., 2017]. It utilizes the continuous variable approximation to make backpropagation possible by

where g is a h0 dimension vector of i.i.d. samples drawn from a Gumbel(0,1) distribution and the continuous distribution converges to expected one-hot samples if we have temperature parameter τ → 0.

3.6 Decoder

The decoder is designed to reconstruct the original input series with the guidance of learned dependence relationship. Recall that zt i,j is an h0 dimension vector. Furthermore, Z t ∈ Rh0×N×N could be considered as the concatenation of all the h0 graph adjacency matrices. We model the reconstruction module by following a Gaussian distribution assumption:

for t′ = t - w + 1, . . . , t - 1. We divide the h0 types of learned dependence into two parts. One of the graphs will be treated as "null dependence", which aims to guarantee the sparsity via prior setting. The other h0 - 1 graphs denote different relationship structures, which will be used as graph adjacency matrix in Diffusion Convolutional Gated Recurrent Unit (DCGRU)[Li et al., 2018] to do the recurrent reconstruction. All the outputs from different graphs will be summed up to obtain the final output.

Mathematically and formally, for any graph dependence structure with adjacency matrix A, DCGRU combines GRU temporal cell and diffusion convolutional GNN:

where the graph convolution operation ⋆A is defined as:

Here, DO and DI are the out-degree and in-degree matrix of graph. ∥ denotes concatenation. For Q = R,U,C, wQ k,1 , wQ k,2 , bQ are learnable parameters and K is the hyperparameter reprensenting the diffusion degree while convolution. At each time step t′ , the hidden state Ht′ g serves as our final reconstruction µt′+1 .

3.7 Objective Function and Training

As a VAE-based model, the objective function could be maximizing the evidence lower bound (ELBO) of VAE:

The first term Eqϕ(Zt |St ) [logpψ (St |Z t )] is the so-called reconstruction loss and KL [qϕ (Z t |St ) ∥pψ (Z t )] indicates the Kullback-Leibler (KL) divergence loss.

In the decoder, we could decompose the first term in a recurrent way:

Combining Equation (1) and (2), for each sample St , the total reconstruction loss could be estimated by

The second term in L, i.e., the KL divergence loss, measures how well the dependence relationship we learn from encoder matches the prior distribution.

Given each sensor has influence or could be influenced by only a small proportion of others, we would assume that the dependence relationship is sparse among all these sensors. Therefore, the type denoting "null dependence" should have a larger prior. Among all these h0 graphs, only h0 -1 graphs will be used in decoder for reconstruction. The graph which has not been used while decoding would be regarded as "null structure". If two series have connections in "null structure", they will not have dependence in practical because of the discarding.

Denote our prior multinomial distribution parameters as ˜ q = (˜ q1, ˜ q2, ..., ˜ qh0 ). Mathematically, the KL divergence loss could be computed by

where ˜ ptk = ∑Ni=1 ∑Nj=1 zt i,j (k) is the sum of probability that the learned dependence between sensor i and sensor j belongs to type k.

4 Experiments

4.1 Datasets

We conduct our proposed model on four real-word datasets: SWaT (Secure Water Treatment Testbed) [Goh et al., 2016], WADI (Water Distribution Testbed) [Ahmed et al., 2017], SMD (Server Machine Dataset) [Su et al., 2019], and PSM (Pooled Server Metrics) [Abdulaal et al., 2021]. Normal data would be divided into training data (80%) and validation data (20%). Training data is used for model training, and validation helps model selection. Anomaly only lies in testing data. More descriptions about datasets are shown in Table 2.

In particular, SWaT, WADI and PSM have only one entity to be monitored. SMD has 28 different entities, where each entity has 38 sensors. Different entities need to be trained and tested independently.

4.2 Baselines and EvaluationMetrics

To illustrate the superiority and effectiveness of our proposed model, we compare our model with the following traditional methods IF: Isolated Forest [Liu et al., 2008]; AAD: Active Anomaly Discovery [Das et al., 2016], and deep-learningbased state-of-the-art baselines: LSTM-VAE [Park et al., 2018]; MAD-GAN [Li et al., 2019]; OmniAnomaly [Su et al., 2019]; USAD [Audibert et al., 2020]; GDN [Deng and Hooi, 2021]; RCoders [Abdulaal et al., 2021].

The evaluation metrics we consider here include Precision, Recall and F1 Score (F1). A commonly used point adjustment strategy [Su et al., 2019; Deng and Hooi, 2021; Audibert et al., 2020] is involved to assure that if at any timestamp of an anomaly's occurring period the anomaly is detected, we regard it as an accurate detection.

4.3 Experimental Setup

In our experiments, we set c1 = 64, h0 = 4 and all the hidden dimension in DCGRU cell equal to 64. The number of DCRNN layers is 2 in SWaT and 1 in other datasets. The batch size is 32. The method is implemented by Pytorch using the Adam [Kingma and Ba, 2015] optimizer with the learning rate 1 × 10-3 . All the samples are generated with a sliding window w = 30. 100 epochs are used for training. Our prior multinomial distribution parameters setting is ˜ q = [0.91, 0.03, 0.03, 0.03].

4.4 Anomaly Score Establishment

We compare two ways of anomaly score establishment in this work. First we utilize the loss function, a common type of anomaly score regarding reconstruction-based methods. More specifically, we use KL divergence loss L2 here (denoted by GReLeN Loss). By experiments we find that L2 is more stable and comprehensive than reconstruction loss L1.

However, the KL divergence loss mainly focuses on global deviation of the dependence structure from its prior distribution. As mentioned before, the specific dependence relationship between sensors also reveals the anomaly and includes more local information. Therefore, in dependence graphs, we use the sum of all the sensors' in/out degree as another anomaly score (denoted by GReLeN Degree). Consider only sudden changes should be detected as anomaly, while the gradual shift is acceptable for normal operation, we use a moving average filter strategy for noise smoothing. By experiments, we find that the sudden, anomalous changes could be enlarged and the gradual, noisy shift could be removed via this strategy.

Mathematically we have, GReLeN Loss = L2, GReLeN Degree = ∑Ni=1 din i + dout i . Here, din i , dout i are moving filtered in/out degree of node i in the learned graph respectively. In our composite score, we consider the bi-direction dependence changes for all the time series. We believe that our anomaly score could provide a comprehensive metric to reveal the changes in the state of system.

4.5 Results

We evaluate our proposed model with 5 runs. Consistent with [Audibert et al., 2020; Abdulaal et al., 2021; Su et al., 2019], we use grid search to get possible anomaly thresholds for every model and report the results with the highest F1 score.

As shown in Table 3, our proposed GReLeN model achieves the best F1 performance on average, with the highest F1 on three datasets and comparable results on SMD. Comparing with other SOTA baselines, GReLeN has better and more balanced performance between Precision and Recall. Furthermore, WADI has low anomaly rate and consequently bring large challenges to other works. Yet our GReLeN could still achieve a high Recall performance with low missing alarm. Also, our model achieves balanced performance on all these datasets, which illustrates the model robustness.

Unsatisfactory performance in other reconstructed-based models (LSTM-VAE, MAD-GAN, OmniAnomaly, USAD) highlights our superiority to involve explicit dependence relationship inference in data reconstruction. As to predictionbased GDN, it uses node embedding and GAT to model node correlation, and improves the prediction performance. Yet it has intrinsic limitation due to only using naive prediction error to be anomaly detection metric. For our GReLeN models, though GReLeN Loss could detect potential anomaly, its global loss based anomaly score is too coarse to detect minor anomaly. In contrast, our designed dependence relationship structure based anomaly score (GReLeN Degree) is more powerful to describe and capture system anomaly.

4.6 Effect of Graph Relational Learning

We further demonstrate the superiority of our dependence relationship learning strategy compared with previous similar work GDN [Deng and Hooi, 2021]. They claim to use explicit graph to capture spatial correlation, but only "top-k" neighbors have be retained. Consequently its GAT can only provide dynamic weights to the retained neighbors. In this way, if two nodes' dependence relationship is not that intense and get discarded while training, they cannot provide further information in testing phase.

In contrast, our proposed GReLeN will retain all the potential neighbors, and use a proper prior to help control the sparsity of learned graph in a data-driven way. Hence it can avoid the information loss occurred in "top-k" sparsity.

We validate the above discussions by a comparison, where we retain only "top-k" neighbors for GReLeN. The sparsity parameter k is 15, 30, 10, 8, which is consistent with the settings in GDN experiments. As shown in Table 3, comparing GReLeN Degree and GReLeN Topk, a better F1 and Recall performance comes from our GReLeN Degree model.

Meanwhile, as mentioned in Section 3.5, we use the Gumbel-Softmax categorical reparameterization trick to generate the discrete latent variables. Fig.4 shows the reparameterized probability of the dependence relationship not belonging to "null structure". The middle figure of Fig.4 shows the expected results of the prior, while the other two plots show the latent variables learned in abnormal and normal cases. We find that the anomaly could be detected by comparing the latent variables. First, since in the normal case, we assume all the dependence relationships of nodes are sparse, hence we can see besides the "null structure", all the other learnt graph structures are quite sparse. However, in the abnormal case, we can observe that the dependence structure of nodes is no longer sparse due to the right tail in the top figure, which means the system with anomaly tends to have more intensive dependencies.

4.7 Model Interpretability and Anomaly Diagnosis

The dependence structure we learned by latent variable makes our proposed model have better interpretability. In this subsection, we discuss the relational learning and anomaly diagnosis ability.

Fig.5(a) shows the relational graph adjacency matrix learned on SWaT training data, which could be partly interpreted by human knowledge. Three types of graphs are involved to capture different views of dependence relationships. According to the second type of graph, one may conclude the first four sensors have strong dependence with other sensors. Referring to the process diagram of this water treatment system [Goh et al., 2016], the first four sensors consist of the first treatment process phase P1 and should be considered as the fundamental elements of system, which leads to much closer relationship with other sensors. The third type of graph indicates that a small group of sensors have strong dependence among each other. In fact, the sensors circled in purple belong to the same treatment process phase P3. It is also acceptable that these sensors have strong dependence because of their same process phase. In Fig.5(b), the system in normal condition shows similar graph pattern to what we have learned in the training phase, while anomaly may lead to a totally different structure (Fig.5(c)).

Our model is also helpful for model diagnosis. A closer look at Fig.5(c) finds two significant abnormal hubs as highlighted by rectangles. The abnormal hub in block1 shows an actuator has an anomalous setting, which has been recorded in the operation log. The abnormal hubs in block2 lies in process stage P5, i.e., the downstream process of the attacked point in block1. It indicates that our method can not only detect the recorded anomaly, but could also report the relative anomalies happening on downstream process beyond the operation log.

Fig.6 shows a case study for WADI. The WADI system has three sub-processes. The sub-process that each sensor belongs to is marked by the first number of the sensor's name. The two most suspicious sensors, whose anomaly scores are exhibited in Fig.6(c)(d), are consistent with the ground-truth operation log. The dependence graph learned in normal condition (Fig.6(a)) has sparse connections for sensors. But in the anomalous state, the attacked sensors(shown in orange color) have much intenser connections with other sensors in the same sub-process, which is consistent with our domain knowledge.

The above case studies show our model is powerful to diagnose anomaly, and can provide crucial information and guidance for daily operation.

5 Conclusion

In our work, we proposed to detect anomaly in multivariate time series from the perspective of graph relational learning. Our Graph Relational Learning Network (GReLeN) model combines VAE structure as well as a graph dependence structure learning strategy for anomaly detection in multivariate time series in a reconstructed way. The latent variable in VAE captures the dependence relationships between sensors explicitly, based on which a well-designed anomaly score is constructed. The experiments show our superiority over both prediction-based and reconstruction-based methods. The interpretability helps us with anomaly diagnosis and provides in-time guidance to daily operation.

Acknowledgements

This work was supported by the RGC GRF 16216119, Foshan HKUST Projects FSUST20-FYTRI03B, in part by the NSFC Grant 71901131, 71932006 and in part by the BNSF Grant 9222014.