Anomaly Detection of Time Series with Smoothness-Inducing Sequential Variational Auto-Encoder

Due to different reasons such as faulty sensors or extreme events, the real-world time series data are often contaminated with anomalies. Consider we have time series data $\mathbf{x}_{1:T}=(\mathbf{x}_{1},\mathbf{x}_{t},\dots,\mathbf{x}_{t})$ contaminated with anomalies, and each $\mathbf{x}_{t}\in\mathbb{R}^{M}$ (indicating $T$-dimensional time series) is modeled by:

where $\bm{f}(t)$ is the true signal, $\bm{\epsilon}_{t}$ is non-stationary noise of observations, $\bm{I}$ is a $M\times T$ binary matrix, which is the indicator matrix of anomalies. The value $\bm{I}_{m,t}=1$ when there is an anomaly at observation $\mathbf{x}_{m,t}$, and $\bm{e}_{t}$ is the anomaly vector, and $\odot$ is element-wise multiplication. Our goal is to recover the true signal $\bm{f}(t)$ and detect anomaly $\bm{I}$ from $\mathbf{x}_{1:T}$.

For the unsupervised learning problem, to recover true signal $\bm{f}(t)$ from anomaly-contaminated data, it is important to impose inductive bias or prior knowledge to the model as suggested in [17, 18]. For time series modeling, we first introduce smoothness prior modeling of time series as described in [31]. Consider the time series $y_{1:T}$ (for convenience, we deal with the uni-variate case), it is assumed to consist of the sum of a function $f$ and observation noise $\epsilon$:

where $f\in\mathcal{F}$ is an unknown smooth function, and $\epsilon\sim\mathcal{N}(0,\sigma^{2})$. The problem is to estimate $f$ from the observations. A common approach is to use a parametric model such as polynomial regression or neural network to approximate $f$. The quality of estimation is dependent upon the appropriateness of the assumed model class. The idea is to balance a trade-off of goodness-of-fit to the data, and a smoothness criterion. The solution is to minimize the objective with an appropriately chosen smoothness trade-off parameter $\lambda$:

where $\nabla^{k}$ is a $k$-th order differential constraint on the solution $f$, with $\nabla f(t)=f_{t}-f_{t-1}$, $\nabla^{2}=\nabla(\nabla(f(t)))$, etc.

The solution can be solved from a Bayesian perspective: consider a homoscedastic (variance is constant over time) time series $y_{t}\sim\mathcal{N}\mathcal{(}f(t),\sigma^{2})$ with constant variance $\sigma^{2}$. The solution becomes maximizing:

Its Bayesian interpretation can be written by:

where $\pi(f|\lambda,\sigma^{2},k)$ is the prior distribution of $f$, and $p(y|\sigma^{2},f)$ is data distribution conditioned on $\sigma^{2}$ and $\pi(f|y,\lambda,\sigma^{2},k)$ is the posterior of $f$.

Consider a dataset $\mathbf{x}=\{\mathbf{x}^{(i)}\}_{i=1}^{N}$ consisting of $N$ i.i.d data points $\mathbf{x}^{(i)}$. From a latent variable modeling perspective, we assume the data are generated by some random processes that capture the variations in the observed variables $\mathbf{x}$. Thus we have the integral of marginal likelihood:

where $\theta$ is generative model parameters, and $p_{\theta}(\mathbf{z})$ is the prior distribution over latent variable $\mathbf{z}$, and $p_{\theta}(\mathbf{x}|\mathbf{z})$ is the conditional distribution that generates data $\mathbf{x}$. Variational Auto-Encoder (VAE) [27] models the conditional probability $p_{\theta}(\mathbf{x}|\mathbf{z})$ with a neural network as a flexible approximator. However, introducing a complex non-linear mapping for the generation process from $\mathbf{z}$ to $\mathbf{x}$ results in intractable inference of posterior $p_{\theta}(\mathbf{z}|\mathbf{x})$. Hence VAE uses a variational approximation $q_{\phi}(\mathbf{z}|\mathbf{x})$ of the posterior, i.e. the evidence lower bound (ELBO) of the marginal likelihood as written by:

where the first R.H.S. term is the Kullback-Leibler divergence between two probability distributions $Q$ and $P$, and $\phi$ denotes variational parameters. The generative model and inference model are trained jointly by maximizing the ELBO w.r.t their parameters. The KL-divergence $D_{KL}(q_{\phi}(\mathbf{z}|\mathbf{x})\rVert p_{\theta}(\mathbf{z}))$ of Eq. 7 can often be integrated analytically. However, the gradient of the second term is problematic because latent variable $\mathbf{z}\sim q_{\phi}(\mathbf{z}|\mathbf{x})$ is stochastic which blocks backpropagation. The Stochastic Gradient Variational Bayes (SGVB) [32] expresses the random variable $\mathbf{z}$ as a deterministic variable $\mathbf{z}=g_{\phi}(\bm{\epsilon},\mathbf{x})$, where $\bm{\epsilon}$ is an auxiliary variable with independent marginal distribution $p(\bm{\epsilon})$. Then we have the SGVB estimator $\tilde{\mathcal{L}}(\theta,\phi;\mathbf{x})\simeq\mathcal{L}(\theta,\phi;\mathbf{x})$ of the lower bound:

Learning is performed by optimizing the SVGB estimator with backpropagation. VAE [27] uses centered isotropic multivariate Gaussian $p_{\theta}(\mathbf{z})=\mathcal{N}(\mathbf{z};\bm{0},\bm{I})$ as prior of latent variable $\mathbf{z}$. The variational approximate posterior $q_{\phi}(\mathbf{z}|\mathbf{x})$ is a multivariate Gaussian with diagonal covariance matrix: $\mathcal{N}(\bm{\mu},\text{diag}(\bm{\sigma}^{2}))$, and its mean $\bm{\mu}$ and variance $\bm{\sigma}^{2}$ are parameterized by a neural network. For $p_{\theta}(\mathbf{x}|\mathbf{z})$, a Bernoulli distribution for binary data or multivariate Gaussian distribution (as adopted in our model) for continuous data is often used.

For time series modeling, many prior publications have extended VAE in a sequential manner [33, 34, 16, 35]. The main idea is to use a deterministic Recurrent Neural Network (RNN) as a backbone, represented by the sequence of RNN hidden states. At each timestamp, the RNN hidden state is conditioned on the previous observation $\mathbf{x}_{t-1}$ and stochastic latent variable $\mathbf{z}_{t-1}$, such that the model is able to capture sequential features from time series data.

Let $\mathbf{x}_{1:T}=(\mathbf{x}_{1},\mathbf{x}_{t},\dots,\mathbf{x}_{t})$ denote a multidimensional time series, such as observation data collected by multiple sensors of T time slots. Our goal is to learn the data distribution of the training time series $p(\mathbf{x}_{1:T})$. The model is composed of two parts, the generative model and the inference model.

As discussed before, we seek to devise suitable inductive bias for anomaly detection of time series in the sequential VAE framework. To this end, we begin with the factorization of the existing loss function. Specifically, we substitute the numerator and denominator in Eq. 14 with Eq. 9 and Eq. 15 respectively, the objective of Sequential VAE becomes:

The first term is the loss of inference model, which can be viewed as a regularizer, encouraging the model to learn disentangled feature by putting a prior over latent variable [37]. The second term is the reconstruction loss of the generative model. Suppose $x_{m,t}=\bm{f}_{m}(t)+\epsilon_{m,t}+\bm{e}_{m,t}$ is an anomaly, and its associated reconstruction distribution is $\hat{x}_{m,t}=\mathcal{N}(\mu_{m,t},\sigma_{m,t})$. The generative model tries to reconstruct the contaminated $x_{m,t}$, rather than the true signal $\bm{f}_{m}(t)$, otherwise the reconstruct loss, which is the negative log likelihood $-p_{\theta}(x_{m,t}|\hat{x}_{m,t})$ will be high. To minimize the objective loss function from the data with anomalies, the model will either learn a biased mean $\mu_{m,t}$, or an over-estimated variance $\sigma_{m,t}$. As a result, the Sequential VAE can be vulnerable to anomalies.

