GRAPH-AUGMENTED NORMALIZING FLOWS FOR ANOMALY DETECTION OF MULTIPLE TIME SERIES
Year: 2022

Abstract:
Anomaly detection is a widely studied task for a broad variety of data types; among them, multiple time series appear frequently in applications, including for example, power grids and traffic networks. Detecting anomalies for multiple time series, however, is a challenging subject, owing to the intricate interdependencies among the constituent series. We hypothesize that anomalies occur in low density regions of a distribution and explore the use of normalizing flows for unsupervised anomaly detection, because of their superior quality in density estimation. Moreover, we propose a novel flow model by imposing a Bayesian network among constituent series. A Bayesian network is a directed acyclic graph (DAG) that models causal relationships; it factorizes the joint probability of the series into the product of easy-to-evaluate conditional probabilities. We call such a graph-augmented normalizing flow approach GANF and propose joint estimation of the DAG with flow parameters. We conduct extensive experiments on real-world datasets and demonstrate the effectiveness of GANF for density estimation, anomaly detection, and identification of time series distribution drift.

heading: 1 INTRODUCTION
text: Anomaly detection (Pimentel et al., 2014; Ruff et al., 2018) is the task of identifying unusual samples that significantly deviate from the majority of the data instances. It is applied in a broad variety of domains, including risk management (Aven, 2016), video surveillance (Kiran et al., 2018), adversarial example detection (Grosse et al., 2017), and fraud detection (Roy & George, 2017). Representative classical methods for anomaly detection are one-class support vector machines (Schölkopf et al., 2001) and kernel density estimation (Parzen, 1962; Kim & Scott, 2012). These methods rely on handcrafted features and often are not robust for high-dimensional data (e.g., images, speech signals, and time series). In recent years, inspired by the success of deep learning for complex data, many deep anomaly detection methods were proposed and they are remarkably effective in applications (Ruff et al., 2018; Sabokrou et al., 2018; Goyal et al., 2020).
Apart from these demonstrated applications, increasing demand exists for the anomaly detection of even more complex data; i.e., multiple time series. They contain a set of multivariate time series that often interact with each other in a system. A prominent example of the source of multiple time series is the power grid, where each constituent series is the grid state over time, recorded by a sensor deployed at a certain geographic location. The grid state includes many attributes; e.g., current magnitude and angle, voltage magnitude and angle, and frequency. Time series readings from sensors at nearby locations are often correlated and their behavior may be causal under cascading effects. Anomaly detection amounts to timely identifying abnormal grid conditions such as generator trip and insulator damage.
Anomaly detection of multiple time series is rather challenging, due to high dimensionality, interdependency, and label scarcity. First, a straightforward approach is to concatenate the constituent series along the attribute dimension and apply a detection method for multivariate time series. However, when the system contains many constituents, the resulting data suffer high dimensionality. Second, constituent series bear intricate interdependencies, which may be implicit and challenging to model. When an explicit graph topology is known, graph neural networks are widely used to digest the rela-
∗This work was done while E. Dai was an intern at MIT-IBM Watson AI Lab, IBM Research. †To whom correspondence should be addressed.
tional information (Seo et al., 2016; Li et al., 2018b; Yu et al., 2018; Zhao et al., 2019). However, a graph may not always be known (because, for example, it is sensitive information) and hence graph structure learning becomes an indispensable component of the solution (Kipf et al., 2018; Wu et al., 2020; Shang et al., 2021; Deng & Hooi, 2021). Third, labeling information is often limited. Even if certain labels are present, in practice, many anomalies may still stay unidentified because labeling is laborious and expensive. Hence, unsupervised approaches are the most suitable choice. However, albeit many unsupervised detection methods were proposed (Ruff et al., 2018; Sabokrou et al., 2018; Malhotra et al., 2016; Hendrycks et al., 2019), they are not effective for multiple time series.
In this work, we explore the use of normalizing flows (Dinh et al., 2016; Papamakarios et al., 2017) for anomaly detection, based on a hypothesis that anomalies often lie on low density regions of the data distribution. Normalizing flows are a class of deep generative models for learning the underlying distribution of data samples. They are unsupervised and they resolve the label scarcity challenge aforementioned. An advantage of normalizing flows is that they are particularly effective in estimating the density of any sample. A recent work by Rasul et al. (2021) extends normalizing flows for time series data by expressing the density of a series through successive conditioning on historical data and applying conditional flows to learn each conditional density, paving ways to build sophisticated flows for multiple time series.
We address the high dimensionality and interdependency challenges by learning the relational structure among constituent series. To this end, Bayesian networks (Pearl, 1985; 2000) that model causal relationships of variables are a principled choice. A Bayesian network is a directed acyclic graph (DAG) where a node is conditionally independent of its non-descendents given its parents. Such a structure allows factorizing the intractable joint density of all graph nodes into a product of easyto-evaluate conditional densities of each node. Hence, learning the relational structure among constituent series amounts to identifying a DAG that maximizes the densities of observed data.
We propose a novel framework, GANF (Graph-Augmented Normalizing Flow), to augment a normalizing flow with graph structure learning and to apply it for anomaly detection. There are nontrivial technical problems to resolve to materialize this framework: (i) How does one inject a graph into a normalizing flow, which essentially maps one distribution to another? (ii) How does one learn a DAG, which is a discrete object, inside a continuous flow model? The solution we take is to factorize the density of a multiple time series along the attribute, the temporal, and the series dimensions and use a graph-based dependency encoder to model the conditional densities resulting from factorization. Therein, the graph adjacency matrix is a continuous variable and we impose a differentiable constraint to ensure that the corresponding graph is acyclic (Zheng et al., 2018; Yu et al., 2019). We propose a joint training algorithm to optimize both the graph adjacency matrix and the flow parameters.
In addition to resolving the high dimensionality and interdependency challenges, an advantage of modeling the relational graph structure among constituent series is that one can easily observe the dynamics of the data distribution from the graph. For time series datasets that span a long period, one naturally questions if the distribution changes over time. The graph structure is a useful indicator of distribution drift. We will study the graph evolution empirically observed.
We highlight the following contributions of this work:
• We propose a framework to augment a normalizing flow with graph structure learning, to model interdependencies exhibited inside multiple time series.
• We apply the augmented flow model to detect anomalies in multiple time series data and perform extensive empirical evaluation to demonstrate its effectiveness on real-world data sets.
• We study the evolution of the learned graph structure and identify distribution drift in time series data that span a long time period.

