# GRAPH-AUGMENTED NORMALIZING FLOWS FOR ANOMALY DETECTION OF MULTIPLE TIME SERIES


**Enyan Dai[∗]**
Pennsylvania State University
emd5759@psu.edu


**Jie Chen[†]**
MIT-IBM Watson AI Lab, IBM Research
chenjie@us.ibm.com

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.

1 INTRODUCTION

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 (Sch¨olkopf
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.

2 RELATED WORK

**Anomaly Detection. Anomaly detection is a widely studied subject owing to its diverse applica-**
tions. 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.

3 PRELIMINARIES

We first recall key concepts and familiarize the reader with notations to be used throughout the paper.

3.1 NORMALIZING FLOWS
Let x ∈ R[D] be a D-dimensional random variable. A normalizing flow is a vector-valued invertible
mapping f (x) : R[D] _→_ R[D] 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 au-_
toregressive 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 ∈ R[d] with a possibly different dimension.
Such a flow is a conditional flow and is denoted by f : R[D] _× R[d]_ _→_ R[D]. 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
written as:of length T, where xt ∈ R[D]. Through successive conditioning, the density of the time series can be
_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).

3.2 BAYESIAN NETWORKS

Let X _[i]_ denote a general random variable, either scalar valued, vector valued, or even matrix valued.
A Bayesian network of n variables (X [1], . . ., X _[n]) is a directed acyclic graph of the variables as_
nodes. Let A denote the weighted adjacency matrix of the graph, where Aij = 0 if X _[j]_ is the parent
_̸_
of X _[i]. A Bayesian network describes the conditional independence among variables. Specifically,_
a node X _[i]_ is conditionally independent of its non-descendents given its parents. In other words, the
density of the joint distribution of (X [1], . . ., X _[n]) is_


_p(X_ [1], . . ., X _[n]) =_


_p(X_ _[i]| pa(X_ _[i])),_ (5)
_i=1_

Y


where pa(X _[i]) =_ _X_ _[j]_ : Aij = 0 denotes the set of parents of X _[i]._
_{_ _̸_ _}_

4 PROBLEM STATEMENT

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 =
(X[1], X[2], . . ., X[n]) where X[i] _∈_ R[T][ ×][D]. We use a Bayesian network (DAG) to model the relational
structure of the constituent series X[i] and augment a normalizing flow to compute the density of X
through a factorization in the form (5). Let A ∈ R[n][×][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}i[|D|]=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._

5 METHOD

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.

5.1 FACTORIZATION

Figure 1 shows a toy example of a Bayesian network as a DAG. Based on (5), the density of a
multiple time series X = (X[1], X[2], . . ., X[n]) can be computed as the product of p(X[i]| pa(X[i])) for
all nodes, where recall that pa(X[i]) denotes the set of parents of X[i]. Then, following Rasul et al.
(2021), we further factorize each conditional density along the temporal dimension. Specifically, for
a time step t, x[i]t [depends on its past history as well as its parents in the DAG. We write]


_p(X[i]| pa(X[i])) =_
_i=1_ _i=1_

Y Y


_p(x[i]t_ 1:t 1[)][,] (6)
_t=1_ _[|][ pa(][x][i][)][1:][t][,][ x][i]_ _−_

Y


_p(X_ ) =


-----

Figure 1: Illustration of a Bayesian network and the proposed framework GANF.

where x[i]1:t−1 [denotes the history of node][ i][ before time][ t][, and][ pa(][x][i][)][1:][t] [=][ {][x][j]1:t [:][ A][ij][ ̸][= 0][}][. In the]
next subsection, we will parameterize each conditional density p(x[i]t 1:t 1[)][ by using]
a graph-based dependency encoder. Note that so far the factorization (6) has been based on the[|][ pa(][x][i][)][1:][t][,][ x][i] _−_
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 x[i]t [will not]
depend on nodes other than its parents and itself.