To address account for such a bias caused by the anomalies, one natural idea is to induce smoothness prior to the generative model. Recall in Section III-B, the smoothness criterion refers to the accumulative difference of signal $\sum_{t=1}^{T}\nabla f_{t}$. One might be tempted to simply add such a smoothness regularization to the objective. However, since the Sequential VAE is under probabilistic framework, such a forced combination is not mathematically coherent. However, the classical approach treat observation as deterministic variable, however, in sequential VAE, we treat each observation $\mathbf{x}_{t}$ as a random variable, and we estimate the probability distribution for each observation $p(\mathbf{x}_{t})\sim\mathcal{N}(\bm{\mu}_{t},\text{diag}(\bm{\sigma})_{t}^{2})$, such that the classical regularizer failed to work in probabilistic framework, which calls for more comprehensive method.

We extend the smoothness regularization to probabilistic framework, by following the common assumption that the probability density function over time would vary smoothly over time. Among all possible reconstructed series $\hat{\mathbf{x}}_{1:T}\in\mathcal{X}$, the true signal should have a smooth transition of distribution. To construct a valid regularizer, we need to measure the smoothness of a time-varying probability distribution.

For two consecutive data points of a single time series: $x_{m,t-1}$ and $x_{m,t}$, associated with reconstructions $p_{\theta}(\hat{x}_{m,t-1})$ and $p_{\theta}(\hat{x}_{m,t})$. Let $d(\cdot,\cdot)$ to be a distance metric between two distributions, we propose accumulative transition cost of time-varying probability distributions:

In this paper, we choose KL-divergence as the distance metric between distributions. For Isotropic Gaussian case in this paper, the transition cost of two Normal distributions is:

The final objective involves: 1) encoding loss for inference model; 2) expected reconstruction error for generative model; 3) the proposed variational smoothness reguarlizer. We apply stochastic gradient variational Bayes (SGVB) estimator to approximate the expected reconstruction error $\mathbb{E}_{q_{\phi}(\mathbf{z}_{\leq T}|\mathbf{x}_{\leq T})}[\log p_{\theta}(\mathbf{x}_{t}|\mathbf{z}_{\leq t},\mathbf{x}_{<t})]$. The learning objective is:

where $\mathbf{z}_{t}=\bm{\mu}_{z,t}+\bm{\sigma}_{z,t}\odot\bm{\epsilon}$ and $\bm{\epsilon}\sim\mathcal{N}(0,\bm{I})$, and $\lambda$ is a smoothness regularization hyper-parameter. We set $p_{\theta}(\hat{\mathbf{x}}_{t}|\mathbf{z}_{\leq t-1},\mathbf{x}_{<t-1})=p_{\theta}(\hat{\mathbf{x}}_{t}|\mathbf{z}_{\leq t},\mathbf{x}_{<t})$ when $t=1$.

The next step is to compute anomaly score for each data point $\mathbf{x}_{m,t}$. For a trained model, input a chunk of time series $\mathcal{X}=\{\mathbf{x}_{t}\}_{1:W}$, each $\mathbf{x}$ is a $M$-dimensional vector, the model first encode it into $\mathbf{z}_{1:W}$, then decode to a sequence of random variables $\hat{\mathcal{X}}=\{\mathcal{N}(\bm{\mu}_{t},\text{diag}(\bm{\sigma}_{t}))\}_{1:W}$. The whole encoding and decoding process is called probabilistic reconstruction.

Recall that in Eq. 1 of Section III-A, we use $\bm{I}$ to indicate the indicator matrix of anomaly labels. In this step, we use $\bm{A}$ to indicate the matrix of anomaly scores.

One natural choice is to use absolute reconstruction error, like anomaly detection with deterministic auto-encoders [40]. The anomaly score $e_{m,t}$ for observation $x_{m,t}$ is computed as:

where $\mu_{m,t}$ is the $m$-th dimension of reconstructed mean $\bm{\mu}_{t}$. Intuitively, a large reconstruction error indicates anomaly behavior of time series, which is a common assumption for many reconstruction based anomaly detection algorithms. The underlying assumption of absolute reconstruction error is that the variance is stationary, which means that the variance is constant over time. However, in many real-world data, the variance is non-stationary, and it varies over time. We illustrate this using an example as shown in Fig. 3. We generate some time series data with Gaussian process as true signal, denoted as the real blue line centering scatters, then we add some Gaussian noise $\mathcal{N}(0,\sigma_{t})$, where $\sigma_{t}$ increases over time. The black dots is the noisy observation, and blue transparent regions increasing level of noise. Then, we add two points, point A and point B, they both deviate from the true signal with distance 1. We see that, point A is more likely to be an anomaly point, while point B is more likely to be a normal point that located at some noisy period. In this case, if we use absolute reconstruction error, we may fail to detect anomaly point A, or otherwise report normal point B as an anomaly.

Due to the non-stationary nature of real-world time series, we need a robust method for computing anomaly scores. Previous VAE-based anomaly detection methods have introduced reconstruction probability based anomaly score, here we devise a reconstruction probability for sequential Auto-encoders. Specifically, the probability density of input $\mathbf{x}$ is:

Draw $L$ samples from $q_{\phi}(\mathbf{z}|\mathbf{x})$, and the anomaly score is negative of average of $L$ probability. For numerical stability, we often work with $\log$ probability, and we have

However, for Sequential VAE, this can be more challenging. Because $\mathbf{z}$ is not independent, but has sequential dependencies. We choose to adopt Sequential Monte Carlo (SMC) [41], which is widely used in state-space models, such as Bayesian filters [42], and dynamical systems [43].

We conclude the SMC for computing anomaly scores and detecting anomalies in Algorithm 2. For the full multivariate time series matrix, $\bm{X}\in\mathbb{R}^{M\times T}$, we slice them into non-overlapped chunks, and reconstruct each of them individually.

Performance Metrics. For time series data $\mathbf{X}\in\mathbb{R}^{M\times T}$, the model outputs anomaly score matrix $\bm{A}\in\mathbb{R}^{M\times T}$, and each item $A_{m,t}$ is the likelihood of data point $x_{m,t}$ being an anomaly. For detection, we set a threshold $\alpha$ to obtain the final anomalies label indicator matrix $\bm{I}_{\alpha}\coloneqq\bm{A}>\alpha$. Then we compute precision, recall, true positive rate (TPR), false positive rate (FPR) at time series point level.

Trade-off is needed as the detection characteristics vary by different threshold $\alpha$. When threshold is tight, the precision will be higher and recall will be lower. We compute area under receiver operating characteristic (AUROC), area under precision recall curve (AUPRC) and the best F1-score. AUROC quantifies the FPR and TPR with different $\alpha$, and AUPRC quantifies the precision and recall with different $\alpha$. Note that all three metrics are independent of specific threshold.

There are inherent differences between AUROC and AUPRC. For anomaly detection tasks where class (anomaly or normal data points) distributions are heavily imbalances, the AUROC metric is less insensitive to false positives than the AUPRC metric. As discussed in [44, 45], the AUPRC score is more discriminative than AUROC score. Previous work may choose one of the metrics to evaluate the performance. To comprehensively evaluate the model’s performance, we report both AUROC and AUPRC scores, but we mainly focus on analysis AUPRC in this paper.

Peer methods. Recent deep learning based methods are compared, including Donut and STORN.

• •

  Donut [3]. The model proposed very recently is based on VAE, which takes a fixed length sliding windows of an uni-variate time series as input. To capture ‘normal patterns’ of data, the modified ELBO technique is devised to leverage generative model $p_{\theta}(\mathbf{x}|\mathbf{z})$ and latent regularization $p_{\theta}(\mathbf{z})$ of VAE, such that the model is encouraged to learn only normal patterns.

• •

  STORN [30]. As a sequential version of VAE, the model is based on stochastic recurrent neural network (STORN). Different from SISVAE, the sequences of latent variables $\mathbf{z}_{t}$ are generated independently. The model uses lower bound output as anomaly score.

We also include three traditional statistical models.

• •

  History Average (HA). We use absolute deviation from history average value as anomaly score. We include this model in comparison to show how a very simple model perform within different datasets.