heading: 2 RELATED WORK
text: Anomaly Detection. Anomaly detection is a widely studied subject owing to its diverse applications. Recently, inspired by the success of deep learning, several deep anomaly detection methods are proposed and they achieve remarkable success on complex data, such as individual time series (Malhotra et al., 2016), images (Sabokrou et al., 2018), and videos (Ionescu et al., 2019). These
methods generally fall under three categories: deep one-class models, generative model-based methods, and transformation-based methods. Deep one-class models (Ruff et al., 2018; Wu et al., 2019) treat normal instances as the target class and identify instances that do not belong to this class. In generative model-based methods (Malhotra et al., 2016; Nguyen et al., 2019; Li et al., 2018a), an autoencoder or a generative adversarial network is used to model the data distribution. Then, an anomaly measure is defined, such as the reconstruction error in autoencoding. Transformationbased methods (Golan & El-Yaniv, 2018; Hendrycks et al., 2019) are based on the premise that transformations applied to normal instances can be identified while anomalies not. Various transformations such as rotations and affine transforms have been investigated. On the other hand, anomaly detection of multiple time series is under explored. Recently, Deng & Hooi (2021) study the use of graph neural networks in combination with structure learning to detect anomalies. Our method substantially differs from this work in that the learned structure is a Bayesian network, which allows density estimation. Moreover, the Bayesian network identifies conditional dependencies among the constituent series and induces a better interpretation of the graph as well as the data distribution.
Normalizing Flows. Normalizing flows are generative models that normalize complex real-world data distributions to “standard” distributions by using a sequence of invertible and differentiable transformations. Dinh et al. (2016) introduce a widely used normalizing flow architecture— RealNVP—for density estimation. Various extensions and improvements are proposed (Papamakarios et al., 2017; Hoogeboom et al., 2019; Kingma & Dhariwal, 2018). For example, Papamakarios et al. (2017) view an autoregressive model as a normalizing flow. To model temporal data, Rasul et al. (2021) use sequential models to parameterize conditional flows. Moreover, graph normalizing flows are proposed to handle graph structured data and improve predictions and generations (Liu et al., 2019). In contrast, normalizing flows for multiple time series are rarely studied in the literature. In this work, we develop a graph-augmented flow for density estimation and anomaly detection of multiple time series.

heading: 3 PRELIMINARIES
text: We first recall key concepts and familiarize the reader with notations to be used throughout the paper.

