# Efficient pattern-based anomaly detection in a network of multivariate
devices

Len Feremans
Department of Computer Science
University of Antwerp
Belgium
<EMAIL_ADDRESS>
&Boris Cule
Department of Cognitive Science and Artificial Intelligence
Tilburg University
Netherlands
<EMAIL_ADDRESS>
&Bart Goethals
Department of Computer Science
University of Antwerp
Belgium
<EMAIL_ADDRESS>

###### Abstract

Many organisations manage service quality and monitor a large set devices and
servers where each entity is associated with telemetry or physical sensor data
series. Recently, various methods have been proposed to detect behavioural
anomalies, however existing approaches focus on multivariate time series and
ignore communication between entities. Moreover, we aim to support end-users
in not only in locating entities and sensors causing an anomaly at a certain
period, but also explain this decision. We propose a scalable approach to
detect anomalies using a two-step approach. First, we recover relations
between entities in the network, since relations are often dynamic in nature
and caused by an unknown underlying process. Next, we report anomalies based
on an embedding of sequential patterns. Pattern mining is efficient and
supports interpretation, i.e. patterns represent frequent occurring behaviour
in time series. We extend pattern mining to filter sequential patterns based
on frequency, temporal constraints and minimum description length. We collect
and release two public datasets for international broadcasting and X from an
Internet company. BAD achieves an overall F1-Score of 0.78 on 9 benchmark
datasets, significantly outperforming the best baseline by 3%. Additionally,
BAD is also an order-of-magnitude faster than state-of-the-art anomaly
detection methods.

_K_ eywords Anomaly detection $\cdot$ frequent pattern mining $\cdot$ pattern-
based embedding $\cdot$ timeseries similarity

##  1 Introduction

A large stream of information is produced by networks of connected devices in
the current Internet of Things world. Many companies are recording temporal
data streams. However, given a large network of devices it becomes infeasible
to inspect sensor logs manually and traditional expert-based supervision
methods are prone to false positives where errors are often detected late or
not at all [2]. Anomaly detection is an active research topic with
applications in event logs, time series, multivariate data and networks [23,
25, 28, 30]. We focus on _contextual_ anomaly detection in devices consisting
of multivariate time series where operators are interested in detecting if a
device behaves unexpectedly during a time period [4]. Since it is difficult to
collect labels we study this as an _unsupervised_ problem. Our research is
motivated by datasets collected from a company that monitors a _heterogeneous
network_ of different devices each consisting of several sensors which we make
publicly available. It occurs often that data is sent to multiple devices and
it is expected that because of load-balancing and broadcasting connected
entities exhibit anomalies at the same time.

Recent methods such as OmniAnomaly and Usad detect behavioural anomalies in a
single entity consisting of multiple univariate time series based on recurrent
neural networks and autoencoders [30, 2, 5]. However these methods ignores the
network context. Related work, that takes the network into account, is studied
specific to a single domain, such as social networks [28] or intrusion
detection [33].

We argue that including additional information for an entity in a network,
such as time series of connected entities, leads to an increase in anomaly
detection accuracy. For identifying intra-device similarities, we compare with
algorithms (and representations) that support fast similarity search in time
series databases. Piecewise Aggregate Approximation (Paa) and Product
Quantization reduce time series dimensions and scale to billions of time
series using GPU processing [19, 18]. For _heterogeneous times databases_ ,
_feature-based_ similarity measures have been proposed [17, 3] where time
series of different types are clustered based on multiple aspects such as the
distribution, peaks, entropy, change points and patterns, instead of computing
the similarity of point values using a distance metric. We focus on an
efficient algorithm using _feature-based_ similarity to tackle the high
variety in time series types. Our similarity score consider features based on
the overall trend, peaks and the distribution of a time series. After
construction the similarity matrix we could apply _domain knowledge_ to filter
spurious relations. However, domain knowledge related to sensor types is
generally not available. Hence, we reverse engineer domain knowledge thereby
computing patterns or _frequent relation types_ , to further spurious
relations. We also propose the usage of _density-based fingerprints_ , which
is a novel approximation of the _maximal information coefficient_ [21].

For anomaly detection _efficiency_ is important, especially for large networks
with many multivariate entities, where the training cost of learning a
separate neural network for each entity is prohibitive and impedes real-time
monitoring. Secondly, _interpretability_ is crucial and we want to assist
operators in troubleshooting the root-cause and not only to enable them to
locate the device and sensor causing an anomaly during a certain time period.
For anomaly detection in data series we compute a high anomaly score if _few
recurring patterns occur_. Fpof, is one of the first _pattern-based anomaly
detection methods_ [14] and recently extended by MivPod [15] and Pbad [8]. All
three methods apply a sliding window, discover frequent patterns and compute
an outlier or anomaly score based on pattern occurrences in each window.
However, there are major differences in: (1) the type of patterns being mined,
(2) how to select an interesting pattern, (3) how to compute the anomaly
score. Existing pattern-based anomaly detection algorithm focus on frequent
itemsets [14], maximal sequential patterns [8] or minimal infrequent itemsets
[15]. Different to existing methods we mine frequent sequential patterns
subject to both _temporal_ and _information-theoretical_ constraints. We use a
constraint on _minimum length_ to avoid short sequential patterns and on
_relative duration_ , or cohesion, to ensure that any occurrence of a
candidate sequential pattern $X$ contains few gaps [7]. For instance, both
$(a,b,c)$ and $(d,e,f)$ occur in sequence $(a,b,c,x,d,x,x,e,x,f)$, but we only
consider the occurrence of pattern $(a,b,c)$. Additionally, this constraint
results in a smaller search space and fewer patterns. A second problem with
existing pattern-based anomaly detection methods, is that an anomaly
prediction is not always explainable. While a pattern and its occurrences can
be inspected if an anomaly is reported using the previously mentioned pattern-
based anomaly detection algorithms, it is arguably _not explainable_ if there
are thousands of patterns. Hence, we propose to filter frequent sequential
patterns using the _Minimum Description Length (MDL)_ principle [12]. That is,
we compute how _many bits are saved_ by storing the sequence database using
patterns including the description of the patterns themselves [32]. An
overview of our context-aware Pattern-Based Anomaly Detection (Bad) method is
shown in Figure 1.

The problem we are tackling is defined as follows: Given a network of
multivariate time series (or devices) identify periods of abnormal behaviour.
Our main contributions are:

  * •

A novel pattern-based anomaly detection method. We search for patterns that
are frequent and filter interesting patterns based on temporal and
information-theoretical constraints. Since, the set of patterns is relatively
small and sparse we facilitate human oversight.

  * •

A novel method to reverse engineering a graph of connected devices. We
represent each time series using a fingerprint such that comparing pairwise
similarities is both fast and accurate. Thereafter, we prune spurious
relations based on patterns on the type of sensor of common occurring
relations.

  * •

For multivariate time series, we construct a pattern-based embedding for
univariate time series and report anomalies by considering the joint
probability of patterns occurring in each time series using an isolation
forest or Fpof for interpretability [25, 14].

  * •

Experimentally we demonstrate the interpretability, efficiency and high
anomaly detection performance of the proposed method on large-scale
heterogeneous networks.

Figure 1: Overview of the proposed Bad method: 1) we search for connected time
series, e.g. sensor $S_{1}$ of device $D_{1}$ is similar to $S_{3}$ of device
$D_{2}$; 2) we extend $D_{1}$ and include time series $S_{3}$ from $D_{2}$ as
an additional signal; 3) we transform each time series to a discrete sequence
using Sax; 4) we search interesting sequential patterns; 5) we assign an
anomaly score based on pattern occurrences, e.g $P_{1}=adda$ is absent in the
3th window in $E_{1}$ and $E_{3}$

The remainder of this paper is organised as follows. In Section 2, we
introduce the necessary preliminaries. In Section 3, we present a novel
pattern-based anomaly detection method. In Section 4, we present a method for
identifying connected devices based on sensor time series. In Section 5, we
combine both steps for detecting anomalies in multivariate time series. In
Section 6, we present an experimental evaluation of our method and compare it
with state-of-the-art methods. We discuss the related work in Section 7 and
conclude our work in Section 8.

##  2 Preliminaries

###  2.1 Time series data

A _time series_ , or sensor, is defined as a sequence of measurements
$S=(x_{1},x_{2},\ldots,x_{n})$ where $x_{i}\in\mathbb{R}$ is sampled at
regular time intervals. We denote the type of a sensor using
$\mathit{type}(S)$.

A time series _window_ , or interval, is a contiguous subsequence of a time
series $S_{i:j}=(x_{i},x_{i+1},\ldots x_{j})$ where $1\leq i<j\leq n$.

Given a window length $l$ and increment $t$ we can transform each time series
to a set of _sliding windows_ , i.e., $S_{1:l},S_{1+t:l+t},\ldots,S_{n-l:n}$
into a continuous sequence database $D$ is, that is

$D=\\{S_{i:j}\ |\ S_{i:j}\in\mathit{sliding\\_window}(S,l,t)\\}$

The size of the sequence database is denoted by $|D|=(|S|-l+1)/t$ where $l$ is
the window length and $t$ an increment of the sliding window transform. In
practice, we select $t=1$.