• •

  Auto-regressive Moving Average (ARMA) [46]. ARMA is a classic time series analysis model [47]. We use absolute fitting residuals as anomaly score.

• •

  Linear Dynamical System (LDS). This model is also known as Kalman filter. Different from deep latent variable based state space models, the transition function and emission probability are linear (see Table II).

We further include some variants of our model, to study the contribution and behavior of each component in our model.

• •

  SISVAE-p Our recommended model, which uses the sliding windows of a time series matrix as input, similar to Donut. This model is more applicable for long sequences.

• •

  SISVAE-e. Different from SISVAE-p, it takes absolute reconstruction error as anomaly score. We use this model to show the performance gap by different detection criteria.

• •

  SISVAE-0 To show the effect of our smoothness regularizer, we let SISVAE-0 to be SISVAE-p for $\lambda=0$ .

• •

  VAE-s. To test the contribution of the proposed smoothness prior technique, we add the regularization to vanilla VAE as a comparison named VAE-s. The model can be viewed as our SISVAE model without latent temporal structure. We use this model to show the effectiveness of introducing latent temporal structure.

In particular, we compare the characteristics of peer methods in Table II, whereby both traditional methods i.e. History Average (HA), Autoregressive Moving Average (ARMA), Kalman Filter (KMF) and deep latent model i.e. Donut, STORN, SISVAE, SISVAE-mb are included.

Implementation details For preprocessing, we set window size $W$ to be 120 for each dataset. The sequential VAE based models, STORN, and SISVAE have a similar network architecture, which is defined by two varying hyper-parameters, h_dim for dimension of hidden layer, and z_dim for dimension of latent variable $\mathbf{z}$. We set h_dim to be 200, and z_dim to be 40 for all datasets. We empirically set the smoothness regularization hyper-parameter $\lambda=0.5$ for SISVAE during training process. The encoder, decoder, and RNN of models use the same hidden dimension. The three models are implemented using PyTorch [48] 0.4, and trained on a single Nvidia RTX 2080ti GPU. We train the three models 200 epochs using Adam [49] with learning rate 0.001. Note for Donut, the M-ELBO technique proposed by Donut may require some anomaly labels to learn ‘normal’ patterns. For fair comparison, we provide no anomaly label to the model, and the experiments are conducted under fully unsupervised setting. We use the code open sourced by the authors. We set h_dim to be 200, and z_dim to be 40, same as our model. The model is fundamentally designated for uni-variate time series and it is nontrivial to extend it to the multi-dimensional case. We test two settings: 1) train a model for each sequence individually, and compute anomaly score individually; 2) flatten the multi-dimensional time series data after preprocessing to accommodate the model, the anomaly score is computed together. We find the second protocol exhibits better performance. This may be because the normal patterns are common along different sequences. With more sequences as input, the model can be less influenced by anomalies. We adopt the second protocol for Donut in the experiments. For all deep latent variable models, we use negative reconstruction probability as anomaly score output. For all three traditional statistical models, we use the same preprocessing method. For history average (HA), the anomaly score is absolute observation value since all series have zero mean. For ARMA and LDS, we use absolute reconstruction error as anomaly score.

In this section, we aim to gain insights into the behavior of our proposed SISVAE model. We investigate the anomaly detection performance of our model for different proportions of anomalies. Furthermore, we are interested in the anomaly detection performance under different settings of smoothness regularization hyper-parameter $\lambda$.

Referring the synthetic data generation protocol in [24], we generate 100 time series sequences of length 200 from a multi-output Gaussian process with a multiple output kernel [50], implemented by a python package named GPy [51]. Sequences are correlated through a random sampled positive definite coregionalization matrix, and we add random non-stationary Gaussian noise on the 100-dimensional sequence of length 200. The location of anomalies is simulated by a binary mask with Bernoulli distribution of probability {0.5%, 1%, 2%, 5%, 10%} for different fractions of anomalies. We insert Poisson noise where the binary mask is equal to one as anomalies, and the intensity $\lambda_{t}$ of each anomaly is proportional to the amplitude of inserted place $x_{t}$. The plus and minus sign is random sampled from Bernoulli distribution with $p=0.5$.

