A new distributional treatment for time series anomaly detection

Abstract

Time series is traditionally treated with two main approaches, i.e., the time domain approach and the frequency domain approach. These approaches must rely on a sliding window so that time-shift versions of a sequence can be measured to be similar. Coupled with the use of a root point-to-point measure, existing methods often have quadratic time complexity. We offer the third $\mathbb {R}$ domain approach. It begins with an insight that sequences in a stationary time series can be treated as sets of independent and identically distributed (iid) points generated from an unknown distribution in $\mathbb {R}$. This $\mathbb {R}$ domain treatment enables two new possibilities: (a) The similarity between two sequences can be computed using a distributional measure such as Wasserstein distance (WD), kernel mean embedding or isolation distributional kernel ($\mathcal {K}_I$), and (b) these distributional measures become non-sliding-window-based. Together, they offer an alternative that has more effective similarity measurements and runs significantly faster than the point-to-point and sliding-window-based measures. Our empirical evaluation shows that $\mathcal {K}_I$ is an effective and efficient distributional measure for time series; and $\mathcal {K}_I$-based detectors have better detection accuracy than existing detectors in two tasks: (i) anomalous sequence detection in a stationary time series and (ii) anomalous time series detection in a dataset of non-stationary time series. The insight makes underutilized “old things new again” which gives existing distributional measures and anomaly detectors a new life in time series anomaly detection that would otherwise be impossible.

Appendix

The machine used in the experiments has one Intel Core i7-12700K CPU with 128 GB memory and one NVIDIA GeForce RTX 3060 GPU with 12 GB memory.

1.1 Experiment settings for anomalous sequence detection

Sequences in each time series are preprocessed with z-score normalization. All final scores output by detectors are normalized in [0, 1]. For detectors that rely on randomization, we report the average result of 10 trials on each time series.

[Datasets] Synthetic datasets noisy_sine and ARMA are originated from a previous work [23]. Real-world datasets include MIT-BIH Supraventricular Arrhythmia Database (MBA) [18, 32], and other datasets from various domains have been studied in earlier works [23, 25, 48].

Some datasets have two versions, e.g., ann_gun and stdb_308; and each version uses one of the two variables. When either version produces similar AUC for most detectors, we have chosen to use one only. Some datasets are trivial, e.g., chfdb_chf0175 and qtdbsel102; and all detectors have the perfect result ($\hbox {AUC}=1$), so we do not show them in Table 7.

We labeled anomalous periods for each time series following the previous work [23, 25, 48]. Details are given in Table 17. Positions of anomalies in MBA datasets can be seen in folder “MBA_Annotation."

Table 17 Locations of anomalous periodic sequences in each dataset in terms of index i in $X_i$

The period of some datasets varies slightly at different time steps in the series; but it has no effect on the detection accuracy of all algorithms. Our algorithm works well when the sequence length is set to be roughly the length of the period.

Brief descriptions of some datasets are given as follows.

dutch_pwrdemand: There are a total of 6 anomalous weeks. Some papers [5, 23, 25] use this dataset with fewer anomalies because they treat continuous anomalous weeks as one anomaly.

ann_gun: It has only one anomalous period when it was first used in Keogh’s work [25], as shown in Fig. 11a. Other anomalous periods in this dataset were later identified [5], and they are shown in Fig. 11b.

Patient_respiration: Like the previous work [23], we use the subset that begins at 15,500 and ends at 22,000 from the nprs44 dataset [25]. There are one apparent anomaly and one subtle anomaly in this dataset as shown in Fig. 12a.

TEK: Following the previous work [23], we also concatenate dataset TEK14, TEK16 and TEK17 as TEK of length 15,000. In Keogh’s work [25], a total of 4 anomalies are marked. But TEK14 has 2 anomalous snippets belonging to the same period as shown in Fig. 12b. Since we regard each anomaly as an anomalous sequence of one complete period, it is treated as one anomalous periodic sequence of length 1000. So there are a total of 3 anomalous sequences in our annotations of this dataset.

MBA803, MBA805, MBA806, MBA820 and MBA14046: These datasets are subsets of the full MBA dataset, as used in the previous work [5].

Fig. 11

a One anomaly period; and b additional anomalous periods as identified by [5] on the ann_gun dataset. The diagrams are extracted from [5, 25], respectively

Fig. 12

Anomalies in a Patient_respiration; b a period of TEK14. The diagrams are extracted from [25]

[Algorithms] The STOMP [66] implementation of MP is used; NormA is from http://helios.mi.parisdescartes.fr/~themisp/norma/;

$\mathcal {K}_I$-based detectors are our implementations based on [59]; and WFD is from https://github.com/GAMES-UChile/Wasserstein-Fourier. Others are from scikit-learn.org. All are in Python.

The parameter search ranges of all algorithms used are given in Table 18.

As for the 1Line method, we use one of the following five types of basic vectorized primitive functions in Matlab as an anomaly score for each sliding window of size $\omega$:

• I. ±diff(X): the difference between the current point and the previous point. Here $\omega =1$.

• II. ±movmax(X, ${{\varvec{\omega }}}$)