An _entity_ , or device, consists of $m$ time series, i.e.,
$D=\\{S^{1},S^{2},\ldots,$ $S^{m}\\}$ where all time series start and stop at
the same time. We denote the type of a device using $\mathit{type}(D)$. We
assume that each device of the same type has the same amount and types of
sensors.

A multivariate window consists of $m$ subsequences (or dimensions) denoted
$D_{i:j}=\\{S^{1}_{i:j},S^{2}_{i:j},\ldots,S^{m}_{i:j}\\}$.

A _network_ , $G=(V,E)$ consists of $k$ entities $V(G)=\\{D_{1},D_{2},\ldots,$
$D_{k}\\}$ and a set of edges between entities $E(G)$. We assume that the set
of edges is unknown as edges are often dynamic in nature and caused by an
unknown underlying data generating process.

###  2.2 Sequential pattern mining

We use _frequent pattern mining_ for discovering recurring behaviour. For
pattern mining we requires a database of discrete sequences. Therefore, we
transform each time series into a discrete sequence using an extension of
Symbolic Aggregate Approximation (Sax) [24].

A _sequence database_ $\mathcal{D}$ is created by applying a sliding window on
a time series $S$ and transforming each continuous window $S_{i:j}$ to a
discrete sequence $\hat{S}_{i:j}$ using Sax, where

$\mathcal{D}=\\{\hat{S}_{i:j}\ |\ S_{i:j}\in
D:\hat{S}_{i:j}=\textsc{Sax}(S_{i:j},\mathit{paa\\_win},\mathit{no\\_bins})\\}.$

Given a sequence database we can discover frequent unordered or ordered
patterns. Many algorithms have been proposed to efficiently mine (non-
redundant) frequent patterns and are publicly available in toolkits such as
Spmf [34, 10].

A _sequential pattern_ $X=(e_{1},e_{2},\ldots e_{p})$ is an ordered list of
one or more discretised values or items. We use $|X|$ to denote the length of
a pattern.

A sequential pattern $X$ _occurs_ in a discretised sequence, i.e.
$X\prec\hat{S}_{i:j}$, if it is a subsequence of $\hat{S}_{i:j}$ thereby
allowing gaps. The _cover_ of a sequential pattern $X$ in $\mathcal{D}$ is
defined as

$\mathit{cover}(X,\mathcal{D})=\\{\hat{S}_{i:j}\ |\
\hat{S}_{i:j}\in\mathcal{D}\wedge X\prec\hat{S}_{i:j}\\}.$

For instance, $X=(a,d,d,a)$ occurs in
$\hat{S}=(c,c,\underline{a},\underline{d},\underline{d},c,\underline{a})$.

The _support_ of a pattern is defined as
$\mathit{support}(X,\mathcal{D})=|\mathit{cover}(X,\mathcal{D})|$. The
relative support is defined as
$\mathit{rsupport}(X,\mathcal{D})=|\mathit{cover}(X,\mathcal{D})|/|\mathcal{D}|$.
A sequential pattern is _frequent_ if its support is higher than the user-
defined minimal support threshold.

More recently, authors suggest to mine the _top- $k$ most frequent patterns_
directly, which is more intuitive for end-users since choosing minsup depends
on a database characteristics which are unknown [9, 20].

An issue with baseline sequential pattern mining algorithms is that they
discover many _redundant_ and uninteresting patterns. We impose additional
_temporal constraints_ based on minimal length and the number of gaps to find
more interesting patterns in data series [6].

A constraint on the _minimal length_ defines that for any discovered frequent
sequential pattern $X$, the pattern length should be higher than a user-
defined threshold, i.e. $|X|>\mathit{min\\_len}$ [1]. By default we set
$\mathit{min\\_len}=3$ to filter sequential patterns consisting of 1 or two
items, since they only cover a small fraction of a sequence and are often
frequent due to chance.

A constraint on the _relative duration_ determines the cohesion of sequential
pattern occurrences [7]. We impose a constraint on the duration relative to
the length of pattern proposed, i.e.

$\mathit{duration}(X,\hat{S}_{i:j})/|X|\geq\mathit{min\\_relative\\_duration}.$

For instance, given a sequence $(a,b,c,x,d,x,x,e,x,f)$ the relative duration
of sequential pattern $(a,b,c)$ is $\frac{3}{3}$, while for $(d,e,f)$ it is
$\frac{6}{3}$. With a relative duration of 1.0 we only count sequential
patterns that occur consecutively in each sequence window. By default we set
minimal relative duration to 1.2. We remark that by incorporating gap
constraints during the mining process, we reduce the search space [26, 7].

###  2.3 Symbolic Aggregate Approximation

Figure 2: Illustration of the Sax transform where a time series (or time
series window) with length $128$ is transformed into the discrete sequence
$(b,a,a,b,c,c,b,c)$ of length $8$ using $\mathit{no\\_bins}=$ and
$\mathit{paa\\_win}=16$.

We tranform each continuous window $S_{i:j}$ to a discrete sequence
$\hat{S}_{i:j}$ using Sax and segment a time series into $n$ pieces, down-
sample using averaging and then discretise into $k$ bins. An example of Sax is
shown in Figure 2. We use non-overlapping sliding windows to compute average
values. This is controlled using parameter $\mathit{paa\\_win}$, which
determines how many raw values are averaged, e.g., in the example
$\mathit{paa\\_win}$ is $16$. Parameter $\mathit{no\\_bins}$ controls how many
symbols or bins are used for discretisation, e.g., in the example
$\mathit{no\\_bins}$ is $3$.

####  2.3.1 Binning methods

For binning several options are possible, i.e., equal-length or equal-width
bins or bins selected assuming a normal distribution. The default binning
method used in our experiments is normal where we use the same equal-length
bins for all windows in one time series, i.e., we compute $\mathit{no\\_bins}$
equal-length bins between the lowest and highest value in the entire time
series. The second strategy is local, where we use local bins for each time
series window, i.e., we create bins $[a,b,c]$ (or low, medium, high) relative
to the minimal and maximal value in one window. This is the strategy
originally proposed by Sax [24]. As a third alternative, we use k-means
clustering to bin all time series values into $k$ clusters, where $k$ is
chosen to be equal to $\mathit{no\\_bins}$.

Each binning approach has its advantages and disadvantages, i.e., using local
binning we match the same shape, such as an increase in slope, irregardless of
the absolute values. The main difference with the normal strategy is that we
compare relative values of symbols and not absolute values. However, using the
normal strategy a certain pattern $X=(a,b,c)$ meaning low, medium, high value
(assuming 3 bins) is matched in the entire time series in the same manner and
therefore easier to interpret and visualise. Finally, k-means, arguably, leads
to a superior binning strategy, compared to equal-length bins. However, we
face the problem that clustering is stochastic and can produce different bins
in highly similar time series.

##  3 Pattern-based anomaly detection

###  3.1 Definition of interesting patterns

Given a time series $S$ is transformed into a discrete sequence database
$\mathcal{D}$ using a sliding window and discretised using Sax. We discover a
set of _interesting_ sequential patterns $X\in\mathcal{P}$, defines as:

  1. 1.

$|X|>\mathit{min\\_len}$

  2. 2.

$\mathit{bits\\_saved}(X,\mathcal{D})>0$

  3. 3.

$X$ is ranked in the top-$k$ most frequent sequential patterns

  4. 4.

The frequency of the pattern $X$ is computed w.r.t. relative duration, i.e.
$rank(X)<rank(Y)\leftrightarrow|\mathit{cover_{rdur}}(X,\mathcal{D})|>|\mathit{cover_{rdur}}(Y,\mathcal{D})|$

####  3.1.1 Temporal constraints

The cover of a pattern $X$ w.r.t. relative duration is defined using:

$\displaystyle\mathit{cover_{rdur}}(X,\mathcal{D})=\\{\hat{S}_{i:j}\ |$
$\displaystyle\ \hat{S}_{i:j}\in\mathcal{D}:X\prec\hat{S}_{i:j}\ \wedge$
$\displaystyle\frac{\mathit{duration}(X,\hat{S}_{i:j})}{|X|}\geq\mathit{rdur}\\}.$

We remark that the maximal length of $X$ is constrained by the Sax
representation, i.e. given a sliding window $l$ and $\mathit{paa\\_win}$,
$|\hat{S}_{i:j}|=l/\mathit{paa\\_win}$.

####  3.1.2 Minimal description length

We are inspired by _MDL_ as an additional criterion to individually select
sequential patterns. That is, we keep a sequential pattern if its frequent,
conforms to temporal constraints, and compresses the data series, i.e., if we
need fewer symbols to encode the data series using the pattern (including the
pattern itself) than we need to encode the data series literally. The
resulting set of patterns $\mathcal{P}$ should describe all regularities in
the time series. Formally, given a sequential pattern $X\in\mathcal{P}$ and
its cover $\mathcal{D}_{X}=\mathit{cover}({X,\mathcal{D})}$, we define the
number of bits saved using:

$\begin{gathered}\mathit{bits\\_saved}(X,\mathcal{D})=\sum_{\hat{S}_{i:j}\in\mathcal{D}_{X}}{\mathit{DL}(\hat{S}_{i:j})}-(|X|\cdot
log(n)+\sum_{\hat{S}_{i:j}\in\mathcal{D}_{X}}{\mathit{DL}(\hat{S}_{i:j}|\
X)}).\end{gathered}$