Effects of anomaly proportion. Model performance variation with respect to the proportion of anomalies is shown in Fig. 4(a). Since all the models are unsupervised, they are inevitably influenced by the anomalies, because a full reconstruction of both normal and anomaly data is encouraged in the objective function. The proposed variational smoothness regularizer places penalties at consecutive outputs of the generative model that violate the temporal smoothness assumption. As observed from the result, as the anomaly proportion grows, the AUPRC scores of two un-regularized model SISVAE-0 and STORN drop significantly, while the regularized models’ performance almost keeps maintained. SISVAE-0, the proposed SISVAE model shows superior performance over competing methods in all settings of anomaly proportion.

Effects of regularization hyper-parameter. The variational smoothness hyper-parameter $\lambda$ in the objective by Eq. 17 determines the penalty of non-smooth output of the generative model. We analyze the detection performance of SISVAE model under different settings of hyper-parameters. We train the model with $\lambda$ = {0.01,0.05,0.1,0.2,0.5,0.8,1.0} using synthetic dataset with 2% anomaly proportion. Figure 4(b) shows the change of AUPRC performance w.r.t $\lambda$, where the model achieves best performance when $\lambda$ is set to be 0.5.

Convergence study. We analyze the contribution of our regularizer to convergence property. On the synthetic 2% anomaly proportion dataset, we train three models: i) SISVAE-p without regularization ($\lambda=0$). ii) classical smoothness regularizer of deterministic model where only the mean is regularized which is formulated as: $\lambda\sum_{t=1}^{T}\sum_{m=1}^{M}(\nabla^{2}_{t}\hat{x}_{m,t})^{2}$. iii) proposed KL regularizer. We record the AUPRC score over training. As shown in Fig. 4(c), at the first 100 iterations, all the three models are trained very fast, and the detection performance increases stably. This suggests that the models are learning the major ‘normal’ patterns of the data, while the anomalies are temporally neglected. However, after 100 iterations, the performance of the un-regularized and mean-regularized models fluctuate severely between 0.45 and 0.85, and our KL-regularized model (i.e. SISVAE-p) also fluctuates but in a relatively small interval. We think the fluctuation is caused by contradiction between generative (decoder) model and recognition (encoder) model: the generative model is trying hard to reconstruct every data point perfectly, while the inference model aims to map the observed anomaly-contaminated data points into a simpler low-dimensional latent space. The decoder and encoder are playing a seesaw battle in the training procedure, thus causing a sawtooth shape AUPRC score. Under the probabilistic framework, classical mean-regularizer which is in its nature deterministic, cannot work well. Its convergence curve is akin to the un-regularzied model. In contrast, the proposed variational smoothness regularizer weakens the reconstruction objective, and guides the decoder to generate sequences that obey temporal smoothness assumption, preventing the model from overfitting ‘anomaly’ patterns. As shown in Fig. 4(d), the reconstruction loss and inference loss of regularized model converge stably, while the training of un-regularized model is unstable. Hence the KL-regularized SISVAE model reaches a higher AUPRC score.

Datasets This experiment involves the following datasets, we conclude statistics of datasets in Table III.

• •

  Yahoo’s S5 Webscope Dataset: This dataset is a part of Yahoo’s Webscope program²²2https://webscope.sandbox.yahoo.com/. It consists of 367 time series. This dataset consists of four different classes A1/A2/A3/A4 with 65/100/100/100 time series respectively. Class A1 has real data collected from computational systems, which records real production traffic of Yahoo’s website, and the anomalies are labelled by human expert manually. Classes A2, A3 and A4 contain synthetic data, whereby anomalies are constructed and are inserted at random positions.

• •

  Open $\mu$PMU power network dataset: This dataset is collected from a power distribution equipped in Alameda, CA, USA with advanced smart meters called phasor measurement units ($\mu$PMU)³³3https://powerdata.lbl.gov [52]. There are three $\mu$PMU devices, each recording three-phase voltage and current magnitude and phase angle. Anomalies caused by voltage disturbance are labeled by human experts. We collect measurement data and anomaly labels during the whole day of June 3 with 20 seconds resolution. We use 36-dimensional time series of length 1080.