5.2 NEURAL NETWORK PARAMETERIZATION
According to Sec. 3.1, conditional densities p(x[i]t 1:t 1[)][ can be learned by using con-]
ditional normalizing flows. However, the conditional information[|][ pa(][x][i][)][1:][t][,][ x][i] _− pa(x[i])1:t and x[i]1: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
which is equivalent toa fixed length vector d p[i]t(x[∈][i]t[|][R][ pa(][d][. Then, a conditional normalizing flow is used to evaluate][x][i][)][1:][t][,][ x][i]1:t 1[)][.] _[ p][(][x]t[i][|][d][i]t[)][,]_
_−_

**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 x[i]1:t[, the]
recurrent model abstracts it into a hidden state h[i]t

_[∈]_ [R][d][ through the following recurrence]

**h[i]t** [=][ RNN][(][x]t[i][,][ h]t[i] 1[)][,] (7)
_−_

where h[i]t [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 x[i]t [is all summarized in][ {][h]t[j] [:][ A][ij] _t_ 1[}][. Inspired]
by the success of GCN (Kipf & Welling, 2016) in node representation learning through neighbor-[̸][= 0][} ∪{][h][i]−
hood aggregation, we design a graph convolution layer to aggregate hidden states of the parents for
dependency encoding. This layer produces dependency representations Dt = (d[1]t _[, . . .,][ d][n]t_ [)][ for all]
constituent series at time t:

**Dt = ReLU(AHtW1 + Ht** 1W2) **W3,** (8)
_−_ _·_

whereWnode’s historical information, respectively; whileimprove the dependency representation.2 ∈ HRt = ([d][×][d] are parameters to transform the aggregated representations of the parents and theh[1]t _[, . . .,][ h][n]t_ [)][ contains all the hidden states at time] W3 ∈ R[d][×][d] is an additional transformation to[ t][.] Here, W1 ∈ R[d][×][d] and

**Density Estimation. With the dependency encoder, we obtain the representations d[i]t** [of the condi-]
tional information. Then, a normalizing flow f : R[D] _× R[d]_ _→_ R[D] conditioned on d[i]t [is applied to]
model each p(x[i]t[|][ pa(][x][i][)][1:][t][,][ x][i]1:t 1[)][. Similar to the computation of the hidden states, the param-]
_−_
eters of the conditional flow are also shared among nodes, to avoid overfitting. Based on (3), the
conditional density of x[i]t [can be written as:]

log p(x[i]t 1:t 1[) = log][ p][(][x][i]t _t[) = log][ q][(][f]_ [(][x][i]t[;][ d][i]t[)) + log][ |][ det][ ∇]x[i]t **[f]** [(][x]t[i][;][ d][i]t[)][|][,] (9)

_[|][ pa(][x][i][)][1:][t][,][ x][i]_ _−_ _[|][ d][i]_

where q(z) is chosen to be the standard normal N (z|0, I) with z ∈ R[D]. 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 : _n_ _T_


log q(f (x[i]t[;][ d][i]t[)) + log][ |][ det][ ∇]x[i]t **[f]** [(][x]t[i][;][ d][i]t[)][|] _._ (10)
i


log p(X ) =


_i=1_


_t=1_


**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 entireT _X_, the computation also produces conditional densities log p(X[i]| pa(X[i])) =
_t=1_ [log][ p][(][x]t[i][|][d][i]t[)][ for each constituent series][ X][i][. We use this conditional density as the anomaly]
tiesmeasure for constituent series. A low densityP _p(X[i]| pa(X[i])) in the Bayesian network, suggesting that abnormal behaviors could be traced to p(X_ ) is caused by one or a few low conditional densiindividual series.

5.3 JOINT TRAINING
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(e[A][◦][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 (likeli-**
hood) 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

1 _|D|_
min (A, θ) = log p( _i),_
**A,θ** _L_ _−_ _X_ (11)

_|D|_ _i=1_

X

s.t. _h(A) = tr(e[A][◦][A]) −_ _n = 0,_

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) =
(e[A][◦][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

_c =_ (A, θ) + λh(A) + _[c]_ (12)
_L_ _L_ 2 _[|][h][(][A][)][|][2][,]_

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:

_ηc[k]_ if _h(A[k])_ _> γ_ _h(A[k][−][1])_ ;
**A[k], θ[k]** = arg min _λ[k][+1]_ = λ[k] + c[k]h(A[k]); _c[k][+1]_ = _|_ _|_ _|_ _|_
**A,θ** _c[k]_ else,

_[L][c][k]_ [;] 
(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.

6 EXPERIMENTS

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?


-----

6.1 SETTINGS

**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( _σ[2]_ ), where ti is the starting time of the
_{_ _−_ [(][t][−][t][i][)][2] _}_

_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.

6.2 PERFORMANCE OF ANOMALY DETECTION AND DENSITY ESTIMATION

Table 1: AUC-ROC (%) of anomaly detection.

Dataset EncDecAD DeepSVDD ALOCC DROCC DeepSAD GANF

PMU-B 55.6±1.8 55.6±3.3 62.9±2.2 58.6±3.0 63.7±0.9 **67.5±0.8**

PMU-C 53.7±0.5 56.9±0.9 60.9±1.3 61.9±2.7 60.1±1.4 **70.6±3.3**

SWaT 76.5±0.7 68.8±2.0 75.4±2.3 73.3±1.6 75.4±1.2 **79.6±0.9**


-----

1.0

0.8

0.6

0.0

0.0 0.2 0.4 0.6 0.8 1.0

GANF
DeepSAD
DROCC
ALOCC
DeepSVDD
EncDecAD

False Positive Rate


1.0

0.8

0.6

0.0

|Col1|Col2|Col3|
|---|---|---|
||GANF DeepSAD DROCC ALOCC DeepSVDD EncDecAD||


0.0 0.2 0.4 0.6 0.8 1.0

GANF
DeepSAD
DROCC
ALOCC
DeepSVDD
EncDecAD

False Positive Rate


1.0

0.8

0.6

0.0

0.0 0.2 0.4 0.6 0.8 1.0

GANF
DeepSAD
DROCC
ALOCC
DeepSVDD
EncDecAD

False Positive Rate


(a) PMU-B (b) PMU-C (c) SWaT

Figure 2: ROC curves of anomaly detection on various datasets.


18

17

16

15


10[4]

10[3]

10[2]

10[1]

10[0]


14 15 16 17 18

Anormaly Measure = log(density)


11-03 11-05 11-07 11-09

|Col1|L|abeled|Anoma|ly|Col6|Score|of T|im|estam|p|
|---|---|---|---|---|---|---|---|---|---|---|
||||||||||||


Labeled Anomaly Score of Timestamp

Time


(a)


(a)

(b)

Figure 3: Qualitative evaluation of GANF on PMU-C. (a) Distribution of log-densities on the test
set (note in log scale). (b) Anomaly detection results for a week in the test set.


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. Distribu-**
tions 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.


6.3 ABLATION STUDY

Table 2: Performance of variants of the proposed method.


Dataset Metrics GANF\G GANF\D GANF\T GANFRNVP GANF

AUC-ROC 0.641 0.643 0.653 0.661 **0.678**
PMU-B
Log-Density 15.31 8.70 15.09 15.90 **16.22**

AUC-ROC 0.630 0.544 0.688 0.703 **0.705**
PMU-C
Log-Density 15.55 8.94 15.70 **17.06** 16.98


-----

(a) Initial graph (b) One-month shift (c) Two-month shift (d) Three-month shift

Figure 4: Evolution of the learned DAG on PMU-B over time.

Figure 5: Evolution of edge weights in the DAG learned by GANF over time (PMU-B).

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_ ) = _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[Q] p[n](X ) = _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

[Q][n] _\_
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.

6.4 EVOLUTION OF THE DAG STRUCTURE
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.

7 CONCLUSIONS

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.

REFERENCES

Terje Aven. Risk assessment and risk management: Review of recent advances on their foundation.
_European Journal of Operational Research, 253(1):1–13, 2016._

Dimitri P. Bertsekas. Nonlinear Programming. Athena Scientific, 2nd edition, 1999.

Ailin Deng and Bryan Hooi. Graph neural network-based anomaly detection in multivariate time
series. In AAAI, 2021.

Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real nvp. arXiv
_preprint arXiv:1605.08803, 2016._

Jonathan Goh, Sridhar Adepu, Khurum Nazir Junejo, and Aditya Mathur. A dataset to support
research in the design of secure water treatment systems. In International conference on critical
_information infrastructures security, pp. 88–99. Springer, 2016._

Izhak Golan and Ran El-Yaniv. Deep anomaly detection using geometric transformations. arXiv
_preprint arXiv:1805.10917, 2018._

Sachin Goyal, Aditi Raghunathan, Moksh Jain, Harsha Vardhan Simhadri, and Prateek Jain. Drocc:
Deep robust one-class classification. In International Conference on Machine Learning, pp.
3711–3721. PMLR, 2020.

Kathrin Grosse, Praveen Manoharan, Nicolas Papernot, Michael Backes, and Patrick McDaniel. On
the (statistical) detection of adversarial examples. arXiv preprint arXiv:1702.06280, 2017.

Dan Hendrycks, Mantas Mazeika, Saurav Kadavath, and Dawn Song. Using self-supervised learning
can improve model robustness and uncertainty. arXiv preprint arXiv:1906.12340, 2019.

Sepp Hochreiter and J¨urgen Schmidhuber. Long short-term memory. Neural computation, 9(8):
1735–1780, 1997.

Emiel Hoogeboom, Rianne Van Den Berg, and Max Welling. Emerging convolutions for generative
normalizing flows. In International Conference on Machine Learning, pp. 2771–2780. PMLR,
2019.

Radu Tudor Ionescu, Fahad Shahbaz Khan, Mariana-Iuliana Georgescu, and Ling Shao. Objectcentric auto-encoders and dummy anomalies for abnormal event detection in video. In Proceed_ings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 7842–7851,_
2019.

JooSeuk Kim and Clayton D Scott. Robust kernel density estimation. The Journal of Machine
_Learning Research, 13(1):2529–2565, 2012._

Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint
_arXiv:1412.6980, 2014._

Diederik P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1x1 convolutions.
_arXiv preprint arXiv:1807.03039, 2018._


-----

Thomas Kipf, Ethan Fetaya, Kuan-Chieh Wang, Max Welling, and Richard Zemel. Neural relational
inference for interacting systems. In ICML, 2018.

Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. arXiv preprint arXiv:1609.02907, 2016.

B Ravi Kiran, Dilip Mathew Thomas, and Ranjith Parakkal. An overview of deep learning based
methods for unsupervised and semi-supervised anomaly detection in videos. Journal of Imaging,
4(2):36, 2018.

Dan Li, Dacheng Chen, Jonathan Goh, and See-kiong Ng. Anomaly detection with generative
adversarial networks for multivariate time series. arXiv preprint arXiv:1809.04758, 2018a.

Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. Diffusion convolutional recurrent neural network: Data-driven traffic forecasting. In ICLR, 2018b.

Jenny Liu, Aviral Kumar, Jimmy Ba, Jamie Kiros, and Kevin Swersky. Graph normalizing flows.
_arXiv preprint arXiv:1905.13177, 2019._

Pankaj Malhotra, Anusha Ramakrishnan, Gaurangi Anand, Lovekesh Vig, Puneet Agarwal, and
Gautam Shroff. Lstm-based encoder-decoder for multi-sensor anomaly detection. arXiv preprint
_arXiv:1607.00148, 2016._

Duc Tam Nguyen, Zhongyu Lou, Michael Klar, and Thomas Brox. Anomaly detection with
multiple-hypotheses predictions. In International Conference on Machine Learning, pp. 4800–
4809. PMLR, 2019.

George Papamakarios, Theo Pavlakou, and Iain Murray. Masked autoregressive flow for density
estimation. arXiv preprint arXiv:1705.07057, 2017.

Emanuel Parzen. On estimation of a probability density function and mode. The annals of mathe_matical statistics, 33(3):1065–1076, 1962._

Judea Pearl. Bayesian networks: A model of self-activated memory for evidential reasoning. In
_Proceedings of the 7th Conference of the Cognitive Science Society, 1985._

Judea Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, 2000.

Marco AF Pimentel, David A Clifton, Lei Clifton, and Lionel Tarassenko. A review of novelty
detection. Signal Processing, 99:215–249, 2014.

Kashif Rasul, Abdul-Saboor Sheikh, Ingmar Schuster, Urs Bergmann, and Roland Vollgraf. Multivariate probabilistic time series forecasting via conditioned normalizing flows. In ICLR, 2021.

Riya Roy and K Thomas George. Detecting insurance claims fraud using machine learning techniques. In 2017 International Conference on Circuit, Power and Computing Technologies (IC_CPCT), pp. 1–6. IEEE, 2017._

Lukas Ruff, Robert Vandermeulen, Nico Goernitz, Lucas Deecke, Shoaib Ahmed Siddiqui, Alexander Binder, Emmanuel M¨uller, and Marius Kloft. Deep one-class classification. In International
_conference on machine learning, pp. 4393–4402. PMLR, 2018._

Lukas Ruff, Robert A. Vandermeulen, Nico G¨ornitz, Alexander Binder, Emmanuel M¨uller, KlausRobert M¨uller, and Marius Kloft. Deep semi-supervised anomaly detection. In International
_[Conference on Learning Representations, 2020. URL https://openreview.net/forum?](https://openreview.net/forum?id=HkgH0TEYwH)_
[id=HkgH0TEYwH.](https://openreview.net/forum?id=HkgH0TEYwH)

Mohammad Sabokrou, Mohammad Khalooei, Mahmood Fathy, and Ehsan Adeli. Adversarially
learned one-class classifier for novelty detection. In Proceedings of the IEEE Conference on
_Computer Vision and Pattern Recognition, pp. 3379–3388, 2018._

Mohammad Sabokrou, Mahmood Fathy, Guoying Zhao, and Ehsan Adeli. Deep end-to-end oneclass classifier. _IEEE transactions on neural networks and learning systems, 32(2):675–684,_
2020.


-----

Bernhard Sch¨olkopf, John C Platt, John Shawe-Taylor, Alex J Smola, and Robert C Williamson.
Estimating the support of a high-dimensional distribution. Neural computation, 13(7):1443–1471,
2001.

Youngjoo Seo, Micha¨el Defferrard, Pierre Vandergheynst, and Xavier Bresson. Structured sequence
modeling with graph convolutional recurrent networks. arXiv:1612.07659, 2016.

Chao Shang, Jie Chen, and Jinbo Bi. Discrete graph structure learning for forecasting multiple time
series. In ICLR, 2021.

Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez,
Łukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in neural information
_processing systems, pp. 5998–6008, 2017._

Peng Wu, Jing Liu, and Fang Shen. A deep one-class neural network for anomalous event detection
in complex scenes. IEEE transactions on neural networks and learning systems, 31(7):2609–
2622, 2019.

Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chengqi Zhang. Connecting the dots: Multivariate time series forecasting with graph neural networks. In KDD, 2020.

Bing Yu, Haoteng Yin, and Zhanxing Zhu. Spatio-temporal graph convolutional networks: A deep
learning framework for traffic forecasting. In IJCAI, 2018.

Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural
networks. In ICML, 2019.

Ling Zhao, Yujiao Song, Chao Zhang, Yu Liu, Pu Wang, Tao Lin, Min Deng, and Haifeng Li.
T-GCN: A temporal graph convolutional network for traffic prediction. IEEE Transactions on
_Intelligent Transportation Systems, 2019._

Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. DAGs with NO TEARS: Continuous optimization for structure learning. In NeurIPS, 2018.


-----

A TRAINING ALGORITHM

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 A[k] 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 d[i]t [is defined in (8), the hidden state][ h]t[i] [is defined in (7), and the]
conditional flow f is RealNVP or MAF


4: Update Lagrange multiplier λ ← _λ + ch(A[k])_

5: **if k > 0 and |h(A[k])| > γ|h(A[k][−][1])| then**

6: _c ←_ _ηc_

7: **end if**

8: **if h(A[k]) == 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))

B CODE

[Code is available at https://github.com/EnyanDai/GANF.](https://github.com/EnyanDai/GANF)

C ADDITIONAL DETAILS OF EXPERIMENT SETTINGS

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/](https://github.com/chickenbestlover/RNN-Time-series-Anomaly-Detection)
[RNN-Time-series-Anomaly-Detection.](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/](https://github.com/lukasruff/Deep-SVDD-PyTorch)
[Deep-SVDD-PyTorch.](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](https://github.com/khalooei/ALOCC-CVPR2018)
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.](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/](https://github.com/lukasruff/Deep-SAD-PyTorch)
[Deep-SAD-PyTorch.](https://github.com/lukasruff/Deep-SAD-PyTorch)

All hyperparameters of the baselines are tuned based on the validation set to make fair comparisons.

D TIME COMPLEXITY ANALYSIS

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 = (X[1], X[2], . . ., X[n]) where X[i] _∈_ R[T][ ×][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(n[2]T ). 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) is O(n[3]). Thus, the overall time complexity of one training iteration is O(n(BTD+
_BTn + n[2]))._

E RESULTS FOR METR-LA

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.

(a) Average speed of the main roads. (b) Estimated log-density.

Figure 6: Density estimation for METR-LA.

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.


-----

ADDITIONAL ANOMALY DETECTION RESULTS ON PMU DATASETS


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.

17

16

Anormaly Measure


15

2017-11-08 2017-11-10 2017-11-12 2017-11-14

|Col1|Col2|Col3|Labe|Col5|Col6|Col7|Col8|Col9|Col10|Col11|Col12|
|---|---|---|---|---|---|---|---|---|---|---|---|
|||||Labe|led Anomaly S|cor|e of Times|tamp||||
|||||||||||||


Time


17

16


2017-11-15 2017-11-17 2017-11-19 2017-11-21

|Col1|Col2|Col3|Col4|Col5|Col6|Col7|f Timestamp|Col9|Col10|Col11|
|---|---|---|---|---|---|---|---|---|---|---|
||||Labeled A|no|maly|Score o|f Timestamp||||
||||||||||||
||||||||||||


Time

Figure 7: Additional anomaly detection results on the test set of PMU-C.


16

14

12

2016-12-03 2016-12-05 2016-12-07 2016-12-09

|Col1|Col2|Col3|Col4|Col5|Col6|Col7|Col8|Col9|Col10|Col11|Col12|Col13|Col14|Col15|Col16|Col17|Col18|Col19|Col20|Col21|
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|||||||Labeled|Anomaly||S|c|ore|of Ti|mesta|mp|||||||
||||||||||||||||||||||
||||||||||||||||||||||


Time


16

14

12

2016-12-10 2016-12-12 2016-12-14 2016-12-16

|Col1|Col2|Col3|Col4|Col5|Col6|Col7|Col8|Col9|Col10|Col11|Col12|Col13|Col14|Col15|Col16|Col17|Col18|Col19|Col20|Col21|Col22|
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|||||||Labele|d|Anomaly||Score|of|Ti|me|st|a|mp||||||
|||||||||||||||||||||||


Time


16

14

12

2016-12-16 2016-12-18 2016-12-20 2016-12-22

|Col1|Col2|Col3|Col4|Col5|Col6|Col7|Col8|Col9|Col10|Col11|Col12|Col13|Col14|Col15|Col16|
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|||||||Labeled A|no|maly|S|core of Ti|mestamp|||||
|||||||||||||||||


Time


Figure 8: Anomaly detection results on the test set of PMU-B.


-----

ADDITIONAL RESULTS FOR DAG EVOLUTION

(a) Four-month shift (b) Five-month shift

(c) Six-month shift (d) Seven-month shift

(e) Eight-month shift (f) Nine-month shift

(g) Ten-month shift (h) Eleven-month shift

Figure 9: Evolution of the learned DAG on PMU-B over time.


-----