Here $\mathit{DL}(\hat{S}_{i:j})$ is the description length to store the
covered sequences sequence computed using Huffman coding111We remark that
Huffman coding computes an optimal table of varying-length codes for all
sequences in $\mathcal{D}_{X}$. For computing description length we
concatenate all covered sequences and ignore the storage cost of the Huffman
table.. The term $|X|\cdot log(n)$ is the number of bits used to store the
pattern, and $\mathit{DL}(\hat{S}_{i:j}|X)$ is computed by replacing the
pattern occurrence of $X$ in $\hat{S}_{i:j}$ and then computing the
description length of the remaining symbols. For example, assuming $X=(a,b,c)$
and $\hat{S}_{i:j}=(x,y,a,b,c,z)$ we compute the description length of
$\mathit{DL}(\hat{S}_{i:j}|X)=\mathit{DL}((x,y,*,z))$ where $*$ is an
additional symbol to demarcate the location(s) of pattern $X$. In case the
occurrence contains gaps we keep gap symbols. We remark that the number of
bits saved is maximal for patterns that are both long and frequent, while
frequency is always maximal for short patterns. In preliminary experiments, we
found that by filtering patterns where
$\mathit{bits\\_saved}(X,\mathcal{D})>0$, we have an order-of-magnitude less
patterns.

###  3.2 Efficient algorithm

Input : Time series $S$, sliding window interval ($l$), Sax parameters
($\mathit{paa\\_win}$ and $no\\_bins$), number of patterns ($k$), min length
pattern (_min_len_) and min relative duration (_rdur_)

Result: Interesting sequential patterns

1
$\mathcal{D}\leftarrow\\{\hat{S}_{i:j}|S_{i:j}\in\mathit{sliding\\_window(S,l,1)}:\hat{S}_{i:j}=\textsc{Sax}(S_{i:j},\mathit{paa\\_win},\mathit{no\\_bins})\\}$

$\mathit{stack}\leftarrow\left[\langle\emptyset,\mathcal{D},\Omega\rangle\right]$;

// Contains current pattern, projection, candidate suffix items

$\mathcal{P}\leftarrow\mathit{make\\_heap}(k)$;

// Max $k$ patterns sorted on support

2 $\mathcal{P}_{\mathit{final}}\leftarrow\emptyset$

3 while  _$\mathit{stack}\neq\emptyset$_ do

4 $\langle X,\mathcal{D}_{X},Y\rangle\leftarrow\mathit{stack}.\mathit{pop}()$

5 if  _$Y=\emptyset$_ then // Stop condition: check if pattern is interesting

6

7 if  _$\mathit{support}_{\mathit{rdur}}(X,\mathcal{D})
>\mathit{min\\_heap}(\mathcal{P})\ \mathbf{and}\ |X|\geq\mathit{min\\_len}$ _
then

8 if  _$\mathit{bits\\_saved}(X,\mathcal{D}) >0$_ then

9
$\mathcal{P}_{\mathit{final}}\leftarrow\mathcal{P}_{\mathit{final}}\cup\\{X\\}$

10

11 $\mathcal{P}.\mathit{push\\_pop}(X)$

12

13 else // Branch and bound: add super-pattern

14 if  _$\mathit{support}_{\mathit{rdur}}(X,\mathcal{D})
<\mathit{min\\_heap}(\mathcal{P})$_ then

15 $\mathbf{continue}$

16

17 $s_{k+1}\leftarrow Y[0]$

18 $X_{k+1}\leftarrow X\oplus s_{k+1}$

19 $\textit{stack}.\mathit{push}(\langle
X,\mathcal{D}_{X},Y\setminus\\{s_{k+1}\\}\rangle)$

20
$\mathcal{D}_{X_{k+1}}\leftarrow\textsc{Project}(\mathcal{D}_{X},X_{k+1},\mathit{rdur})$

21 $Y_{X_{k+1}}\leftarrow\textsc{Candidates}(\mathcal{D}_{X_{k+1}})$

22 $\mathit{stack}.\mathit{push}(\langle
X_{k+1},\mathcal{D}_{X_{k+1}},Y_{X_{k+1}}\rangle)$

23

$\mathbf{return}\ \mathcal{P}_{\mathit{final}}$

Algorithm 1 Mip($S$, $l$, $\mathit{paa\\_win}$, $no\\_bins$, $k$,
$\mathit{min\\_len}$, $\mathit{rdur}$) Discover sequential patterns that are
frequent, satisfy both temporal and MDL constraints

We now present and algorithm to mine interesting patterns. The algorithm is an
extension of the algorithm proposed by [7] for mining sequential patterns with
respect to temporal constraints. Major differences are that we compute
multiple pattern-occurrences for each window add select patterns using MDL.

Mine_Interesting_Patterns (Mip), is shown in Algorithm 1. We start by creating
an empty stack. During depth-first recursion we maintain candidate sequential
pattern prefixes $X$ and a set of candidate items $Y$. At each iteration we do
a divide-and-conquer step, i.e. assume $X=(s_{1},\ldots,s_{k})$ and
$Y=\\{s_{a},s_{b},s_{c}\\}$ we generate candidate supersequences $X_{k+1}$ by
adding every item $s_{k}\in Y$, i.e. $X_{k+1}=X\oplus
s_{k+1}=(s_{1},\ldots,s_{k},s_{k+1})$. If no more candidate items exists we
evaluate $X$ and check the frequency, minimum length and MDL constraints. The
space of possible varying-length sequential patterns is exponential, thus
pruning is extremely important. We limit candidates using three different
ways.

  * •

First, if $X$ is infrequent, we know all super-sequences are infrequent and
stop early. We use the minimal support in the heap to prune $X$ instead of a
fixed value for $\mathit{min\\_sup}$ (if fewer than $k$ candidates are in the
heap, this is 0). As the number of iterations increase, so does the minimal
support of the worst pattern 222Since the minimal support threshold
monotonically increases, we trivially guarantee that the final heap contains
indeed the top-$k$ most frequent patterns..

  * •

Second, we use prefix-projected pattern-growth [13]. Hence, we make use of an
additional data-structure, i.e. the _prefix-projection_ , which we compute for
each candidate during recursion. The projection of the database $\mathcal{D}$
on $X$ contains all sequences covered by $X$ and a pointer to the first and
last occurring item $s_{k}$. In case $X$ occurs multiple times in a sequence,
we keep a reference to all positions. We update the prefix-projected database
for each new super-sequential pattern, instead of having to scan the entire
database, which is more efficient, especially for longer (or low-support)
patterns where $|\mathcal{D}_{X}|\ll|\mathcal{D}|$. Additionally, we compute
the set of candidates $Y$ to make super-sequential patterns based on the
projection.

  * •

Third, we further limit the number of candidates based on relative duration.
If the minimal relative duration is 1, the set of candidates consists of the
next item after $s_{k}$ in each covered sequence in $\mathcal{D}_{X}$. If the
relative duration is higher we have an upper bound of
$\lfloor\mathit{max\\_len}\times duration\rfloor-\mathit{max\\_len}$ on the
number of gaps and consequently limit the search to a only a few items after
$s_{k}$ in $\mathcal{D}_{X}$.

####  3.2.1 Complexity

The algorithm is memory-friendly, i.e. given the depth-first search strategy
of searching for frequent patterns, we maintain at most $\mathit{max\\_len}$
candidates. Secondly, we do not actually materialise (or copy) the original
sequence database in memory to compute the prefix-projection of the database,
but rather compute a _pseudo_ prefix-projection thereby storing only indices
(or pointers) to covered subsequences and the first and last item position.

The algorithm generates fewer candidates then typically in sequential pattern
mining. To enumerate all sequential patterns of size $3$ with $n$ symbols
without pruning would generate $n^{3}$ candidates333We remark that repeating
the same symbol is allowed which is not always the case in sequential pattern
mining literature. Traditional level-wise or column-wise algorithms, would
first generate $n$ singletons. Assuming $n^{\prime}$ singletons are frequent,
we would then generate $n^{\prime}\times n$ candidates of size two. Next,
assuming $n^{\prime\prime}$ frequent patterns of size 2, we would then
generate $n^{\prime\prime}\times n$ candidates of length 3. In contrast, we do
not consider $n$ super-sequences for each candidate pattern of a certain size,
but generate candidates based on items that occur at least once (after the
last occurring item of $X$) based on the prefix-projected database.

It is also interesting to compare the complexity with consecutive sequence,
motif, or _string mining_ algorithms [27]. To generate all frequent strings of
length 3 in a discrete sequence of length $k$ we have only $k-3+1$ candidates,
i.e. we would traverse the sequence from left the right using a sliding window
of size 3 and an increment of 1. To generate all frequent string of length 3
with at most 1 gap, we have $k-4+1$ windows and need to check $\binom{4}{3}=4$
combinations for each window, or in general, for a pattern of size $l$ with
$g$ gaps, we have $(k-(l+g)-1)\times\binom{l+g}{l}$ candidates which is
inefficient. We limit the number of candidates in each suffix based on
relative duration and leverage the anti-monotonicity of support, thereby
combining optimisations from both string and sequential pattern mining.

###  3.3 Create an embedding

Given a set of patterns $\mathcal{P}=\\{X_{1},X_{2}\ldots,X_{k}\\}$ we
construct an embedding $E_{i:j}\in\mathbb{R}^{k}$ for each discretised window
$\hat{S}_{i:j}$, using