• III. ±movmin(X, ${{\varvec{\omega }}}$)

• IV. ±movmean(X, ${{\varvec{\omega }}}$)

• V. ±movstd(X, ${{\varvec{\omega }}}$)

where X is the time series; and the maximum, minimum, mean or standard deviation is computed for each window of ${{\varvec{\omega }}}$ points.

We run these 5 one-liners on each dataset and report the median AUC (out of the five values) in Table 7. Low median values indicate that the datasets are hard to detect using the 1Line method; otherwise, the datasets have anomalies that can be easily detected.

Table 18 Parameter search ranges

[Measures] The detection accuracy of an anomaly detector is measured in terms of AUC (area under ROC curve). As all the anomaly detectors are unsupervised learners, all models are trained with the given datasets with no labels. Only after the trained models have made predictions, the ground truth labels are used to compute the AUC for each dataset.

Given a periodic time series X of length n and period length m, a sequence $X_i$ of X is a subset of contiguous values of length m, for $i=1,\ldots ,s$, where $s= \lfloor \hbox {length}(X)/m\rfloor$. A distribution-based (non-sliding-window) algorithm outputs a score of each periodic sequence $X_{i}$. Then AUC can be calculated based on scores $\alpha _i$ for $X_{i}\ \forall i=1,\ldots ,s$.

An anomaly detector using the sliding-window size $\omega$ produces a total of $\hbox {length}(X)-\omega +1$ sequences from X. When calculating AUC, scores of the sliding sequences are transformed into periodic sequence scores as follows: Let $S_h$ be the anomaly score of sequence $X_h$, where $1 \le h \le (\hbox {length}(X)-\omega +1)$. The final score corresponds to a periodic sequence $X_i$ is the maximum score of $S_h\ \forall h$ such that at least half of $X_{h}$ is included in $X_{i}$.

1.2 Experiment settings for anomalous time series detection

[Datasets] Out of the initial 109 time series datasets (used in [11]) in UCR time series classification archive [10], we remove the datasets having the length of time series less than 200 and choose the top 20 datasets with the largest number of time series for our evaluation. (The datasets MixedShapesSmallTrain and FreezerSmallTrain are not included because they are the subsets of the two datasets chosen, namely, MixedShapesRegularTrain and FreezerRegularTrain, respectively).

New labeled datasets (normal vs. anomalous) are created as follows.

Each dataset in the archive contains a training set and a testing set. For each of the combined training-and-testing time series dataset having k classes, we take $\lfloor k/2 \rfloor$ largest classes and treat them to be the normality. Then, we sample 2% of time series from the other classes with an initial seed number 10 and treat them as anomalous time series.

The characteristics of the generated datasets are shown in Table 19. Note that the previous work [2, 3] chose only 1 class as the normality which makes the task trivial. The way we produce datasets allowing several classes as the normality makes this task more interesting and challenging.

Table 19 Characteristics of the 20 datasets used in the experiment of anomalous time series detection

[Algorithms] $\mathcal {K}_I$-based detectors are our implementations based on [59]. Mini-Rocket is from https://github.com/angus924/minirocket. TS2Vec is from https://github.com/yuezhihan/ts2vec. WTK is from https://github.com/BorgwardtLab/WTK. NeuTraL is from https://github.com/boschresearch/NeuTraL-AD. All the above are in Python. DOTS is a Scala implementation from https://github.com/B-Seif/anomaly-detection-time-series; Hyndman is a R implementation from https://github.com/robjhyndman/anomalous-acm. DTW and SBD are from the Python package tslearn [54].

The warping window of DTW is set to 5% of the length of the time series of each dataset, because the detection accuracy becomes worse when we attempt other window sizes or no windows.

For TS2Vec, we use its experimental settings of classification task because we require the instance-level representations.

For Mini-Rocket and NeuTraL, we use their default parameter configurations suggested in their papers.

The parameter search ranges of all the algorithms that need parameter tuning are given in Table 20.

Table 20 Parameter search ranges

As for the 1Line method, we use one of the following five types of basic vectorized primitive functions in Matlab as an anomaly score for each time series $T_j, j=1,\ldots ,n$, where n is the number of time series in the dataset.

• I. ±max($T_j$): the maximum value of time series

• II. ±min($T_j$): the minimum value of time series

• III. ±mean($T_j$): the mean value of time series

• IV. ±std($T_j$): the standard deviation of time series

• V. ±max(diff($T_j$)): the maximum of the difference between the current point and the previous point

We run these five 1-Line methods on each dataset and report the maximum AUC in Table 10. Low maximum values indicate that the datasets are hard to detect using any of the 1Line methods; otherwise, the datasets have anomalies which can be easily detected.

Table 21 Actual runtime (in CPU seconds) when the number of time series increases from $2^{10}$ to the total number of 50,000 on the InsectSound datatset

Table 22 Actual runtime (in CPU seconds) when the length of time series increases from $2^{14}$ to the full length (236,784) on the DucksAndGeese datatset

1.3 Additional experimental results for anomalous time series detection are given in Tables 21 and 22.