Preprocessing In real application, different time series may have different scales, e.g. temperature and air pressure sensors have different scales. So we follow the common practice to normalize each sequence to standard normal distribution.

Threshold based anomaly detection. First, we conduct threshold based anomaly detection experiment. For each model, we collect anomaly score matrix $\bm{A}$ for each model, then we compute AUROC, AUPRC and F1-score for each possible threshold $\alpha$. We select the highest F1-score the model could achieve as the final F1-score. We compare all the peer methods in this experiment, whose results are shown in Table IV. SISVAE-p outperforms other models in most datasets and metrics. For SISVAE-p and SISVAE-e, we find that our proposed SMC based reconstruction probability for computing anomaly score notably improves the performance. Among all the five datasets, we find SISVAE-p and SISVAE-e perform similar on A2 dataset. We further find that the A2 dataset is in general stationary regarding with variance, which means that modeling the time-varying variance is important for anomaly detection. For traditional statistical models, we find the simple method HA shows decent performance in some datasets, however, the potential of HA is very limited. We are surprised that ARMA failed on A3, A4 and $\mu$PMU dataset. Our further study uncovers that these three datasets are highly non-stationary with both mean and variance, which violate the basic stationary assumption of ARMA. In contrast, we find LDS performs competitively on A3, A4 and $\mu$PMU comparing to the best model, however, it fails on dataset A1. This suggests that linear transition and linear emission probability of LDS may not be able to capture complex non-linear relationship. For VAE-s, the proposed smoothness prior regularizer works even without a temporal structure model. In our analysis, the regularizer constructs inductive bias of the true signal, leading to robustness to unknown anomaly.

As discussed, we find models’ anomaly detection performance varies across datasets, such as LDS and ARMA. To make the overall performance comparison more intuitionistic, we compute the rank of models for each dataset and performance metric, and compute the mean and standard deviation of ranks for each models. The plots of rankings are shown in Fig. 5 showing our model outperforms other models notably.

Ranking based anomaly detection. Due to limited budget, sometimes we can only deal with a limited number of detected anomalies hence the precision need be high. In this context, precision@K is a useful metric, which is widely used in information retrieval [53]. For each model, we select data points that are assigned with top $K$ highest anomaly score, and precision@K is computed by ${V}/{K}$, where $V$ is the number of true anomalies in $K$ predictions. We examine models with $K={10,50,200}$, and the results are shown in Fig. 8. We can see that, when $K$ is 10, STORN, LDS and SISVAE-p perform well. However, Donut fails on A2, A3 and A4 dataset. The reason is that Donut places very high anomaly score on some normal data points. Overall, SISVAE-p’s precision keeps stable as $K$ increases. This means SISVAE-p makes more precise detection when the number of detection trials is limited.

Model characteristics. The detection threshold $\alpha$ is often set according to specific needs. In order to investigate the detailed detection characteristic, we plot the precision-recall curve (PRC), and receiver operating characteristic curve (ROC), as shown in Fig. 6 and Fig. 7. Comparing to LDS and STORN, SISVAE-p outperforms at most thresholds. For precision and recall curve in Yahoo A2, A3 and A4 dataset, note that Donuts performs strangely, where the precision and recall curves do not monotonically decrease as expected, but increase first and then decrease. We find this is because of Donut places very large anomaly scores to some normal points, while the other three models do not have this problem. This means modeling latent temporal dependencies is important for time series anomaly detection tasks.

The trade-off between recall and precision. Despite the fact that SISVAE-p is better than baselines in most cases, however, the preference of detection metric also affects the choice of model. It is important to know when to use which model. The cost of validating a detected anomaly can be expensive which calls for model’s high precision. For large-scale power grid system, it is required the precision need be strictly above 0.95. In this case, SISVAE-p outperforms other methods in recall. In contrast, for fraud detection, recall can be more important. If we set a hard bar for recall to 0.95, we find Donut is also competitive, it outperforms SISVAE-p in A2 and A4 dataset w.r.t. precision.