$E_{i:j}=[f(X_{1},\hat{S}_{i:j}),\ldots,f(X_{k},\hat{S}_{i:j})]\quad\text{where}$
$f(X_{n},\hat{S}_{i:j})=\big{\\{}\begin{array}[]{ c l
}\mathit{rsupport}(X_{n},\mathcal{D})&\quad\textrm{if
}X_{n}\prec\hat{S}_{i:j}\\\ 0&\quad\textrm{otherwise}\end{array},$

where we use the _relative support_ of a pattern as a feature as opposed to a
simpler binary function i.e., 1 if $X_{k}\prec\hat{S}_{i:j}$ and 0
otherwise444We remark that a a distance-weighted similarity score is applied
by [8]..

A naïve algorithm for creating the embedding would search for occurrences for
each pattern $X$ in all sequences in $\mathcal{D}$ and has a complexity of
$O(|\mathcal{D}|\times|\mathcal{P}|)$. In preliminary experiments we found
that the naïve algorithm has a longer runtime than pattern mining itself, so
its optimisation is also crucial. We adopt the prefix-projected data-
structure, but for _computing occurrences_ of patterns in any data series. For
a pattern $X=(s_{1},\ldots,s_{k})$ we compute the projection for each prefix
$(s_{1}),(s_{1},s_{2}),\ldots$, $(s_{1},\ldots,s_{k})$ incrementally. The
number of sequences matching in each successive projection grows smaller
(since support is anti-monotonic), thereby _early abandoning_ a sequence when
if fails to match a prefix pattern. This method creates the embedding more
efficiently by leveraging the _sparsity_ of the embedding and temporal
constraints.

###  3.4 Compute the anomaly score

####  3.4.1 Isolation forest

Given the pattern-based embedding for each window
$E_{i,j}\in\mathbb{R}^{|\mathcal{P}|}$ we compute anomaly scores using an
_isolation forest_ [25]. An embedding $E_{i:j}$ has a low anomaly score if
many frequent patterns occur and the (average) depth of random trees is
relatively high, or conversely a high anomaly score if few frequent patterns
occur in the embedding vector for each window.

####  3.4.2 Frequent pattern outlier factor

To improve interpretability we compute anomalies alternatively using the
_frequent pattern outlier factor_ (Fpof) [14]. The anomaly score in for window
$\hat{S}_{i:j}$ is computed using

$\mathit{fpof}(\hat{S}_{i:j})=\frac{1}{|\mathcal{P}|}\sum_{X\in\mathcal{P}\wedge
X\prec\hat{S}_{i:j}}\mathit{rsupport}(X,\mathcal{D}).$

Here, the outlier factor is high (close to 1) if many frequent patterns
(weighted by the relative support) match the current window and low otherwise.
For ranking anomaly scores we use $1-\mathit{fpof}(\hat{S}_{i,j})$.

##  4 An efficient algorithm to identify inter-device similarities

In the previous section we explained a method to detect anomalies based on
interesting patterns. A second goal is to determine edges between different
devices based on similar time series. We assume a heterogeneous set of $k$
entities $D_{1},D_{2},\ldots,D_{k}$ where each entity is of certain type, e.g.
a server for computing or storage, an industrial device, a client-oriented
device, etc. Moreover, each entity consist of $m$ sensors or time series, each
of a certain type, e.g. CPU or RAM or disk usage, in going or outgoing network
traffic etc. Our aim is to discover relations, or edges, between device
sensors, i.e. two devices having a similar CPU load or network transmission.
This is challenging because of following reasons.

  * •

Give a large network of devices and sensors, let $N$ denote the total number
of time series and $|S|$ the length of time series. It is not scalable to
compare each pair of time series, since the runtime is $O(N^{2}\cdot|S|)$.

  * •

There is _domain knowledge_ related to type , i.e. time series of the same
type are often similar to each other, e.g. “device-A/CPU load”. Vice versa, if
only single pair of time series “CPU temperature” is similar to “network bit
rate” this is probably spurious.

  * •

Similar to [17] we find that common point-based similarity measures, such as
Euclidean distance, Pearson correlation or Dynamic Time Warping often result
in _false positives_ in a heterogeneous time series databases. Which is
aggravated since time series often exhibit a similar trend, being mostly flat,
monotonically increasing, or exhibiting periodical concept drift corresponding
to day and night activity.

  * •

[21] report that common similarity scores, except mutual information, are not
_equitable_. Since, we have heterogeneous sensor types (with high variety in
both distribution and trend) this means that a single threshold on the
similarity score, will result in a skew in the degree of connections based on
the type of sensors.

###  4.1 Density-based fingerprint

First, we generate _density-based fingerprint_. We create a grid of $n\times
k$ cells for each time series thereby computing the _density_ , i.e., the
number of samples $x_{i}$ in each cell. Additionally, we use _logarithmic
scaling_ to give more weight to low-density cells such as _peaks_. Formally,
the distance between two time series $S^{i},S^{j}$ using fingerprinting is
defined as

$\begin{gathered}\mathit{dist\\_fp}(S^{i},S^{j},n,k)=\\\
\sum_{x=1}^{n}\sum_{y=1}^{k}\|\mathit{log}(I_{x,y}(S^{i})+1)-\mathit{log}(I_{x,y}(S^{j})+1)\|_{2}\end{gathered}$
(1)

where $I_{x,y}(S)$ denotes the number of sensor values in the cell $(x,y)$.
Note that in the limit, i.e., for high values of $n$ and $k$, the fingerprint
distance converges towards the Euclidean distance. Because we are dealing with
time series, we usually set $n$ using an _interval_ and discretise the sensor
values into a small number of _bins_. Fingerprints are inspired by _the
maximal information coefficient_ , which is a similarity metric defined using
the mutual information on a discretised grid of values which is _equitable_ ,
i.e. the similarity score is comparable between different types of time series
[21].

Figure 3: Example of two continuous time series collected over 98 days.
Fingerprints are created by counting values in each cell using an interval of
$24$ hours and $4$ bins resulting in $103\times 4$ values instead of $28\ 224$
values

An example is shown in Figure 3, where we consider a fingerprint as a low
dimensional visual summary for each time series. We highlighted the _concept
drift_ in green which foils common point-based distance measures, such as
Pearson correlation and Dynamic Time Warping. In this example, the distance
between the two fingerprints is low as there is a large correspondence in
_peaks_ , i.e., light grey cells, and the overall _trend_ , i.e., black and
white cells.

###  4.2 Histogram distances

The histogram is invariant to phase changes and complementary to distance
measures in the time domain. We use histogram distances as a second feature to
avoid false positives. Given time series $S^{i},S^{j}$ we define the distance
between histograms as

$\mathit{dist\\_hist}(S^{i},S^{j},l)=\sum_{k=1}^{l}|H^{k}(S^{i})-H^{k}(S^{j})|$

where $l$ is the number of equal-length bins and $H^{k}(S)$ denotes the number
of sensor values in bin $k$ divided by $|S|$.

###  4.3 Feature-based distance

We define two time series are similar, i.e., $S^{i}\sim S^{j}$, if:

  1. 1.

$\mathit{dist\\_fp}(S^{i},S^{j},n,k)<t_{f}$

  2. 2.

$\mathit{dist\\_hist}(S^{i},S^{j},l)<t_{h}$

  3. 3.

$(\mathit{type}(S^{i}),\mathit{type}(S^{j}))\in\mathit{frequent\\_relations}$

where $t_{f}$ is a threshold for fingerprint-based distance and $t_{h}$ is a
threshold for histogram-based distance. First we compute the distance using
fingerprinting and histograms between all $\frac{N\times(N-1)}{2}$ pairs to
create a symmetric similarity matrix $\mathbf{M}\in\mathbb{R}^{N\times N}$.
Next, we search $\mathbf{M}$ to compute the frequency of a relation between
two types of sensors as discussed next.

Input : A list of time series filenames $\mathbf{S}$; interval, number of bins
and threshold for fingerprint-based and histogram-based distances; coverage
for finding frequent relation types

Result: Set of edges between devices with similar time series

// 1\. Create feature-based representation

1 $\mathbf{F}\leftarrow\emptyset$; $\mathbf{H}\leftarrow\emptyset$

2 foreach  _$S^{i}\in\mathbf{S}$_ do

3 $S^{i}\leftarrow\mathit{min\\_max\\_norm}(\mathit{load\\_ts}(S^{i}))$

4
$\mathbf{F}_{i}\leftarrow\mathit{create\\_fingerprint}(S^{i},\mathit{interval},\mathit{bins_{f}})$

5
$\mathbf{H}_{i}\leftarrow\mathit{create\\_histogram}(S^{i},\mathit{bins_{h}})$

6

// 2\. Compute similarity matrix

7 $\mathbf{M}\leftarrow 0^{|\mathbf{F}|\times|\mathbf{F}|}$

8 for  _$i\leftarrow 1\ \textnormal{{to}}\ |\mathbf{F}|$_ do

9 for  _$j\leftarrow i+1\ \textnormal{{to}}\ |\mathbf{F}|$_ do

10 $\mathit{d\\_fingerprint}\leftarrow\|\mathbf{F}_{i}-\mathbf{F}_{j}\|_{2}$

11 $\mathit{d\\_hist}\leftarrow\|\mathbf{H}_{i}-\mathbf{H}_{j}\|_{1}$

12 if  _$\mathit{d\\_fingerprint}
<\mathit{t_{f}}\wedge\mathit{d\\_hist}<\mathit{t_{h}}$_ then

13 $\mathbf{M}_{i,j}\leftarrow\mathit{d\\_fingerprint}$

// 3\. Find frequent relation types

14$\mathit{freq\\_relations}\leftarrow\textsc{FindFrequentRelationTypes}(\mathbf{M},\mathit{coverage})$

// 4\. Create network of connected devices

15 $E(G)\leftarrow\\{(D_{i},D_{j})|\ \exists S^{i}\in D_{i},S^{j}\in
D_{j}:\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\mathbf{M}_{i,j}>0\
\wedge(\mathit{type(S^{i})},\mathit{type(S^{j})})\in\mathit{freq\\_relations}\\}$

$\mathbf{return}\ E(G)$

Algorithm 2 BuildNetwork($\mathbf{S}$, $\mathit{interval}$,
$\mathit{bins_{f}}$, $\mathit{t_{f}}$, $\mathit{bins_{h}}$, $\mathit{t_{h}}$,
$\mathit{coverage}$) Find connected devices in a network

###  4.4 Find frequent relation types

We propose to discover _frequent relation types_. Thereby, we make the
assumption that, if the data series between sensors of a certain type are
often similar, or frequent, this is a _pattern_ in the domain. We remark that
the resulting list of frequent relation types also has value for the end-user,
i.e. it provides a _comprehensible summary_ of the edges in a heterogeneous
network.

We represent the similarity matrix $\mathbf{M}\in\mathbb{R}^{N\times N}$ by an
_attributed network_ , i.e., $G(V,E)$ where each nodes represents a sensor
instance labelled by its type and $E(G)=\\{(S^{k},S^{l})\ |\ S^{k}\sim
S^{l}\\}$. Formally, we define the _support_ of different types of edges
between devices as:

$\begin{gathered}\mathit{support_{RT}}(G,\mathit{type_{i}},\mathit{type}_{j})=|\\{(D_{x},D_{y})\
|\exists S^{k}\in D_{x},S^{l}\in D_{y}:\\\ (S^{k},S^{l})\in
E(G)\wedge\mathit{type(S^{k})}=\mathit{type}_{i}\wedge\mathit{type(S^{l})}=\mathit{type}_{j}\\}|.\end{gathered}$

The algorithm for discovering frequent relation types is inspired on
_separate-and-conquer rule learning_ [11]. First, we identify the relation
type between sensors with the highest support. Next, we remove all devices
connected through this relation. In the remaining database, we identify the
next most frequent relation type and continue recursively, until we cover, for
instance, $90\%$ of time series.

###  4.5 Algorithm

Input : A similarity matrix $\mathbf{M}$ between time series;
$\mathit{coverage}$ (%) parameter, i.e., required percentage of relations in
$\mathbf{M}$ to be covered

Result: Set of frequent relation types

// 1\. Create database

1 $\mathbf{T}\leftarrow\emptyset$

2 foreach  _$(i,j)\in\mathbf{M}$_ do

3 if  _$\mathbf{M}_{i,j}
>0\wedge\mathit{device}(S^{i})\neq\mathit{device}(S^{j})$_ then

4
$\mathbf{T}\leftarrow\mathbf{T}\cup(\mathit{device}(S^{i}),\mathit{device}(S^{j}),\mathit{type}(S^{i}),\mathit{type}(S^{j}))$

// 2\. Separate-and-conquer

5 $\mathit{freq\\_relations}\leftarrow\emptyset$

6 $\mathit{min\\_cover}\leftarrow|\mathbf{T}|*(1-\mathit{coverage})$

7 while  _$|\mathbf{T}| >\mathit{min\\_cover}$_ do

8
$\mathit{freq\\_rel}\leftarrow\underset{(t_{1},t_{2})\in\pi_{\mathit{type}_{1},\mathit{type}_{2}}(\mathbf{T})}{\mathit{argmax}}\mathit{support_{RT}}(T,t_{1},t_{2})$

9
$\mathit{freq\\_relations}\leftarrow\mathit{freq\\_relations}\cup\mathit{freq\\_rel}$

10
$\mathit{device\\_pairs}\leftarrow\pi_{dev_{1},dev_{2}}(\sigma_{\mathit{type}_{1}=t_{1}\wedge\mathit{type}_{2}=t_{2}}(\mathbf{T}))$

11 $\mathbf{T}\leftarrow\mathbf{T}\ \setminus\ \\{r\ |\
r\in\mathbf{T}:(r.dev1,r.dev2)\in\mathit{device\\_pairs}\\}$

$\mathbf{return}\text{ }\mathit{freq\\_relations}$

Algorithm 3 FindFrequentRelationTypes($\mathbf{M}$, $\mathit{coverage}$) Find
frequent relation types using separate-and-conquer

The main BuildNetwork algorithm is shown in Algorithm 2. We load time series
one by one, pre-process them and compute a memory-friendly representation
consisting of one fingerprint and histogram for each time series (line 2-2).
Next, we compute distances between all $N^{2}/2$ candidates and create a
similarity matrix $\mathbf{M}$ (line 2-2).

In Algorithm 3 we show the algorithm for discovering frequent relation types
based on _coverage_. First, we join all pairs of discovered similar time
series with metadata, i.e., the device and type of both time series. Next, we
identify the relation type between time series with the highest support (line
3). Next, we identify all pairs of devices connected through this frequent
relation type and remove all edges between these devices. We keep on mining
the next most frequent relation types recursively until we cover the required
number of the discovered edges between time series as determined by the
$\mathit{coverage}$ hyper-parameter. We remark that a pair of devices is often
connected based on multiple correlated time series, however by removing pairs
of devices at each iteration we focus on _non-redundant_ relation types by
computing the support of a relation type after filtering the device pairs
connected through higher support relation types.

Finally, we _build the network_ of connected devices based on the similarity
matrix and after filtering on the discovered frequent relations as shown in
Algorithm 2 (line 2).

####  4.5.1 Complexity

The main bottleneck is the computation of the distance function between all
$N^{2}/2$ candidates to create a symmetric similarity matrix
$\mathbf{M}\in\mathbb{R}^{N\times N}$. The proposed feature-based similarity
representation has a linear space complexity of $O(n\times k+l)$. A grid of
$n\times k$ values for the fingerprint (as determined by the hyper-parameters
$interval$ and $bins_{f}$), and $l$ bins for the histogram (as determined by
$bins_{h}$). Clearly, keeping $n\times k+l\ll|S|$ values in memory is
efficient. The algorithm for discovering frequent relation types requires
$O(N^{2})$ steps in in the worst case, however, in practice the runtime is
only a few seconds.

##  5 Pattern-based anomaly detection in a network of multivariate time series

Given the methods described in Section 3 and 4, we are now able to tackle the
main problem of detecting context-aware anomalies at the entity level. The
main Bad algorithm is shown in Algorithm 4. We characterise anomaly detection
based on the following three tasks where we compute an anomaly score for:

  1. 1.

a window $S_{i:j}$ in a _univariate time series_

  2. 2.

a window $D_{i:j}$ in a _multivariate time series_ or _entity_

  3. 3.

a window $D_{i:j}$ in an entity with regards to the _network context_.

For instance, for the last case we could discover an edge if the transmitting
network bit rate $S^{k}$ of device $D_{p}$ is highly correlated to the
receiving network bit rate $S^{l}$ of device $D_{q}$ and report an anomaly at
$D^{q}_{i:j}$ (or $D^{k}_{i:j}$) if this is not the case.

Input : A set of entities $V$ where each entity $D$ has $m$ time series (or
sensors); $D_{q}$ a target entity for anomaly detection; hyper-parameters for
each component; method to compute anomalies, i.e. Fpof or isolation forest

Result: Anomaly score for each multivariate window $D_{i:j}\in D_{q}$, or for
each univariate window $S^{k}_{i:j}$ in each sensor $S^{k}\in D^{q}$

// 1\. Discover relations between devices

1 $\mathcal{S}\leftarrow\\{S\ |\ \exists\ D\in V:S\in D\\}$

2 $E(G)\leftarrow\textsc{BuildNetwork}(\mathcal{S},\Omega_{network})$

3 $D^{q}_{\mathit{ext}}=D_{q}\cup\\{S^{j}\ |\ \exists\ S^{i}\in D_{q},\
S^{j}\in D_{k}:(D_{q},D_{k})\in E(G)\wedge S^{i}\sim S^{j}\\}$

// 2\. Create a pattern-based embedding

4 for  _$\text{sensor}\ S^{k}\in D^{q}_{\mathit{ext}}$_ do

$\mathcal{P}^{k}\leftarrow\textsc{Mip}(S,\Omega_{discrete},\Omega_{mining})$ ;

// 2.1 Mine interesting patterns

$E^{k}\leftarrow\textsc{CreateEmbedding}(\mathcal{P},S^{k})$ ;

// 2.2 Create embedding per sensor

5

// 2.3 Compute anomalies

6 if  _$\mathbf{not}\ \mathit{use\\_fpof}$_ then

7 $\mathbf{E}^{1,\ldots,m}\leftarrow\emptyset$

8 for  _$\text{window}\ (i,j)\in D^{q}_{\mathit{ext}}$_ do

9 $\mathbf{E}^{1,\ldots,m}_{i:j}=\mathit{concatenate}(\\{E^{k}_{i:j}\ |\
S^{k}\in D^{q}_{\mathit{ext}}\\})$

10 $\mathbf{return}\ \textsc{IsolationForest}(\mathbf{E})$

11else

12 $A\leftarrow\text{init matrix of anomaly scores for each sensor and
window}$

13 foreach  _$\text{sensor}\ S^{k}\in D^{q}_{\mathit{ext}}$_ do

14 $A[k]\leftarrow\textsc{Fpof}(\mathcal{P}^{k},E^{k})$

15 $\mathbf{return}\ A$

16

Algorithm 4 Bad($V$, $D_{q}$, $\Omega_{network}$, $\Omega_{discrete}$,
$\Omega_{mining}$, $\mathit{use\\_fpof}$) Efficient method to discover
anomalies in a multivariate network

###  5.1 Context-aware entity anomaly detection

First, we identify connected devices using feature-based similarity as
explained in Section 4. After we have identified connected devices we classify
anomalies in a _context-aware_ manner.

For each device $D$, we create a new collection $D_{\mathit{ext}}$ where we
add all time series of the device and similar time series of connected
devices, that is

$\displaystyle D_{\mathit{ext}}=D\cup\\{S^{j}\ |\ \exists\ S^{i}\in D,\
S^{j}\in D_{k}:(D,D_{k})\in E(G)\wedge S^{i}\sim S^{j}\\}.$

Next, we mine patterns and create an embedding for each time series and
concatenate the embedding vectors. For instance, if $S^{i}$ and $S^{j}$ are
similar (e.g., the sending and receiving bit rates are similar) we expect that
frequent patterns of $\mathcal{P}^{i}$ and $\mathcal{P}^{j}$ co-occur causing
a high anomaly score if this is not the case.

###  5.2 Anomaly detection in multivariate time series

For entity-level anomaly detection we discover patterns in each time series
_independently_ as time series are of different types and an entity
$D_{\mathit{ext}}$ (or $D$) is in essence a _heterogeneous_ collection.

####  5.2.1 Isolation forest

Given an entity $D_{\mathit{ext}}=\\{S^{1},\ldots,S^{m}\\}$ we discover
patterns in each time series separately and concatenate the discovered
embedding vectors in each time series $S^{k}$ for each window $(i,j)$ using

$E^{1,\ldots,m}_{i:j}=\mathit{concatenate}(\\{E^{k}_{i:j}\ |\ S^{k}\in D\\}).$

Next, we use an isolation forest thereby predicting a window $D_{i:j}$ as
anomalous if one or more sensors behave unexpectedly, i.e., if the joint
occurrence of patterns is unusual. The advantage of an isolation forest is
that dependencies between patterns are captured, i.e., if
$X_{1}\in\mathcal{P}^{1}$ and $X_{2}\in\mathcal{P}^{2}$ always occur together,
then an anomaly would be flagged if only $X_{1}$ occurs, which is especially
important when combining patterns in multivariate time series.

####  5.2.2 Frequent pattern outlier factor

Alternatively, we compute the anomaly score _independently_ for each time
series in $D_{\mathit{ext}}$. Ignoring interactions between patterns in
different time series can result in a lower accuracy. However, the anomaly
score can be inspected on the level of the individual time series (or
connected time series). For instance, if at window $D_{i:j}$ only time series
$S^{k},S^{n}$ and $S^{l}$ have a high anomaly score while others are normal
this _localisation_ information is useful to the end-user [30].

##  6 Experiments

In this section, we begin by evaluating our method for identifying inter-
device similarities, before comparing our anomaly detection method to state-
of-the-art algorithms on several benchmark datasets. Our implementation is
open-source and publicly
available555https://bitbucket.org/len_feremans/pbad_network/. We also provided
an online tutorial notebook in Python with many examples
666https://lfereman.github.io/.

In Table 1 we show more details on the broadcasting video network dataset.
This network consists of 315 entities and $8\,805$ time series. Each sensor
was collected over 98 days where the mean value is stored every 5 minutes
resulting in $28\,224$ values, i.e. telemetry and physical sensor data
association with a single device such as the CPU usage, CPU temperature,
memory consumption, disk usage and input and output network bitrate. Each of
the time series belong to one of 5 types of entities and there are in total 49
distinct types of time series. For instance, there are 134 entities of type A
where each entity of type A has at most 11 sensors of a certain type.

For pre-processing we cap outlier values based on the 99% quantile and apply
min-max normalisation. Many time series consist of straight lines and are less
of interest. Hence, we detect straight lines by fitting a linear function
using least squares to the time series and remove time series if the Pearson
correlation coefficient is higher 0.98.

###  6.1 Telecommunications network dataset

###  6.2 Parameter selection for the discrete representation

The accuracy of any pattern-based anomaly detection technique is highly
dependent on reasonable values for parameters determining the discrete
representation. In this section we discuss how to select parameters.

Entity type | # entities | # time series
---|---|---
A | 134 | 11
C | 60 | 20
S | 73 | 12
T | 27 | 2
D | 21 | 4
Table 1: The details of the data extracted from a broadcasting video network.
The network consists of $315$ entities of $5$ types and $8\,805$ time series,
or sensors, of $49$ types. Each time series was collected over 100 days.

####  6.2.1 Sliding window parameters

We suggest to always use an increment of 1. An increment of 1 minimising the
risk that patterns do no match a window because of arbitrary boundaries
induced by semi-overlapping sliding windows. The interval $l$ depends on the
domain but the order of magnitude of the interval, i.e., days, hours or
seconds depends on the characteristics of the dataset such as the _sampling
rate_ and the total length of the time series. In the telecom datasets we have
$12$ values logged every hour, hence an interval in hours is reasonable, i.e.,
2, 8, 24 or 48 hours. We remark for setting the interval, a measure of the
granularity in time, there is a trade-off, i.e., if the interval is selected
too course-grained anomaly detection is high, but meaningless since anomaly
scores cover a large period. Vice versa, if it is too fine-grained we have the
risk of too few data points and low accuracy.

####  6.2.2 Sax parameters

For Sax we have two parameters that control how many raw values are averaged
($\mathit{paa\\_win}$) before binning, and how many bins or the alphabet size
($\mathit{no\\_bins}$) are used. The length of the discrete sequence is given
by

$|\hat{S}_{i:j}|=\mathit{no\\_symbols}=|S_{i:j}|/\mathit{paa\\_win}=l/\mathit{paa\\_win}.$

We suggest to vary $\mathit{no\\_bins}$ and $\mathit{no\\_symbols}$ between 5
and 20 and propose defaults values of $\mathit{no\\_bins}=5$ and
$\mathit{no\\_symbols}=10$. The rationale behind this is that long sequential
patterns, i.e. consisting of more than 20 symbols are typically infrequent.
Additionally, we want a good coverage of patterns, i.e., a frequent pattern of
length 5 will cover $50\%$ of the symbols in a discrete sequence of length 10.
If labels are available, we suggest to use grid search. If labels are not
available, we suggest to plot both the time series, time series after
discretisation and the pattern-based embedding and use our human visual
expertise to select reasonable settings. We have a trade-off, i.e., if
$\mathit{no\\_bins}$ is high and $\mathit{paa\\_win}$ is small the discrete
representation will have high fidelity to the original time series, but few
frequent patterns, i.e., most subsequences will occur infrequently, especially
if there is more noise. Vice versa, if too few bins are selected and a high
value for $\mathit{paa\\_win}$ we risk that the discrete sequences are short
and no longer resembling the original time series.

###  6.3 Efficiency and accuracy of BuildNetwork

####  6.3.1 Parameter Selection

Since there is no ground truth we select parameters based on expert knowledge
and visual confirmation. For fingerprinting, we set the interval to 24 hours
and use 5 bins for discretisation. For the histogram creation, we use 100
equal-length bins. We set the threshold for fingerprinting using the 99%
quantile and for histogram distances using the 95% quantile. The coverage
parameter for filtering on frequent relations is set to 95%.

####  6.3.2 Results for frequent relation types

We find $9\,458$ raw pairs of similar time series using these parameter
settings. By filtering on intra-device similarities, edges based on known
groups of devices, and frequent relation types, only $351$ edges remain. In
Table 2 we show the top-6 of the total $28$ frequent relation types which is
much smaller than the $\binom{49}{2}=1176$ theoretically possible
combinations. These frequent relations are a good starting point to
_comprehend_ this complex network of heterogeneous devices and sensors.

Type |  # edges
---|---
A/System load |  91
A/CPU load |  52
A/IP bit rate |  36
A/TS rate |  29
A/User load |  24
A/CPU System - A/CPU system time load |  10
Table 2: Most frequent relation types between devices

####  6.3.3 Runtime performance

Computing the similarity of all $78\cdot 10^{6}$ possible pairs _takes just
$8.5$_ minutes or $59.0$ ms per instance to search in $8\,708$ time series
after excluding straight lines.

####  6.3.4 Comparing the accuracy of similarity measures

We compare the proposed fingerprint-histogram distances against Paa [19] and
the Pearson correlation coefficient after sampling each time series to 2000
values. To evaluate performance we compute precision, recall and F1. Since we
do not have labels of similar time series, we pick 50 time series, search for
the top-5 most similar time series for each of them, and manually assign
labels based on domain knowledge. The results are shown in Table 3. We find
that fingerprinting combined with histogram distances has a _higher precision_
and _F1_ but _lower recall_ compared to Paa and that the performance of
Pearson correlation is relatively low. However, we argue that for context-
aware anomaly detection high precision is of more importance than high recall.

Method | Precision | Recall | F1
---|---|---|---
Fingerprinting and histogram | $\mathbf{0.971}$ | $0.779$ | $\mathbf{0.864}$
Piecewise aggregate approximation | $0.706$ | $\mathbf{0.895}$ | $0.789$
Pearson correlation coefficient | $0.462$ | $0.639$ | $0.536$
Table 3: Accuracy of time series similarity identification in the broadcasting
video network dataset

####  6.3.5 Qualitative performance

In Figure 4 we show the resulting network visualised using the NetworkX
package. We find 4 large dense clusters of mostly “type A” devices which are a
type of card located inside the same “type C”, so, unsurprisingly, they have
multiple similar sensors and a connection with a few “type C” elements. We
also discover interesting smaller clusters, e.g., “type D”, and “type C”
devices are often connected using bit rate. Furthermore, we observe that most
discovered relations are between devices where the identifier is near, e.g.,
device 1268 is connected with 1269. This validates that the proposed algorithm
for identifying edges using only time series data matches with the prior
knowledge that devices with a near human-assigned identifier are often
related.

Figure 4: Discovered network of broadcasting video devices. The colour of each
node represents a known cluster and the shape the type of device

###  6.4 Efficiency and accuracy of pattern-based anomaly detection on
benchmark datasets

Since we have no labels for the network dataset, we evaluate our anomaly
detection method on univariate and multivariate benchmark datasets.

####  6.4.1 Baselines

We compare Bad with two state-of-the-art pattern-based anomaly detection
methods: _Frequent pattern based outlier factor_ (Fpof) computes an outlier
factor based on the number of frequent closed itemsets [14], and _Pattern-
based anomaly detection_ (Pbad) uses an isolation forest trained on the
embedding of raw sensor values, closed frequent itemsets and sequential
patterns [8]. Additionally, we compare with _Isolation forest_ that uses a
forest trained on the original sensor values in each window [25] and Bad-Fpof
which uses the frequent pattern based outlier factor instead of an isolation
forest to facilitate interpretability.

####  6.4.2 Evaluation protocol

Given a time series $S$ with a number of labelled timestamps, we divide the
time series into fixed-sized sliding windows and compute an anomaly score for
each window. Similar to the evaluation protocol of OmnyAnomaly [30], we adopt
the _point-adjust_ approach where we consider an anomaly detected if it occurs
within the boundary of a sliding window. Each method requires the setting of a
threshold and we assume an Oracle that reports the best possible F1 value by
varying the threshold. Since the isolation forest is not deterministic we run
each method 5 times using the same parameters and report the mean values.

Dataset | IsoForest | Fpof | Pbad | Bad | Bad-Fpof
---|---|---|---|---|---
Temp | 0.927 | 0.875 | 0.776 | 0.948 | 0.881
Taxi | 0.782 | 0.718 | 0.807 | 0.851 | 0.630
Latency | 0.827 | 0.653 | 0.864 | 0.901 | 0.845
Water 1 | 0.686 | 0.333 | 0.755 | 0.787 | 0.549
Water 2 | 0.528 | 0.450 | 0.596 | 0.539 | 0.452
Water 3 | 0.791 | 0.321 | 0.615 | 0.567 | 0.447
Water 4 | 0.561 | 0.661 | 0.760 | 0.666 | 0.663
_Avg. rank univariate_ | 3 | 4.4 | 2.1 | 1.6 | 3.8
SMD | 0.949 | n/a | n/a | 0.966 | 0.957
SMAP | 0.708 | n/a | n/a | 0.775 | 0.733
_Avg. rank multivariate_ | 3 | n/a | n/a | 1 | 2
Table 4: F1 score of various methods on $7$ univariate datasets and $2$
multivariate datasets consisting of $28$ and $54$ entities respectively, where
we report the mean result

####  6.4.3 Parameter selection

Each method has the same preprocessing steps which include setting the window
size $l$ and increment $t$ to create sliding windows and parameters for
transforming to a symbolic representation such as the _paa_window_ for
reducing the dimensionality and the number of equal-length _bins_ which we
optimise using _grid search_. For the IsolationForest, we use 500 _trees_. For
Pbad we set $\mathit{min\\_len}$ to 2 and the _minimal support_ to $0.01$. For
Fpof, we set the _minimal support_ to $0.01$. For Bad, we set the number of
patterns $k$ to $10\,000$, $\mathit{min\\_len}$ to 2 and the _relative
duration_ to $1.2$.

####  6.4.4 Datasets

For the univariate test case, we use 7 datasets. Temp, Latency and Taxi are
smaller datasets from the Numenta Anomaly Benchmark [23]. Water 1 to 4 is a
proprietary dataset that contains the average consumption of water in 4
branches of a retail company where the goal is to detect abnormal water
consumption [8]. The water consumption is logged every 5 minutes for two years
resulting in $170\,000$ values. For the multivariate test case, we use the
Server Machine Dataset (SMD) which consists of 28 entities, with 38 sensors
each, containing $56\,960$ values [30]. The semantics are similar to our
broadcasting video network with data such as CPU load, RAM usage, etc. Soil
Moisture Active Passive (SMAP) is a dataset from NASA and consists of 54
entities each monitored by 25 metrics [16].

####  6.4.5 Compare the accuracy with baselines

Table 4 shows the F1 obtained by each method on each dataset. Results in bold
highlight the best performing method on each dataset. We find that Bad has the
lowest average rank on univariate datasets thus outperforming existing state-
of-the-art methods such as IsoForest, Pbad and Fpof. We find that both Bad and
Bad-Fpof perform better than IsoForest on the multivariate datasets. We report
no results for Fpof and Pbad on the multivariate datasets as the run times are
prohibitively long.

We compare with the reported results of recent neural-network based methods,
such as OmniAnomaly and Usad [2]. On the SMD dataset, we find the average F1
value of Bad is $0.966$ $(\pm 0.002)$ which is higher than OmniAnomaly (0.962)
and Usad (0.946). However, the F1 value for SMAP is lower, i.e. $0.775$ for
Bad versus $0.853$ for OmnyAnomaly and $0.863$ for Usad. The relatively worse
performance on SMAP is possibly explained by the special nature of this
dataset, where time series often consist of only peaks alternating between two
values. We remark that all existing methods ignore context when considering
anomaly detection in a network of devices.

####  6.4.6 Runtime performance

Concerning runtime performance we find that Pbad is considerably slower and
requires more than 1.5 hours to complete on each Water dataset. The minimum
support threshold results in an unpredictable and larger set of patterns. In
Bad using default parameters results in fewer patterns (800 instead of 5000
for the Water dataset) and both pattern discovery and the creation of the
embedding are more efficient resulting in a total runtime of less than a few
minutes. On SMD the total time for pre-processing, learning the model and
making predictions was just 37 _minutes_ (on a laptop). This is in stark
contrast to OmniAnomaly, that required more than two hours to train a model
for a _single_ device on the SMD datasets. We conclude that Bad and IsoForest
scale to long time series and collections thereof.

####  6.4.7 Qualitative evaluation

For decision support, we use Bad-Fpof to allow _localisation_ by computing the
Fpof-based anomaly score based on patterns discovered in each time series. In
Figure 5 we visualise the pattern-based embedding of the 7th sensor of the
first device in the SMD dataset where the Fpof score is based on only 6
patterns. The known anomaly is located after the first series of
$X_{1}=[2,2,2,\ldots]$ occurrences and before occurrences of
$X_{6}=[3,3,3,\ldots]$, which is also when few patterns occur causing a high
anomaly score.

Figure 5: Example of a discrete representation of metric-7 from the SMD
dataset where we show the pattern-based embedding. The pattern mining
algorithm discovers only 6 patterns. We show the occurrences of each pattern
and remark that in this case a high anomaly score (red) is correlated with an
absence of pattern occurrences.

In Figure 6 we show 4 time series from the first device in the SMD dataset.
For SMD we have labelled anomalies specific to a subset of time series
highlighted in red. Below each time series we show the discrete representation
computed using an interval of 6 hours, 10 symbols and 10 equal-length bins and
the Fpof score in orange. We remark that there is high correspondence between
known and predicted anomalies.

Figure 6: Example of localised anomaly scores of 4 time series from a single
device in the SMD dataset using Bad-Fpof . We show the time series before and
after discretisation, the known sensor-specific anomalies (in red), and the
predicted anomalies (in orange)

##  7 Related work

We have shown that the proposed similarity measure that considers the trend,
peaks and the distribution of values results in a simple efficient method that
produces fewer false positive than Paa. For comparing time series more
advanced techniques have been proposed, however, few scale to large time
series databases. Finally, the approach for filtering spurious relations based
on the type of sensor is of specific interest to intra-device similarities and
has not been considered before.

For pattern-based anomaly detection we have shown that Bad is more accurate
and efficient compared with state-of-the-art methods, such as Fpof and Pbad.
Related work, such as MivPod [15] was shown to be less accurate than Pbad.

In contrast to prior work in sequential pattern mining, we are interested in
patterns that are both frequent, cohesive and compress the data series, while
previous work considers only one or two aspects, i.e. [31, 22] and [29]
discovers a set of patterns the compresses the dataset best, while [7]
discover frequent and cohesive patterns for data series.

Recently, authors proposed to make complex neural network-based decisions
explainable, thereby providing methods to highlight which segments and sensors
of a device contribute to a predicted anomaly [2]. In contrast, we propose an
_intrinsically_ interpretable white-box model that consists of a reduced set
of interpretable features (or patterns). This enables oversight from human
experts that can inspect patterns to trust decisions both case-by-case and
globally.

##  8 Conclusions

In this paper, we present a method for detecting anomalies in networks of
connected devices. To the best of our knowledge, no existing method is capable
of dealing with this problem. In the first step, we propose a novel algorithm
for detecting similar time series (of potentially different types of sensors)
in different devices. In this way, we identify how devices interact with each
other, as this information is often dynamic or unavailable. Once we have built
a network of devices, we discover frequent cohesive patterns in the time
series to describe normal device behaviour. To improve comprehension, we
select the best patterns using the minimum description length principle. We
then transform the time series into a pattern-based embedding, in which we
search for anomalies. Intuitively, we see an anomaly as a time period in which
frequent patterns are absent. By taking interactions between time series on
the same device and various connected devices into account, we can capture
contextual anomalies that would otherwise potentially be missed. In addition
to being able to handle time series from a network of connected devices, our
experiments show that our anomaly detection method also outperforms existing
methods on both univariate and multivariate benchmark datasets.

## Acknowledgements

The authors would like to thank VLAIO CONSCIOUS (Anomaly Detection for Complex
Industrial Assets) project for funding this research.

## References

  * [1] Komate Amphawan and Philippe Lenca.  Mining top-k frequent-regular closed patterns.  Expert Systems with Applications, 42(21):7882–7894, 2015.
  * [2] Julien Audibert, Pietro Michiardi, Frédéric Guyard, Sébastien Marti, and Maria A Zuluaga.  Usad: unsupervised anomaly detection on multivariate time series.  In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 3395–3404, 2020.
  * [3] Kasun Bandara, Christoph Bergmeir, and Slawek Smyl.  Forecasting across time series databases using recurrent neural networks on groups of similar series: A clustering approach.  Expert systems with applications, 140:112896, 2020.
  * [4] Varun Chandola, Arindam Banerjee, and Vipin Kumar.  Anomaly detection: A survey.  ACM computing surveys, 41(3):1–58, 2009.
  * [5] Rongjun Chen, Xuanhui Yan, Shiping Wang, and Guobao Xiao.  Da-net: Dual-attention network for multivariate time series classification.  Information Sciences, 610:472–487, 2022.
  * [6] Boris Cule, Len Feremans, and Bart Goethals.  Efficiently mining cohesion-based patterns and rules in event sequences.  Data Mining and Knowledge Discovery, 33(4):1125–1182, 2019.
  * [7] Len Feremans, Boris Cule, and Bart Goethals.  Petsc: pattern-based embedding for time series classification.  Data Mining and Knowledge Discovery, pages 1–47, 2022.
  * [8] Len Feremans, Vincent Vercruyssen, Boris Cule, Wannes Meert, and Bart Goethals.  Pattern-based anomaly detection in mixed-type time series.  In Joint European conference on machine learning and knowledge discovery in databases, pages 240–256. Springer, 2019.
  * [9] Philippe Fournier-Viger, Antonio Gomariz, Ted Gueniche, Espérance Mwamikazi, and Rincy Thomas.  Tks: efficient mining of top-k sequential patterns.  In International Conference on Advanced Data Mining and Applications, pages 109–120. Springer, 2013.
  * [10] Philippe Fournier-Viger, Jerry Chun-Wei Lin, Antonio Gomariz, Ted Gueniche, Azadeh Soltani, Zhihong Deng, and Hoang Thanh Lam.  The spmf open-source data mining library version 2.  In Joint European conference on machine learning and knowledge discovery in databases, pages 36–40. Springer, 2016.
  * [11] Johannes Fürnkranz.  Separate-and-conquer rule learning.  Artificial Intelligence Review, 13(1):3–54, 1999.
  * [12] Peter D Grünwald.  The minimum description length principle.  MIT press, 2007.
  * [13] Jiawei Han, Jian Pei, Behzad Mortazavi-Asl, Helen Pinto, Qiming Chen, Umeshwar Dayal, and MC Hsu.  Prefixspan: Mining sequential patterns efficiently by prefix-projected pattern growth.  In proceedings of the 17th international conference on data engineering, pages 215–224, 2001.
  * [14] Zengyou He, Xiaofei Xu, Zhexue Joshua Huang, and Shengchun Deng.  Fp-outlier: Frequent pattern based outlier detection.  Computer Science and Information Systems, 2(1):103–118, 2005.
  * [15] C Sweetlin Hemalatha, Vijay Vaidehi, and R Lakshmi.  Minimal infrequent pattern based approach for mining outliers in data streams.  Expert Systems with Applications, 42(4):1998–2012, 2015.
  * [16] Kyle Hundman, Valentino Constantinou, Christopher Laporte, Ian Colwell, and Tom Soderstrom.  Detecting spacecraft anomalies using lstms and nonparametric dynamic thresholding.  In Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining, pages 387–395, 2018.
  * [17] Rob J Hyndman, Earo Wang, and Nikolay Laptev.  Large-scale unusual time series detection.  In 2015 IEEE international conference on data mining workshop (ICDMW), pages 1616–1619. IEEE, 2015.
  * [18] Herve Jegou, Matthijs Douze, and Cordelia Schmid.  Product quantization for nearest neighbor search.  IEEE transactions on pattern analysis and machine intelligence, 33(1):117–128, 2010.
  * [19] Eamonn Keogh, Kaushik Chakrabarti, Michael Pazzani, and Sharad Mehrotra.  Dimensionality reduction for fast similarity search in large time series databases.  Knowledge and information Systems, 3(3):263–286, 2001.
  * [20] Tung Kieu, Bay Vo, Tuong Le, Zhi-Hong Deng, and Bac Le.  Mining top-k co-occurrence items with sequential pattern.  Expert Systems with Applications, 85:123–133, 2017.
  * [21] Justin B Kinney and Gurinder S Atwal.  Equitability, mutual information, and the maximal information coefficient.  Proceedings of the National Academy of Sciences, 111(9):3354–3359, 2014.
  * [22] Hoang Thanh Lam, Fabian Mörchen, Dmitriy Fradkin, and Toon Calders.  Mining compressing sequential patterns.  Statistical Analysis and Data Mining: The ASA Data Science Journal, 7(1):34–52, 2014.
  * [23] Alexander Lavin and Subutai Ahmad.  Evaluating real-time anomaly detection algorithms–the numenta anomaly benchmark.  In Machine Learning and Applications (ICMLA), 2015 IEEE 14th International Conference on, pages 38–44. IEEE, 2015.
  * [24] Jessica Lin, Eamonn Keogh, Stefano Lonardi, and Bill Chiu.  A symbolic representation of time series, with implications for streaming algorithms.  In Proceedings of the 8th ACM SIGMOD workshop on Research issues in data mining and knowledge discovery, pages 2–11. ACM, 2003.
  * [25] Fei Tony Liu, Kai Ming Ting, and Zhi-Hua Zhou.  Isolation forest.  In 2008 Eighth IEEE International Conference on Data Mining, pages 413–422. IEEE, 2008.
  * [26] Florent Masseglia, Pascal Poncelet, and Maguelonne Teisseire.  Efficient mining of sequential patterns with time constraints: Reducing the combinations.  Expert Systems with Applications, 36(2):2677–2690, 2009.
  * [27] Atif Raza and Stefan Kramer.  Accelerating pattern-based time series classification: a linear time and space string mining approach.  Knowledge and Information Systems, 62(3):1113–1141, 2020.
  * [28] David Savage, Xiuzhen Zhang, Xinghuo Yu, Pauline Chou, and Qingmai Wang.  Anomaly detection in online social networks.  Social networks, 39:62–70, 2014.
  * [29] Mohammad Shokoohi-Yekta, Yanping Chen, Bilson Campana, Bing Hu, Jesin Zakaria, and Eamonn Keogh.  Discovery of meaningful rules in time series.  In Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining, pages 1085–1094, 2015.
  * [30] Ya Su, Youjian Zhao, Chenhao Niu, Rong Liu, Wei Sun, and Dan Pei.  Robust anomaly detection for multivariate time series through stochastic recurrent neural network.  In Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, pages 2828–2837, 2019.
  * [31] Jilles Vreeken, Matthijs Van Leeuwen, and Arno Siebes.  Krimp: mining itemsets that compress.  Data Mining and Knowledge Discovery, 23(1):169–214, 2011.
  * [32] Jilles Vreeken and Kenji Yamanishi.  Modern mdl meets data mining insights, theory, and practice.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 3229–3230, 2019.
  * [33] Miao Xie, Song Han, Biming Tian, and Sazia Parvin.  Anomaly detection in wireless sensor networks: A survey.  Journal of Network and computer Applications, 34(4):1302–1325, 2011\.
  * [34] Mohammed J Zaki and Wagner Meira.  Data mining and analysis: fundamental concepts and algorithms.  Cambridge University Press, 2014.