heading: 3.1 NORMALIZING FLOWS
text: Let x ∈ RD be a D-dimensional random variable. A normalizing flow is a vector-valued invertible mapping f(x) : RD → RD that normalizes the distribution of x to a “standard” distribution (or called base distribution). This distribution is usually taken to be an isotropic Gaussian or other ones that are easy to sample from and whose density is easy to evaluate. Let z = f(x) with probability density function q(z). With the change-of-variable formula, we can express the density of the x, p(x), by:
log p(x) = log q(f(x)) + log |det∇xf(x)|. (1)
In practical uses, the Jacobian determinant in (1) needs be easy to compute, so that the density p(x) can be evaluated. Moreover, as a generative model, the invertibility of f allows drawing new instances x = f−1(z) through sampling the base distribution. One example of such f is the masked autoregressive flow (Papamakarios et al., 2017), which yields z = [z1, . . . , zD] from x = [x1, . . . , xD] through zi = (xi − µi(x1:i−1)) exp(αi(x1:i−1)), (2) where µi and αi are neural networks such as the multilayer perceptron.
A flow may be augmented with conditional information h ∈ Rd with a possibly different dimension. Such a flow is a conditional flow and is denoted by f : RD × Rd → RD. The log-density of x conditioned on h admits the following formula:
log p(x|h) = log q(f(x;h)) + log |det∇xf(x;h)|. (3)
We now consider a normalizing flow for time series. Let X = [x1,x2, . . . ,xT ] denote a time series of length T , where xt ∈ RD. Through successive conditioning, the density of the time series can be written as: p(X) = p(x1)p(x2|x<2) · · · p(xT |x<T ), (4) where x<t denotes all variables before time t. When the conditional probabilities are parameterized, Rasul et al. (2021) propose to model each p(xt|x<t) as p(xt|ht−1), where ht−1 summarizes the
past information x<t. For example, ht−1 is the hidden state of a recurrent neural network before accepting input xt. Then, a conditional normalizing flow can be applied to evaluate each p(xt|ht−1).

heading: 3.2 BAYESIAN NETWORKS
text: Let Xi denote a general random variable, either scalar valued, vector valued, or even matrix valued. A Bayesian network of n variables (X1, . . . , Xn) is a directed acyclic graph of the variables as nodes. Let A denote the weighted adjacency matrix of the graph, where Aij 6= 0 if Xj is the parent of Xi. A Bayesian network describes the conditional independence among variables. Specifically, a node Xi is conditionally independent of its non-descendents given its parents. In other words, the density of the joint distribution of (X1, . . . , Xn) is
p(X1, . . . , Xn) = n∏ i=1 p(Xi|pa(Xi)), (5)
where pa(Xi) = {Xj : Aij 6= 0} denotes the set of parents of Xi.

heading: 4 PROBLEM STATEMENT
text: In this paper, we focus on unsupervised anomaly detection with multiple time series. The training set D consists of only unlabeled instances and we assume that the majority of them are not anomalies. Each instance X ∈ D contains n constituent series with D attributes and of length T ; i.e., X = (X1,X2, . . . ,Xn) where Xi ∈ RT×D. We use a Bayesian network (DAG) to model the relational structure of the constituent series Xi and augment a normalizing flow to compute the density of X through a factorization in the form (5). Let A ∈ Rn×n be the adjacency matrix of the DAG and let F : (X ,A) → Z denote the augmented flow. Because anomaly points tend to have low densities, we propose to conduct unsupervised anomaly detection by evaluating the density of a multiple time series computed through the augmented flow. The problem is formulated as the following.
Problem 1. Given a training set D = {Xi}|D|i=1 of multiple time series, we aim to simultaneously learn the adjacency matrix A of the Bayesian Network that represents the conditional dependencies among the constituent series, as well as the correspondingly graph-augmented normalizing flow F : (X ,A) → Z , which is used to estimate the density of an instance X . Here, Z is a random variable with a “simple” distribution, such as the anisotropic Gaussian.

heading: 5 METHOD
text: In this section, we materialize the graph-augmented normalizing flow F : (X ,A) → Z introduced in the problem statement and use it to compute the density of a multiple time series X . The central idea is factorization: we factorize p(X ) along the series dimension by using a Bayesian network and then factorize along the temporal dimension by using conditional normalizing flows. Then, we employ a novel graph-based dependency encoder to parameterize the conditional probabilities resulting from the factorization. The DAG used for factorization is a discrete object and is usually intractable to learn; however, the discrete structure is reflected in the dependency encoder through a graph adjacency matrix A that is differentiable. Moreover, the requirement that A must correspond to a DAG can be expressed as a differentiable equation. Hence, one can jointly optimize A and the flow components by using gradient based optimization. Once F is learned, the density p(X ) is straightforwardly evaluated for anomaly detection. An illustration of the framework GANF is shown in Figure 1.