We compare the three deep generative models, SISVAE-p, STORN, and Donut with real-world examples from Yahoo-A1 dataset. For each model, we choose to set threshold $\alpha$ when F1-score is maximized. The dataset consists of 65 sequences each with length 1420. In what follows we show some representative time series fragments for discussion.

False alarms on normal time series. Figure 9(a) shows a normal sequence without any anomaly. Since the training data is contaminated with anomalies, the model can be influenced by anomalies and fails to capture normal patterns such that false alarms can be be reported. Figure 9(b) illustrates the density estimation of the three models, and Fig. 9(c) shows the anomaly score generated by three models. Note STORN reports three false alarms, while SISVAE-p and Donut do not. Moreover, we find the anomaly score generated by SISVAE-p is stably low, while Donut is more fluctuating given noisy data. The results demonstrate that SISVAE-p is robust to contaminated data for unsupervised training.

Detection of point-level anomalies. Figure 10(a) includes a sequence from Yahoo-A1 dataset with anomaly at around time point 300. This type of anomaly is often called additive outlier [54] or extreme value [55]. As shown in Fig. 10(b), density estimation of STORN and Donut are corrupted by the point anomaly, such that they place much probability density to the anomaly point. In contrast, SISVAE-p successfully recovers the true underlying density. Figure 10(c) shows the anomaly score generated by the three models, and we can see that all the models identify the anomaly point. However, STORN and Donut report some false alarms. Moreover, the anomaly score of SISVAE-p is more smooth than those of other models, suggesting that SISVAE-p is more robust to the moderate fluctuation of normal data. Overall, SISVAE is capable of detecting point-level anomalies robustly.

Detection of sub-sequence level anomalies. In some cases, the anomalies is in the presence of sub-sequences, such as anomaly cardiac cycle in ECG data [56]. Figure 11(a) shows a sequence with two anomaly sub-sequences, whereby the red shaded region covers the true anomaly. Figure 11(b) shows the density estimation, whereby Donut and STORN are notably influenced by the anomaly sub-sequence. In contrast, SISVAE-p estimates a smooth probability density over time, which is more likely to be the true density. Figure 11(c) shows SISVAE-p successfully detects most of the anomaly sub-sequences, while STORN and Donut suffer from false negative and false positives. The results demonstrate the effectiveness of the proposed smoothness prior regularizer, helping recover the true density and improve detection performance.

Simulation study on uni-variate time series. As discussed in Section V-A, Donut is designated for uni-variate time series. In order to investigate the detection performance in one-dimensional case, we conduct simulation study on data generated from Mackey-Glass equation [57], which serves a benchmark for nonlinear models since the true underlying function is nonlinear [58]. The time series is generated as:

with $\gamma=0.1,\beta=10,\alpha=0.2$. We run the Mackey-Glass equation 5000 iterations, then we collect an uni-variate time series with length 5000. Two types of anomalies are added:

• •

  Point level anomaly: we insert point anomalies at $0.3\%$ time step randomly, each point anomaly is a combination of: i) Poisson noise with $\lambda=1$; ii) Gaussian noise with zero mean and unit variance; iii) 0.2 fixed bias.

• •

  Sub-sequence level anomaly: two sub-sequences anomalies are inserted, starting at random time step, with random length sampled from even distribution of range $10\sim 28$. For each anomaly sub-sequence, true signal is replaced by a sequence sampled from a Gaussian process with RBF kernel of unit variance, and 0.3 length-scale.

The generated Mackey-Glass time series is shown in Fig. 12. To illustrate the generated anomalies, we zoom in three fragments out of the whole time series.

In this experiment, we find Donut shows better performance when z_dim is lower. We enumerate Donut’s performance on $\text{z\_dim}=\{2\dots 10\}$, and we find Donut performs best when $\text{z\_dim}=5$. For SISVAE-p, we find the performance is relatively stable when we change z_dim. Hence we use the same setup as discussed in Section V-A. The ROC and PRC are shown in Fig. 13. In conclusion, benefiting from the sequential network structure of SISVAE, and the proposed smoothness prior regularizer, SISVAE-p exhibits better anomaly detection performance for uni-variate time series.