heading: 5.1 FACTORIZATION
text: Figure 1 shows a toy example of a Bayesian network as a DAG. Based on (5), the density of a multiple time series X = (X1,X2, . . . ,Xn) can be computed as the product of p(Xi| pa(Xi)) for all nodes, where recall that pa(Xi) denotes the set of parents of Xi. Then, following Rasul et al. (2021), we further factorize each conditional density along the temporal dimension. Specifically, for a time step t, xit depends on its past history as well as its parents in the DAG. We write
p(X ) = n∏ i=1 p(Xi| pa(Xi)) = n∏ i=1 T∏ t=1 p(xit | pa(xi)1:t, xi1:t−1), (6)
where xi1:t−1 denotes the history of node i before time t, and pa(x i)1:t = {xj1:t : Aij 6= 0}. In the next subsection, we will parameterize each conditional density p(xit | pa(xi)1:t, xi1:t−1) by using a graph-based dependency encoder. Note that so far the factorization (6) has been based on the discrete structure of the Bayesian network. The dependency encoder we introduce next, however, uses the adjacency matrix A in a differentiable manner, which is sufficient to ensure that xit will not depend on nodes other than its parents and itself.

heading: 5.2 NEURAL NETWORK PARAMETERIZATION
text: According to Sec. 3.1, conditional densities p(xit | pa(xi)1:t, xi1:t−1) can be learned by using conditional normalizing flows. However, the conditional information pa(xi)1:t and xi1:t−1 cannot be directly used for parameterization, because its size is not fixed. Therefore, as is illustrated in Figure 1, we design a graph-based dependency encoder to summarize the conditional information into a fixed length vector dit ∈ Rd. Then, a conditional normalizing flow is used to evaluate p(xit|dit), which is equivalent to p(xit| pa(xi)1:t, xi1:t−1). Dependency Encoder. Since the history has an arbitrary length, we first employ a recurrent neural network (RNN) to map multiple time steps to a vector of fixed length. For a time series xi1:t, the recurrent model abstracts it into a hidden state hit ∈ Rd through the following recurrence
hit = RNN(x i t,h i t−1), (7)
where hit summarizes the time series up to step t. The RNN can be any sequential model, such as the LSTM (Hochreiter & Schmidhuber, 1997) and in a broad sense a transformer (Vaswani et al., 2017). We let the RNN parameters be shared across all nodes in the DAG to avoid overfitting and to reduce computational costs.
With (7), the conditional information of xit is all summarized in {h j t : Aij 6= 0} ∪ {hit−1}. Inspired by the success of GCN (Kipf & Welling, 2016) in node representation learning through neighborhood aggregation, we design a graph convolution layer to aggregate hidden states of the parents for dependency encoding. This layer produces dependency representations Dt = (d1t , . . . ,d n t ) for all constituent series at time t:
Dt = ReLU(AHtW1 +Ht−1W2) ·W3, (8) where Ht = (h1t , . . . ,h n t ) contains all the hidden states at time t. Here, W1 ∈ Rd×d and W2 ∈ Rd×d are parameters to transform the aggregated representations of the parents and the node’s historical information, respectively; while W3 ∈ Rd×d is an additional transformation to improve the dependency representation.
Density Estimation. With the dependency encoder, we obtain the representations dit of the conditional information. Then, a normalizing flow f : RD × Rd → RD conditioned on dit is applied to model each p(xit| pa(xi)1:t, xi1:t−1). Similar to the computation of the hidden states, the parameters of the conditional flow are also shared among nodes, to avoid overfitting. Based on (3), the conditional density of xit can be written as:
log p(xit | pa(xi)1:t, xi1:t−1) = log p(xit |dit) = log q(f(xit;dit)) + log |det∇xitf(x i t;d i t)|, (9)
where q(z) is chosen to be the standard normal N (z|0, I) with z ∈ RD. The conditional flow f can be any effective one proposed by the literature, such as RealNVP (Dinh et al., 2016) and
MAF (Papamakarios et al., 2017). Combining (9) and (6), we obtain the log-density of a multiple time series X :
log p(X ) = n∑ i=1 T∑ t=1 [ log q(f(xit;d i t)) + log |det∇xitf(x i t;d i t)| ] . (10)
Anomaly Measure. Because anomalies deviate significantly from the majority of the data instances, we hypothesize that their densities are low. Thus, we use the density computed by (10) as the anomaly measure, where a lower density indicates a more likely anomaly. Apart from evaluating the density for the entire X , the computation also produces conditional densities log p(Xi|pa(Xi)) =∑T t=1 log p(x i t|dit) for each constituent series Xi. We use this conditional density as the anomaly measure for constituent series. A low density p(X ) is caused by one or a few low conditional densities p(Xi|pa(Xi)) in the Bayesian network, suggesting that abnormal behaviors could be traced to individual series.

heading: 5.3 JOINT TRAINING
text: Learning a Bayesian network is a challenging combinatorial problem, due to the intractable search space superexponential in the number of nodes. A recent work by Zheng et al. (2018) proposes the equation tr(eA◦A) = n that characterizes the acyclicity of the corresponding graph of A, where e is matrix exponential and ◦ denotes element-wise multiplication. We will impose this equation as a constraint in the training of GANF.
Training Objective. Following the training of a usual normalizing flow, the joint density (likelihood) of the observed data is the training objective, which is equivalent to the Kullback–Leibler divergence between the true distribution of data and the flow recovered distribution. Together with the DAG constraint, the optimization problem reads
min A,θ L(A,θ) = 1 |D| |D|∑ i=1 − log p(Xi),
s.t. h(A) = tr(eA◦A)− n = 0,
(11)
where θ contains all neural network parameters, including those of the dependency encoder and the normalizing flow. Here, the DAG constraint h(A) admits an easy-to-evaluate gradient ∇h(A) = (eA◦A)T ◦ 2A, which allows a gradient based optimizer to solve (11). Training Algorithm. Problem (11) is a nonlinear equality-constrained optimization. Such problems are extensively studied and the augmented Lagrangian method (Bertsekas, 1999; Yu et al., 2019) is one of the most widely used approaches. The augmented Lagrangian is defined as
Lc = L(A,θ) + λh(A) + c
2 |h(A)|2, (12)
where λ and c denote the Lagrange multiplier and the penalty parameter, respectively. The general idea of the method is to gradually increase the penalty parameter to ensure that the constraint is eventually satisfied. Over iterations, λ as a dual variable will converge to the Lagrangian multiplier of (11). The update rule at the kth iteration reads the following:
Ak,θk = argmin A,θ Lck ; λk+1 = λk+ ckh(Ak); ck+1 = { ηck if |h(Ak)| > γ|h(Ak−1)| ; ck else, (13) where η ∈ (1,+∞) and γ ∈ (0, 1) are hyperparameters to be tuned. We set η and γ as 10 and 0.5, respectively. The subproblem of optimizing A and θ can be solved by using the Adam optimizer (Kingma & Ba, 2014). The training algorithm is summarized in Appendix A.

heading: 6 EXPERIMENTS
text: In this section, we conduct a comprehensive set of experiments to validate the effectiveness of the proposed GANF framework. In particular, they are designed to answer the following questions: • Q1: Can GANF accurately detect anomalies and estimate densities?
• Q2: Does the proposed graph structure learning help? Is the framework sufficiently flexible to include various normalizing flow backbones?
• Q3: What can one observe for a dataset spanning a long time? E.g., does the graph pattern change?

heading: 6.1 SETTINGS
text: Datasets. To evaluate the effectiveness of GANF for anomaly detection and density estimation, we conduct experiments on two power grid datasets, one water system dataset, and one traffic dataset.
• PMU-B and PMU-C: These two datasets correspond to two separate interconnects of the U.S. power grid, containing time series recorded by 38 and 132 phasor measurement units (PMUs), respectively. We process one-year data at the frequency of one second to form a ten-month training set, one-month validation set, and one-month test set. Each multiple time series is obtained by shifting a one-minute window. Additionally, to investigate distribution drift, we shift a one-month window to obtain multiple training/validation/test sets (12 in total, because of availability of twoyear data). Sparse grid events (anomalies) labeled by domain experts exist for evaluation; but note that the labels are both noisy and incomplete. These datasets are proprietary.
• SWaT: We also use a public dataset for evaluation. The Secure Water Treatment (SWaT) dataset originates from an operational water treatment test-bed coordinated with Singapore’s Public Utility Board (Goh et al., 2016). The data collects 51 sensor recordings lasting four days, at the frequency of one second. A total of 36 attacks were conducted, resulting in approximately 11% time steps as anomaly ground truths. We use a sliding window of 60 seconds to construct series data and perform a 60/20/20 chronological split for training, validation, and testing, respectively.
• METR-LA: This dataset is also public; it contains speed records of 207 sensors deployed on the highways of Los Angles, CA (Li et al., 2018b). No anomaly labels exist however and we use this dataset for exploratory analysis only. Results are deferred to Appendix E.
Evaluation metrics (under noisy labels). For SWaT, which offers reliable ground truths, we use the standard ROC and AUC metrics for evaluation. For the two PMU datasets, however, the resolution of the time series and the granularity of the events result in rather noisy ground truths. Hence, we adapt ROC for noisy labels. We smooth the time point of a “ground truth” event (anomaly) by introducing probabilities to the label. Specifically, the probability that a multiple time series starting at time t is a ground truth anomaly is maxi{exp(− (t−ti) 2
σ2 )}, where ti is the starting time of the ith labeled anomaly. Then, when computing the confusion matrix, we sum probabilities rather than counting 0/1s. The smoothing window σ is chosen to be 6 time steps.
Baselines. We compare with the following representative, state-of-the-art deep methods.
• EncDecAD (Malhotra et al., 2016): In this method, an autoencoder based on LSTM is trained. The reconstruction error is used as the anomaly measure.
• DeepSVDD (Ruff et al., 2018): This method minimizes the volume of a hypersphere that encloses the representations of data. Samples distant from the hypersphere center are considered anomalies.
• ALOCC (Sabokrou et al., 2020): In this GAN-based method, the generator learns to reconstruct normal instances, while the discriminator works as an anomaly detector.
• DROCC (Goyal et al., 2020): This method performs adversarial training to learn robust representations of data and identifies anomalies.
• DeepSAD (Ruff et al., 2020): This method extends DeepSVDD with a semi-supervised loss term for training. We use noisy labels as supervision.
To apply these baselines on multiple time series, we concatenate the constituent series along the attribute dimension (resulting in high-dimensional series) and use LSTM or CNN as the backbones. On the other hand, for the proposed method, we use LSTM as the RNN model and MAF as the normalizing flow. See Appendix C for more implementation details.

heading: 6.2 PERFORMANCE OF ANOMALY DETECTION AND DENSITY ESTIMATION
text: To answer Q1, we evaluate quantitatively and qualitatively on datasets with labels.
Anomaly detection. We compare GANF with the aforementioned baselines in Table 1, where standard deviations of the AUC scores are additionally reported through five random repetitions of model training. The table suggests an overwhelmingly high AUC score achieved by GANF. Observations follow. (i) GANF outperforms generative model-based methods (EncDecAD and ALOCC). Being a generative model as well, normalizing flows augmented with a graph structure leverage the interdependencies of constituent series more effectively, leading to a substantial improvement in detection. (ii) GANF significantly outperforms deep one-class models (DeepSVDD and DROCC), corroborating the appeal of using densities for detection. (iii) GANF also performs better than the semi-supervised method DeepSAD, probably because such methods rely on high quality labels for supervision (especially in the case of label scarcity) and they are less effective facing noisy labels.
Besides a single score, we also plot the ROC curve in Figure 2. One sees that the curve of GANF dominates those of others. This behavior is generally more salient in the low false-alarm regime.
Density estimation. We investigate the densities estimated by GANF, shown in Figure 3. Distributions of the log-densities in the test set are shown in Figure 3a. We use log-density as the anomaly measure; the lower the more likely. Note that the vertical axis is in the log-scale. One sees that a log-density of 16 approximately separates the majority normal instances from the minority anomalies. To cross-verify that the instances with low densities are suspiciously anomalous, we investigate Figure 3b, which is a temporal plot of log-densities for a week, overlaid with given labels. From this plot, one sees that the noisily labeled series generally have low densities or are near a low density time step. Additionally, GANF discovers a few suspicious time steps with low densities undetected earlier. These new discoveries raise interest to power system experts for analysis and archiving.

heading: 6.3 ABLATION STUDY
text: To answer Q2, we conduct an ablation study (including varying architecture components) to investigate impacts of DAG structure learning and the flexibility of the GANF framework. To investigate the power of modeling pairwise relationship, we train a variant GANF\G that factorizes p(X ) = ∏n i=1 p(X
i); i.e., assuming independence among constituent series. To investigate the effectiveness of graph structure learning, we train a variant GANF\D that decomposes the joint density as p(X ) = ∏n i=1 p(X
i|X<i); i.e., a full decomposition without a DAG. It is equivalent to concatenating the series along the attribute dimension and running MAF on the resulting series. To verify the contribution of joint training of A and θ, we train a variant GANF\T where A is separately learned by using NOTEARS (Zheng et al., 2018). To prove the flexibility of GANF, we replace the MAF-based normalizing flow by RealNVP, denoted as GANFRNVP.
Results are presented in Table 2. Apart from AUC-ROC, the log-density is also reported. Observations follow. (i) GANF significantly outperforms GANF\G and GANF\D, corroborating the importance of interdependency modeling among constituent series. Note that GANF\D results in particularly poor performance in general, likely because the high dimensional input (resulting from concatenating too many series) impedes the learning of normalizing flows. (iii) GANF\T is slightly better than GANF\G, because of the presence of relational modeling, but it cannot match the performance of GANFRNVP and GANF that jointly train the DAG and the flow. (ii) These latter two models are the best for both datasets and both metrics. MAF works more often better than RealNVP.

heading: 6.4 EVOLUTION OF THE DAG STRUCTURE
text: To answer Q3, we investigate how the learned DAG evolves by shifting the train/validation/test sets month by month. The graphs within the first three-month shiftings are shown in Figure 4 and more can be found in Appendix F. In addition to the graph structure, we plot in Figure 5 the learned edge weights over time, one column per edge. The appearance and disappearance of edges demonstrate changes of the conditional independence structure among constituent series over time, suggesting data distribution drift (i.e., a change of internal data generation mechanism). It is interesting to observe the seasonal effect. The columns (edges) in Figure 5 can be loosely grouped in three clusters: those persisting the entire year, those appearing in the first half of the year, and those existing more briefly (e.g., within a season). Such a pattern plausibly correlates with electricity consumption, which is also seasonal. Were spatial information of the PMUs known, these identified DAGs would help mapping the seasonal patterns to geography and help planning a more resilient grid.

heading: 7 CONCLUSIONS
text: In this paper, we present a graph-augmented normalizing flow GANF for anomaly detection of multiple time series. The graph is materialized as a Bayesian network, which models the conditional dependencies among constituent time series. A graph-based dependency decoder is designed to summarize the conditional information needed by the normalizing flow that calculates series density. Anomalies are detected through identifying instances with low density. Extensive experiments on real-world datasets demonstrate the effectiveness of the framework. Ablation studies confirm the contribution of the learned graph structure in anomaly detection. Additionally, we investigate the evolution of the graph and offer insights of distribution drift over time.
ACKNOWLEDGMENT AND DISCLAIMER This material is based upon work supported by the Department of Energy under Award Number(s) DE-OE0000910. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

heading: A TRAINING ALGORITHM
text: We summarize the training method in Algorithm 1.
Algorithm 1 Training Algorithm of GANF Input: Training set D, hyperparameters η and γ Output: GANF F and adjacency matrix A of the DAG
1: Initialize c← 0 and initialize λ randomly 2: for k = 0, 1, 2, . . . do 3: Compute Ak and θk as a minimizer of (12) by using the Adam optimizer, where the loss
L and the constraint h are defined in (11), the log-density log p(X ) is defined in (10), the dependency representation dit is defined in (8), the hidden state h i t is defined in (7), and the
conditional flow f is RealNVP or MAF 4: Update Lagrange multiplier λ← λ+ ch(Ak) 5: if k > 0 and |h(Ak)| > γ|h(Ak−1)| then 6: c← ηc 7: end if 8: if h(Ak) == 0 then 9: break
10: end if 11: end for 12: return A and F (including f , the neural network (8), and the RNN (7))

heading: B CODE
text: Code is available at https://github.com/EnyanDai/GANF.

heading: C ADDITIONAL DETAILS OF EXPERIMENT SETTINGS
text: C.1 IMPLEMENTATION DETAILS OF GANF. An LSTM is used as the RNN model in the dependency encoder. For normalizing flows, we use MAF with six flow blocks. All hidden dimensions as set as 32. The initial learning rate is set as 0.001 for the adjacency matrix A and the model parameters θ. Learning rate decay is 0.1. To avoid gradient explosion, we clip the gradients whose values are larger than 1.0.
For hyperparameter tuning, we select the hyperparameters that yield the highest log-density on the validation set. Specifically, we conduct grid search by varying the number of normalizing flow blocks from {1, 2, 4, 6, 8}, the learning rate from {0.003, 0.001, 0.0003, 0.0001}, and the hidden dimension from {16, 32, 64, 128}.
C.2 IMPLEMENTATION DETAILS OF BASELINES.
• EncDecAD (Malhotra et al., 2016): This method is applied for anomaly detection on time series. Thus, we concatenate constituent series along the attribute dimension and adopt the code released by the authors in https://github.com/chickenbestlover/ RNN-Time-series-Anomaly-Detection.
• DeepSVDD (Ruff et al., 2018): To handle time series data, we replace the backbone to an LSTM based on the official implementation https://github.com/lukasruff/ Deep-SVDD-PyTorch.
• ALOCC (Sabokrou et al., 2020): We use the official implementation released by the authors in https://github.com/khalooei/ALOCC-CVPR2018. We replace the two-dimensional convolution to one-dimensional convolution to build a GAN for time series data.
• DROCC (Goyal et al., 2020): Similar to other baselines, this method is proposed for tabular data and image data. We replace the backbone to LSTM to deal with multiple time series by revising the encoder in https://github.com/microsoft/EdgeML/tree/master/pytorch.
• DeepSAD (Ruff et al., 2020): This is a semi-supervised approach, which requires labeling. We utilize the noisy labels in PMU-B and PMU-C as supervision. We use LSTM as the
backbone, based on the official implementation in https://github.com/lukasruff/ Deep-SAD-PyTorch.
All hyperparameters of the baselines are tuned based on the validation set to make fair comparisons.

heading: D TIME COMPLEXITY ANALYSIS
text: The GANF framework involves Bayesian network structure learning, which is known to be highly challenging, owing to the intractable search space superexponential in the number of graph nodes. In this work, we formulate a continuous optimization of the graph structure, so that the training of GANF is more scalable. In what follows, we analyze the time complexity.
Recall that each instance X of the multiple time series dataset contains n constituent series with D attributes and of length T ; i.e., X = (X1,X2, . . . ,Xn) where Xi ∈ RT×D. In the evaluation of the model, the dominant costs appear in running the dependency encoder and the normalizing flow. For the dependency encoder, an RNN is first deployed to map the multiple time series to hidden vectors; the time complexity is O(nTD). Then, graph convolution is conducted to obtain dependency vectors; the convolution cost is O(n2T ). For the normalizing flow module, the time complexity is O(nTD). Therefore, the time complexity of computing log-density of one instance is O(nT (D + n)). If we use a batch size B for training, the time cost of calculating the augmented Lagrangian L(A,θ) in (12) is O(nBT (D + n)). Additionally, the time cost of calculating the constraint h(A) isO(n3). Thus, the overall time complexity of one training iteration isO(n(BTD+ BTn+ n2)).

heading: E RESULTS FOR METR-LA
text: METR-LA contains speed records of 207 sensors deployed on the highways of Los Angles, CA (Li et al., 2018b). The records are in four months at the frequency of five minutes. We shift a one-hour window to obtain multiple time series. The first three months are used for training and the last month is split in halves for validation and testing. No anomaly labels exist however and we use this dataset for exploratory analysis only.
Figure 6a shows the traffic speed on four main highways on June 13, 2012. Each speed is the average over all sensors on the same highway. We observe that despite spatial proximity, the speeds vary significantly around 4PM (rush hour) but they are unanimously high around 5AM and 8PM– 12AM. The estimated densities, shown in Figure 6b, tracks this pattern rather closely, with rush hours corresponding to low density and night traffics corresponding to high density. Note the nature of traffic: speed varies smoothly on the macroscopic level and hence does density, too. Such a phenomenon is in striking contrast to power systems where events are rare and abrupt.

heading: F ADDITIONAL ANOMALY DETECTION RESULTS ON PMU DATASETS
text: See Figure 7 and Figure 8 for additional anomaly detection results on the test sets of PMU-C and PMU-B, respectively. The observations are rather similar to those of Figure 3b.

heading: G ADDITIONAL RESULTS FOR DAG EVOLUTION

Review summaries:
The authors propose a model for anomaly detection on multivariate, high-dimensional time series data, where there are statistical dependencies across the different time series (referred to as "constituents" in the paper). The authors model statistical dependence using a Bayesian network and temporal dependence using an RNN. One key feature of the proposed model is that it simultaneously learns the structure of the Bayesian network AND a representation of the temporal state. The authors combine these two representations, which then gets used as input to a normalized flow technique, which makes computation of the likelihood of the data mathematically and computationally feasible.  The authors show that the proposed model has improved performance for the anomaly detection task compared to recently-proposed neural-network anomaly detection methods.

The paper addresses the task of detecting anomalies in multiple time-series, for the setting where multiple instances of a fixed time-window are available for learning. The authors propose a method that involves using a directed acyclic graph (DAG) to model the (contemporaneous) dependencies between the variables in the time-series. An RNN is applied to form a representative state for each time-series at each point in time. A mapping is applied using the DAG to process the states from the previous time-step and the ancestor states from the current time-step to construct a “dependency representation” vector for each time series. A normalizing flow is then applied to map to a base distribution so that the log-density of the instance (the observed multiple time-series) can be evaluated. The conditional densities for each time series can also be evaluated. The paper includes an evaluation of the performance of the approach on a proprietary dataset, with a comparison to state-of-the-art baselines. Some qualitative results are also included for the task of density estimation for a public traffic dataset. 

This paper proposes to combine learning of the DAG structure between variables and learning of distribution using normalizing flows for unsupervised anomaly detection from multivariate time series. Once the distribution is learned, one can detect anomalies as low density regions. The effectiveness of the proposed method is empirically evaluated on real-world datasets.

This paper focuses on unsupervised anomaly detection in multivariate time series. This is a important problem because anomalies are rare and labeling anomalies for supervised learning is labor intensive. To make progress in this challenging domain, they  use learned graph structure and normalizing flows, which seem well suited for the task. They evaluate their methods on multiple datasets, and compare against multiple modern deep networks.