# High-performance Time Series Anomaly Discovery
on Graphics Processor

Mikhail L. Zymbler111South Ural State University, Chelyabinsk, Russian
Federation Yana A. Kraevasusu

###### Abstract

Currently, discovering subsequence anomalies in time series remains one of the
most topical research problems. A subsequence anomaly refers to successive
points in time that are collectively abnormal, although each point is not
necessarily an outlier. Among a large number of approaches to discovering
subsequence anomalies, the discord concept is considered one of the best. A
time series discord is intuitively defined as a subsequence of a given length
that is maximally far away from its non-overlapping nearest neighbor. Recently
introduced the MERLIN algorithm discovers time series discords of every
possible length in a specified range, thereby eliminating the need to set even
that sole parameter to discover discords in a time series. However, MERLIN is
serial and its parallelization could increase the performance of discords
discovery. In this article, we introduce a novel parallelization scheme for
GPUs, called PALMAD, Parallel Arbitrary Length MERLIN-based Anomaly Discovery.
As opposed to its serial predecessor, PALMAD employs recurrent formulas we
have derived to avoid redundant calculations, and advanced data structures for
the efficient implementation of parallel processing. Experimental evaluation
over real-world and synthetic time series shows that our algorithm outperforms
parallel analogs. We also apply PALMAD to discover anomalies in a real-world
time series employing our proposed discord heatmap technique to illustrate the
results.

###### keywords:

time series, anomaly detection, discord, MERLIN, DRAG, parallel algorithm,
GPU, CUDA

## Introduction

Over past decades, time series data are ubiquitous in diverse spheres of a
human’s activity: industry, healthcare, science, social, and so on. Currently,
discovering anomalies (or outliers as a synonym) in time series remains one of
the most topical research problems. In a time series, point and subsequence
anomalies can serve as the aims to be detected [4]. The former defines a datum
that deviates in a specific time instant when compared either to the other
values in the time series or to its neighboring points. The latter refers to
successive points in time whose collective behavior is unusual, although each
observation individually is not necessarily a point anomaly. Subsequence
anomaly detection is more challenging due to the need to take into account the
subsequence length among other aspects [4].

Among a wide spectrum of analytical and neural network-based approaches to
time series subsequence anomaly detection [4, 10], the discord concept [31] is
considered one of the best [7, 8]. A time series discord is intuitively
defined as a subsequence that is maximally far away from its non-overlapping
nearest neighbor. Discords look attractive for an end-user since they require
the only parameter to be specified, the subsequence length. However,
application of discords is reduced by sensitivity to this single user choice.
A straightforward solution of this problem, namely discovering discords of all
the possible lengths and then selecting the best discords with respect to some
measure looks computationally prohibitive.

Nevertheless, recently introduced the MERLIN algorithm [36] can efficiently
and exactly discover discords of every possible length in a specified range,
being ahead of competitors in terms of accuracy and performance. Thus, MERLIN
allows an end-user for removing the need to set even that above-mentioned sole
parameter to discover discords in a time series. MERLIN employs repeated calls
of the DRAG algorithm [51] that discovers discords of a given length with a
distance of at least $r$ to their nearest neighbors, and adaptive selection of
the parameter $r$. However, multiple calls of the above sub-algorithm result
in calculations that are partially repeated and being executed at once would
increase the performance of MERLIN. Furthermore, the MERLIN algorithm is
serial and looks attractive for parallelization to increase the performance of
discords discovery.

In this study, we address the problem of parallelization the MERLIN algorithm
for discovery of arbitrary length discords on GPU continuing our research on
accelerating various time series mining tasks with parallel architectures and
in-database time series analysis [29, 58, 62, 60, 65, 63, 61, 59]. The
article’s contribution can be summarized as follows. We perform thorough
review of works related to the discord-based approaches to discovering time
series anomalies and their parallelization for diverse hardware platforms.
Next, based on MERLIN [36], we introduce the parallel algorithm PALMAD
(Parallel Arbitrary Length MERLIN-based Anomaly Discovery) for discovery of
arbitrary length discords on a graphics processor. PALMAD employs recurrent
formulas we have derived to avoid redundant calculations, and advanced data
structures for the efficient implementation of parallel processing. Further,
in the extensive experimental evaluation over real-world and synthetic time
series, we show that our algorithm outperforms existent analogs. Finally, we
apply PALMAD to anomaly discovery in a real-world time series employing our
proposed discord heatmap technique to illustrate the results.

The remainder of the article is organized as follows. In Section 1, we discuss
related works. Section 2 contains notation and formal definitions along with a
short description of the original serial algorithm. Section 3 introduces the
proposed parallel algorithm to discover time series discords on GPU. In
Section 4, we give the results and discussion of the experimental evaluation
of our algorithm. Section 5 describes a case study on discord discovery in a
real-world time series. Finally, in Conclusions, we summarize the results
obtained and suggest directions for further research.

##  1 Related Work

The time series discord concept was introduced by Keogh et al. [31], and
currently, is considered as one of the best analytical approaches to
discovering anomalies in time series [7, 8]. Time series discord is
intuitively defined as the subsequence of a time series that is the most
distant to its non-overlapping nearest neighbor. Discords look attractive for
an end-user since they require the only parameter to be specified, the
subsequence length. Below, we consider discord-based research on discovering
anomalies including studies addressed the parallelization of the approaches
above.

In [31], Keogh et al. proposed the HOTSAX (Heuristically Ordered Time series
using Symbolic Aggregate ApproXimation) algorithm for discords discovery in a
time series that can be entirely placed in RAM. HOTSAX employs time series
encoding through the SAX technique [32] and the Euclidean distance. HOTSAX
iterates through all the pairs of subsequences calculating the distance
between them, and finds the maximum among the distances to the nearest
neighbor. The algorithm employs a prefix trie [15] to index the subsequences.
When iterating, unpromising subsequences are discarded without calculating
distances. The subsequence with a neighbor closer than the best-so-far maximum
of distances to all the nearest neighbors is unpromising. HOTSAX exploits a
certain heuristic that allows to discarding more unpromising candidates.
Improvements of HOTSAX include iSAX [41] and HOT-iSAX [6] (indexable SAX), WAT
[5, 16] (application of the Haar wavelets instead of SAX and augmented trie),
HashDD [45] (employing a hash table instead of the prefix trie), HDD-MBR [9]
(application of R-trees), BitClusterDiscord [30] (employing clustering of the
bit representation of subsequences), and HST (HOTSAX Time) [3] (reduction the
size of the discord search space through the warm-up process and the
similarity between subsequences close in time).

In [40], Senin et al. proposed the HOTSAX-based RRA (Rare Rule Anomaly)
algorithm to discover variable-length discords. RRA deals with a time series
discretized with SAX applying grammar-induction procedures. Since symbols
infrequently used in grammar rules are non-repetitive and thus potentially
unusual, the discords correspond to infrequent grammar rules that vary in
length. Taking into account that the lengths of the subsequences vary, the
distance between them is calculated by shrinking the longest subsequence with
the Piecewise Aggregate Approximation (PAA) [23] to obtain subsequences of the
same length. Although the RRA algorithm is a step forward from HOTSAX to
parameter-free discord discovering, as its predecessor, it is limited to RAM-
stored time series.

In [51], Yankov, Keogh et al. presented the DRAG (Discord Range Aware
Gathering) algorithm for discovering discords in a time series stored on a
disk rather than in RAM. DRAG introduces the range discord concept, where such
a discord has a distance of at least $r$ to its non-overlapping nearest
neighbor, and $r$ is a user-defined threshold. The DRAG algorithm performs in
two phases, namely the candidate selection (collecting potential range
discords) and discord refinement (discarding false positives), with each phase
requiring one linear scan through the time series on the disk. In [43], Son
slightly improved DRAG through the employing a hash bucket data structure to
speed up the candidate selection phase. The authors of DRAG proposed the
following procedure to choose the parameter $r$. Through the uniform sampling,
one can obtain a maximum length fragment of the original time series that fits
in RAM. Next, the HOTSAX algorithm discovers a discord in the above-obtained
fragment. Finally, the $r$ threshold is assumed to be equal to the distance of
the discord found to its nearest neighbor.

However, in the DRAG algorithm, the above-described heuristic does not define
a formal way to choose the parameter $r$ to guarantee the efficiency of
discords discovery [34]. Ideally, $r$ should be set in such a way that it is a
little less than the distance between the discord eventually found and its
nearest neighbor [36]. Then the time and space complexity of DRAG is $O(mn)$,
where $n$ is the length of the time series, and $m$ is the discord length. If
the value of $r$ is set significantly less than the above-mentioned distance,
then the algorithm will find the discord, but the time and spatial complexity
will be higher, $O(n^{2})$. Finally, if $r$ is greater than the above
distance, then no discords will be found. In addition, DRAG (like HOTSAX as
its predecessor) is not able to discover _all_ the discords in the following
meaning: the algorithm finds discords of a single specified length, but not
discords of every possible length in a specified range. In the latter case, a
brute-force approach involving a cyclic runs of the DRAG algorithm for a
specified range of the discord length does not work, since at each iteration
of such a loop, we should choose the parameter $r$ from scratch.

Recently introduced by Keogh et al. the MERLIN algorithm [36] overcomes the
above-described limitations of DRAG. MERLIN calls DRAG repeatedly and
adaptively selects the parameter $r$. In the experiments, MERLIN efficiently
and exactly discovers discords of every possible length, being ahead of
competitors both in accuracy and performance [36]. The authors also
empirically showed that MERLIN is able to discover point, contextual, and
collective anomalies according to the taxonomy in [7]. Moreover, the authors
mentioned that despite the recent explosion of the deep learning anomaly
detection methods (e.g., [2, 20, 13, 35]), it is not obvious that they
outperform discord-based approaches, since the former, by its nature, require
many critical parameters to be set whereas the latter are domain-independent
and require one intuitive parameter the need to set which can even be removed
by MERLIN. However, the MERLIN algorithm is still serial, and its
parallelization (for various hardware platforms) could increase the
performance of discords discovery. In addition, let us mention that, in
MERLIN, multiple calls of DRAG result in calculations that are partially
repeated (e.g., normalization of the subsequences in a specified range of
length) and being executed at once would also increase the performance of
MERLIN.

Research addressed the problem of parallelization the discord discovery
include the following. In [65, 58], Zymbler et al. proposed a parallelization
schema for HOTSAX to discover discords with Intel many-core processors or GPUs
through the OpenMP [44] or OpenACC [39] technology, respectively. The
algorithm employs matrix data layout to organize calculations with as many
vectorizable loops as possible. Similarly to its predecessor, the algorithm
differs the following sets of subsequences: ones with the least frequent SAX
words and the rest, and for any subsequence – ones whose SAX words match the
given subsequence’s SAX word and the rest. When iterating all the subsequences
through two nested loops, the algorithm parallelizes separately and
differently for the outer and the inner loops, depending on the number of
running threads and the cardinality of the above-mentioned sets.

In [52] (an expanded version of [51]), Yankov, Keogh et al. discussed the
parallel version of DRAG based on the MapReduce paradigm [11], and the key
idea is as follows. Let the input time series be partitioned evenly across $P$
high-performance cluster nodes. Each node selects candidates in its own
partition with the same parameter $r$ resulting in the local candidate set
$\mathcal{C}_{i}$. Then the global candidate set $\mathcal{C}$ is constructed
as $\mathcal{C}=\cup^{P}_{i=1}\mathcal{C}_{i}$ and sent to each cluster node.
Next, a node refines candidates in its own partition taking the global
candidate set $\mathcal{C}$ as an input, and produces the local refined
candidate set $\tilde{\mathcal{C}}_{i}$. Finally, the global discords set
$\mathcal{D}$ is computed as
$\mathcal{D}=\cap^{P}_{i=1}\tilde{\mathcal{C}}_{i}$. In the experimental
evaluation, the authors, however, just simulated the above-mentioned scheme on
up to eight computers resulting in a close-to-linear speedup.

In [60], Zymbler et al. introduced a parallelization scheme of DRAG for a
high-performance cluster with Intel many-core processors. For the cluster
node, the authors define matrix data structures and employ thread-level
parallelism through the OpenMP technology [44] whereas communication among the
cluster nodes is implemented through MPI (Message Passing Interface) [42]. As
opposed to DRAG parallelization scheme, at each cluster node, the authors
first refine the local candidate set $\mathcal{C}_{i}$ with respect to the
same parameter $r$ resulting in the $\tilde{\mathcal{C}}_{i}$ set, and then
construct the global candidate set as
$\mathcal{C}=\cup^{P}_{i=1}\tilde{\mathcal{C}}_{i}$ relying on the fact that a
candidate is not a true discord if it was pruned by at least one cluster node
during the selection phase. In the experiments [60], the authors showed that
such a technique allows for significant reduction of the global candidate set
and increasing the overall algorithm’s performance. The algorithm also
significantly outperforms the following DRAG-based parallel discords discovery
algorithms for high-performance clusters: DDD (Distributed Discord Discovery)
[50] and PDD (Parallel Discord Discovery) [18]. The above-mentioned
competitors are far behind since the fact that they involve intensive data
exchanges across cluster nodes. However, the above-described parallelization
still cannot efficiently discover discords of every possible length in a
specified range.

In the review, we should also mention the matrix profile (hereinafter MP)
concept proposed by Keogh et al. [53]. For a given time series, MP can
informally be defined as a time series, where the $i$-th element is the
distance from the $i$-th subsequence of the original time series to its non-
overlapping nearest neighbor. MP plays the role of a building block the
solutions of various time series motif discovery related problems are based on
(semantic motifs [21], snippets [22], chains [55], etc.). According to the
above definition, top-$k$ discords can be discovered as a by-product of the MP
calculation since they are the subsequences on which the top-$k$ maximum
values in MP are achieved. However, the time complexity of MP computation is
high, namely $O(n^{2})$ (where $n$ is the time series length) [53, 56], so
straightforward employing of MP in the discords discovery results in low
performance as it has been evaluated in the following experiments. The serial
SCRIMP algorithm [53] is inferior to MERLIN [36]. Parallel MP algorithms for
graphics processor and high-performance cluster, GPU-STAMP [53] and MP-HPC
[38], respectively, are inferior to the parallel discords discovery algorithm
for high-performance cluster with Intel many-core processors proposed in [60].

In [46], Thuy et al. introduced the notion of the $K$-distance discord, namely
a subsequence with the largest sum of distances to its non-overlapping $K$
nearest neighbors. Such an approach aims at solving so-called “twin freak”
problem [48] when a discord fails to discover an anomalous (rare) subsequence
if it occurs more than once in the time series, and is a modification of the
J-distance discord [19] concept, where the distance between a subsequence and
its $k$-th non-overlapping nearest neighbor is employed. The authors also
presented the KBF_GPU (Brute-Force for K-distance discord) algorithm that
accelerates $K$-distance discord discovery on a graphics processor. KBF_GPU
iterates all the subsequences of a given time series through two nested loops
where the inner loop is parallelized and adapted to calculate sum of
distances. In the experiments, the authors, however, compare their algorithm
only with serial HOTSAX [31], and the latter, as expected, is significantly
inferior to KBF_GPU.

In [54], Zhu et al. presented a parallel algorithm to accelerate discords
discovery with GPU. The authors exploit the normalized Euclidean distance and
its efficient calculation through the Pearson correlation by the technique
proposed in [34]. To provide high performance of discord discovery, the
algorithm employs two computational patterns. The first one prescribes the
following two-step procedure. First, calculate the minimum distance between
the discord candidate subsequence and all other subsequences of the time
series that do not overlap the candidate. Then, find a candidate on which the
maximum distance among all the candidates is achieved. The second pattern
assumes an early stop of calculations in the pattern above when the distance
between the candidate and a certain subsequence is less than the best-so-far
distance. In such a case, both the candidate and the subsequence are obviously
not discords, and we do not need to calculate distances from the candidate to
other non-overlapping subsequences. In the experiments [54], the proposed
algorithm outruns SCAMP [57] that is currently the fastest parallel algorithm
for calculating the matrix profile. However, the proposed computational
patterns limit the result to a single (albeit the most important) discord of
the time series, whereas the above-described algorithms are based on the range
discord concept and able to discover top-$k$ discords, where the parameter $k$
is prespecified by an expert in the subject domain.

Concluding our overview of related work, it can be seen that, currently, the
MERLIN algorithm [36] based on the range discord concept [51] is one of the
most promising approaches to discover anomalies in time series. Moreover,
being analytical and agnostic, MERLIN is at least competitive with deep
learning methods. However, parallelization of MERLIN could increase the
performance of discord discovery. Such parallelization is a topical issue
since, to the best of our knowledge, no research has addressed the
accelerating discovery of discords of every possible length with GPU or any
other parallel hardware architecture.

##  2 Preliminaries

Prior to detailing the proposed parallel algorithm for discords discovery, in
Sections 2.1 and 2.2, we introduce basic notation and formal definitions
according to [36, 51] and give an overview of the original serial algorithms
MERLIN and DRAG our development is based on, respectively.

###  2.1 Notation and Definitions

A _time series_ is a chronologically ordered sequence of real-valued numbers:

$T=\\{t_{i}\\}_{i=1}^{n},\quad t_{i}\in\mathds{R}.$ (1)

The length of a time series, $n$, is denoted by $\lvert T\rvert$. Hereinafter,
we assume that the time series $T$ fit into the main memory.

A _subsequence_ $T_{i,\,m}$ of a time series $T$ is its subset of $m$
successive elements that starts at the $i$-th position:

$T_{i,\,m}=\\{t_{k}\\}_{k=i}^{i+m-1},\quad 1\leq i\leq n-m+1,\quad 3\leq m\ll
n.$ (2)

We denote the set of all $m$-length subsequences in $T$ by $S^{m}_{T}$. Let
$N$ denotes the number of subsequences in $S^{m}_{T}$, i.e.,
$N=|S^{m}_{T}|=n-m+1$.

A _distance function_ for any two $m$-length subsequences is a nonnegative and
symmetric function
$\operatorname{Dist}:\mathds{R}^{m}\times\mathds{R}^{m}\rightarrow\mathds{R}$.

Given a time series $T$ and its two subsequences $T_{i,\,m}$ and $T_{j,\,m}$,
we say that they are _non-self match_ to each other at distance
$\operatorname{Dist}(T_{i,\,m},\;T_{j,\,m})$ if $|i-j|\geq m$. Let us denote a
non-self match of a subsequence $C\in S^{m}_{T}$ by $M_{C}$.

Given a time series $T$, its subsequence $D\in S^{m}_{T}$ is said to be the
_discord_ if $D$ has the largest distance to its nearest non-self match.
Formally speaking, the discord $D$ meets the following:

$\forall C\in
S^{m}_{T}\quad\min\bigl{(}\operatorname{Dist}(D,\,M_{D})\bigr{)}>\min\bigl{(}\operatorname{Dist}(C,\,M_{C})\bigr{)}.$
(3)

The definition above is generalized from top-1 to top-$k$ discord as follows:
$D\in S^{m}_{T}$ is said to be the $k$-th discord if the distance to its
$k$-th nearest non-self match is the largest.

Given the positive real number $r$, the discord at a distance at least $r$
from its nearest non-self match is called the _range discord_. That is, the
range discord $D$ with respect to the parameter $r$ meets the following:
$\min\bigl{(}\operatorname{Dist}(D,\,M_{D})\bigr{)}\geq r$.

The MERLIN [36] and DRAG [51] algorithms deal with subsequences of the time
series that previously z-normalized to have mean zero and a standard deviation
of one. Here, z-normalization of a subsequence $X=\\{x\\}_{i=1}^{m}\in
S^{m}_{T}$ is defined as a subsequence $\hat{X}=\\{\hat{x}\\}_{i=1}^{m}$,
where

$\begin{gathered}\hat{x}_{i}=\frac{x_{i}-\mu_{X}}{\sigma_{X}},\quad\mu_{X}=\frac{1}{m}\sum\limits_{i=1}^{m}x_{i},\quad\sigma^{2}_{X}={\frac{1}{m}\sum\limits_{i=1}^{m}x_{i}^{2}-\mu^{2}}.\end{gathered}$
(4)

Both MERLIN and DRAG algorithms employ the Euclidean metric as the
$\operatorname{Dist}(\cdotp,\cdotp)$ function to measure the distance between
subsequences that is defined as follows. Let us have $X,Y\in S^{m}_{T}$, then
the Euclidean distance between the subsequences is calculated as below:

$\operatorname{ED}(X,\;Y)=\sqrt{\sum\limits_{i=1}^{m}(x_{i}-y_{i})^{2}}.$ (5)

In our study, being motivated by the highest possible performance of discord
discovery, we employ the square of the Euclidean metric as a distance
function. For the sake of simplicity, we denote by $\operatorname{ED_{norm}}$
the Euclidean distance between two z-normalized subsequences:
$\operatorname{ED_{norm}}(X,\,Y)=\operatorname{ED}(\hat{X},\,\hat{Y})$. To
compute $\operatorname{ED^{2}_{norm}}$, we further employ the following
technique proposed in [34] that allows for faster calculation than in Equation
5:

$\operatorname{ED^{2}_{norm}}(X,\,Y)=2m\left(1-\dfrac{X\cdotp
Y-m\cdotp\mu_{X}\cdotp\mu_{Y}}{m\cdotp\sigma_{X}\cdotp\sigma_{Y}}\right),$ (6)

where $X\cdotp Y$ denotes the scalar product of vectors
$X,Y\in\mathds{R}^{m}$.

###  2.2 MERLIN and DRAG Algorithms

Algorithm 1 MERLIN (in $T$, minL, maxL, topK; out $\mathcal{D}$)

1:$\mathcal{D}\leftarrow\varnothing$; $r\leftarrow 2\sqrt{\textit{minL}}$;
$\textit{nnDist}_{\textit{minL}}\leftarrow-\infty$

2:while $\textit{nnDist}_{\textit{minL}}<0$ and
$|D_{\textit{minL}}|<\textit{topK}$ do

3: $D_{minL}\leftarrow\hyperref@@ii[alg:DRAG]{\textsc{DRAG}}(T,minL,r)$;
$\mathcal{D}\leftarrow\mathcal{D}\cup D_{minL}$;
$\textit{nnDist}_{\textit{minL}}\leftarrow\min\limits_{d\in
D_{\textit{minL}}}d.\textit{nnDist}$

4: $r\leftarrow 0.5\cdot r$

5:for $i\leftarrow\textit{minL}+1$ to $\textit{minL}+4$ do

6: $\textit{nnDist}_{i}\leftarrow-\infty$

7: while $\textit{nnDist}_{i}<0$ and $|D_{i}|<\textit{topK}$ do

8: $r\leftarrow 0.99\cdot nnDist_{i-1}$

9: $D_{i}\leftarrow\hyperref@@ii[alg:DRAG]{\textsc{DRAG}}(T,i,r)$;
$\mathcal{D}\leftarrow\mathcal{D}\cup D_{i}$;
$\textit{nnDist}_{i}\leftarrow\min\limits_{d\in D_{i}}d.\textit{nnDist}$

10: $r\leftarrow 0.99\cdot r$

11:for $i\leftarrow\textit{minL}+5$ to maxL do

12:
$\mu\leftarrow\operatorname{Mean}(\\{\textit{nnDist}_{k}\\}^{i-5}_{k=i-1})$;
$\sigma\leftarrow\operatorname{Std}(\\{\textit{nnDist}_{k}\\}^{i-5}_{k=i-1})$;
$r\leftarrow\mu-2\sigma$

13: $D_{i}\leftarrow\hyperref@@ii[alg:DRAG]{\textsc{DRAG}}(T,i,r)$;
$\mathcal{D}\leftarrow\mathcal{D}\cup D_{i}$;
$\textit{nnDist}_{i}\leftarrow\min\limits_{d\in D_{i}}d.\textit{nnDist}$

14: while $nnDist_{i}<0$ and $|D_{i}|<\textit{topK}$ do

15: $D_{i}\leftarrow\hyperref@@ii[alg:DRAG]{\textsc{DRAG}}(T,i,r)$;
$\mathcal{D}\leftarrow\mathcal{D}\cup D_{i}$;
$\textit{nnDist}_{i}\leftarrow\min\limits_{d\in D_{i}}d.\textit{nnDist}$

16: $r\leftarrow r-\sigma$

17:return $\mathcal{D}$

Alg. 1 depicts a pseudo code of MERLIN [36] (up to discovering top-$k$
discords of each length instead of all ones in the specified length range).
Hereinafter, let us have an $n$-length time series $T$, and we are to find a
set $\mathcal{D}$ of its discords that have a length in the range
$\textit{minL}..\textit{maxL}$ (where $\textit{minL}\leq\textit{maxL}\ll n$),
so that $\mathcal{D}=\cup_{m=\textit{minL}}^{\textit{maxL}}D_{m}$, where
$D_{m}$ denotes a subset of $m$-length discords. The distance to the
$d\in\mathcal{D}$ discord’s nearest neighbor is denoted by
$d.\textit{nnDist}$.

The algorithm prescribes the following procedure to select the parameter $r$.
Discords are discovered sequentially, starting from a minimum length of the
specified discord range to a maximum one. At each step, MERLIN calculates the
arithmetic mean $\mu$ and the standard deviation $\sigma$ of the last five
distances from the discords found to their nearest neighbors, and then calls
the DRAG algorithm passing it the parameter $r=\mu-2\sigma$. If DRAG has not
found a discord, then $\sigma$ is subtracted from $r$ until DRAG stops
successfully (i.e., a discord will be found). For the first five discord
lengths, the parameter $r$ is set as follows. For discords of minimum length
minL, the parameter is set as $r=2\sqrt{\textit{minL}}$ since it is the
maximum possible distance between any pair of minL-length subsequences, and
then $r$ is reduced by half until DRAG with such a parameter results in
success. To obtain the next four discord lengths, the algorithm takes the
distance from the discord to its nearest neighbor obtained in the previous
step, minus a small value equal to 1%. The subtraction of an additional 1%
proceeds until the discord discovery with such a parameter results in success.
For a detailed explanation of the above-described procedure, we refer the
reader to the original work [36].

Algorithm 2 DRAG (in $T$, $m$, $r$; out $\mathcal{D}$)

Phase 1. Select candidates

1:$\mathcal{C}\leftarrow\\{T_{1,\,m}\\}$

2:for all $s\in S^{m}_{T}\smallsetminus T_{1,\,m}$ do

3: $\textit{isCand}\leftarrow$ TRUE

4: for all $c\in\mathcal{C}\;\textbf{and}\;c\in M_{s}$ do

5: if $\hyperref@@ii[eq:Euclid]{\operatorname{ED}}(s,c)<r$ then

6: $\mathcal{C}\leftarrow\mathcal{C}\smallsetminus c$

7: $\textit{isCand}\leftarrow$ FALSE

8: if isCand then

9: $\mathcal{C}\leftarrow\mathcal{C}\cup s$

10:return $\mathcal{C}$

Phase 2. Refine discords  1:$\mathcal{D}\leftarrow\varnothing$; $\forall
c\in\mathcal{C}\;c.\textit{nnDist}\leftarrow+\infty$ 2:for all $s\in
S^{m}_{T}$ do 3: for all $c\in\mathcal{C}\;\textbf{and}\;c\in
M_{s}\;\textbf{where}\;s\neq c$ do 4:
$dist\leftarrow\operatorname{EarlyAbandon\hyperref@@ii[eq:Euclid]{\operatorname{ED}}}(s,c)$
5: if $dist<r$ then 6: $\mathcal{C}\leftarrow\mathcal{C}\smallsetminus c$ 7:
else 8: $\mathcal{D}\leftarrow\mathcal{D}\cup c$ 9:
$c.\textit{nnDist}\leftarrow\operatorname{min}(c.\textit{nnDist},\textit{dist})$
10:return $\mathcal{D}$

The DRAG algorithm [51] (see Alg. 2) performs in two phases, namely the
candidate selection and discord refinement, where it collects potential range
discords and discards false positives, respectively. At the first phase, DRAG
scans through the time series $T$, and for each subsequence $s\in S^{m}_{T}$
it validates the possibility for each candidate $c$ already in the candidate
set $\mathcal{C}$ to be discord. If a candidate $c$ fails the validation, then
it is removed from this set. In the end, the new $s$ is either added to the
candidates set, if it is likely to be a discord, or it is pruned. At the
second phase, the algorithm initially sets distances of all candidates to
their nearest neighbors to positive infinity. Then, DRAG scans through the
time series $T$, calculating the distance between each subsequence $s\in
S^{m}_{T}$ and each candidate $c$. When calculating $\operatorname{ED}(s,c)$,
the $\operatorname{EarlyAbandonED}$ procedure stops the summation of
$\sum^{m}_{k=1}(s_{k}-c_{k})^{2}$ if it reaches $k=\ell$, such that
$1\leq\ell\leq m$ for which $\sum^{\ell}_{k=1}(s_{k}-c_{k})^{2}\geq
c.nnDist^{2}$. If the distance is less than $r$ then the candidate is false
positive and permanently removed from $\mathcal{C}$. If the above-mentioned
distance is less than the current value of $c.nnDist$ (and still greater than
$r$, otherwise it would have been removed) then the current distance to the
nearest neighbor is updated. The correctness of the above-described procedure
is proved in the original work [51].

##  3 Arbitrary Length Discords Discovery with GPU

Currently, GPU (Graphics Processing Unit) [25] is one of the most popular
many-core hardware platforms. GPU fits well for SIMD (Single Instructions
Multiple Data) computations being composed of symmetric streaming
multiprocessors, each of wich, in turn, consists of symmetric CUDA (Compute
Unified Device Architecture) cores. CUDA API (Application Programming
Interface) makes it possible to assign multiple threads to execute the same
set of instructions over multiple data. In CUDA, all threads form a _grid_
consisting of blocks. In a _block_ , threads are divided into _warps_ ,
logical groups of 32 threads. The block’s threads run in parallel and
communicate with each other through shared memory. A CUDA function is called a
_kernel_. When run a kernel on GPU, an application programmer specifies both
the number of blocks in the grid and the number of threads in each block.

Below, in Section 3.1, we introduce general architecture and data structures
of PALMAD, the parallel algorithm to discover discords of every possible
length in a specified range on GPU that is based on the original serial MERLIN
algorithm [36]. PALMAD employs PD3 (Parallel DRAG-based Discord Discovery)
[28], our parallel version of the original serial DRAG algorithm [51].
Similarly to the original serial algorithm, PD3 performs in two phases, where
each phase is parallelized separately from the other. In Sections 3.2 and 3.3,
we discuss the parallelization of the candidate selection and discords
refinement phases, respectively.

###  3.1 General Architecture and Data Structures

####  3.1.1 Avoiding Redundant Calculations

Basically, our parallel algorithm follows the computational scheme of its
serial predecessor (see Alg. 1). PALMAD employs repeated calls of PD3, our
designed parallel version of the serial DRAG algorithm (see Alg. 2 and lines
3, 9, 13, and 15 in Alg. 1) for a graphics processor. As opposed to the
original algorithm, PALMAD avoids redundant calculations in iterative calls of
DRAG when the subsequence length is one more than at the previous step.
Indeed, to calculate the distance between any two candidate subsequences, we
need to partially repeat alike calculations regarding subsequences that are
one less length self-matches to the candidates above. More formally, for any
$i,\,j$ ($1<i\leq n-m$ and $3\leq m\ll n$), when calculating
$\operatorname{ED^{2}_{norm}}(T_{i,\,m},T_{j,\,m})$ through Equation 6 from
scratch, we partially repeat calculations of both mean values
$\mu_{T_{i,\,m-1}}$ and $\mu_{T_{i,\,m-1}}$, and standard deviations
$\sigma_{T_{i,\,m-1}}$ and $\sigma_{T_{i,\,m-1}}$ through Equation 4.

To avoid the overhead above, we employ the vectors
$\bar{\mu},\bar{\sigma}\in\mathds{R}^{n-minL+1}$, namely the mean values and
standard deviations of all the given time series subsequences of the given
length, respectively. In these vectors, first $n-m+1$ elements are processed,
where $m$ is the given subsequence length ($\textit{minL}\leq
m\leq\textit{maxL}$), and the rest ones are left unattended.

For the minL-length subsequences, these vectors are calculated according to
Equation 4 once before the very first call of PD3 (see line 3 in Alg. 1)
whereas for the rest values of $m$, the vectors $\bar{\mu}$ and $\bar{\sigma}$
are updated before the each further call of PD3 (see lines 9, 13, and 15 in
Alg. 1) according to the following recurrent formulas:

$\displaystyle\mu_{T_{i,\,m+1}}$ $\displaystyle=$
$\displaystyle\dfrac{1}{m+1}\bigl{(}m\mu_{T_{i,\,m}}+t_{i+m}\bigr{)},$ (7)
$\displaystyle\sigma^{2}_{T_{i,\,m+1}}$ $\displaystyle=$
$\displaystyle\dfrac{m}{m+1}\Bigl{(}\sigma_{T_{i,\,m}}^{2}+\dfrac{1}{m+1}\bigl{(}\mu_{T_{i,\,m}}-t_{i+m}\bigr{)}^{2}\Bigr{)}.$
(8)

In order not to overload the article’s text with details, we present the lemma
with the proof of Equations 7 and 8 in Appendix.

Initial calculation and update of the vectors $\bar{\mu}$ and $\bar{\sigma}$
are implemented as CUDA kernels. We form a grid of $N$ threads where the
number of threads in each block is the algorithm’s parameter that is set as a
multiple of the GPU warp size. Each thread calculates elements of the vectors
$\bar{\mu}$ and $\bar{\sigma}$ according to Equations 7 and 8.

####  3.1.2 Data Structures and Segmentation

Figure 1: Data structures

In Fig. 1, we depict basic data structures involved in computations. The real-
valued vector $\textit{nnDist}\in\mathds{R}^{N}$, for a subsequence of the
input time series contains distance to its nearest neighbor:
$\textit{nnDist}(i)=\min\bigl{(}\operatorname{ED^{2}_{norm}}(T_{i,\,m},M_{T_{i,\,m}})\bigr{)}$.
Two boolean-valued vectors $\textit{Cand},\textit{Neighbor}\in\mathds{B}^{N}$
are the bitmaps for the subsequences and their nearest neighbors,
respectively: $\textit{Cand}(i)=\texttt{TRUE}$ (or
$\textit{Neighbor}(i)=\texttt{TRUE}$, respectively) if the subsequence
$T_{i,\,m}$ (or its nearest neighbor, respectively) is a discord, and FALSE
otherwise. These bitmaps are initialized with TRUE values. Further, we employ
element-wise conjunction of the bitmaps above to discard more candidates
during processing relying on the obvious fact that a subsequence that is not a
discord cannot have a nearest neighbor that is a discord.

Figure 2: Data segmentation

To implement the candidate selection and discord refinement phases, we exploit
the data parallelism concept and segment data as depicted in Fig. 2. The time
series is divided into equal-length segments, where each segment is processed
separately by a block of GPU threads. Performing the phase, the thread block
scans the subsequences in chunk-wise manner. The number of elements in a chunk
is equal to the segment length, and the first chunk begins with the $m$-th
element in the segment. Such a technique avoids redundant checks of the fact
that candidates and subsequences in chunks overlap.

The segment length is the algorithm’s parameter to be set as a multiple of the
GPU warp size. To balance the load of threads in the block, we require that
the number of $m$-length subsequences in the time series is a multiple of the
number of subsequences of the specified length in the segment. If this is not
the case, we pad the time series right with dummy positive infinity-valued
elements. Let us denote the segment length and the number of $m$-length
subsequences in the segment by seglen and segN, respectively, then
$\textit{segN}=\textit{seglen}-m+1$. Let the number of dummy elements in the
rightmost segment is denoted by pad, then it is defined as follows:

$pad=\begin{cases}m-1,&\text{$N$ \text{mod} $\textit{segN}$ = 0}\\\
\lceil\dfrac{N}{\textit{segN}}\rceil\cdotp\textit{segN}+2(m-1)-n,&\text{otherwise}\end{cases}.$
(9)

###  3.2 Parallelization of Candidate Selection

Algorithm 3 PD3select (in $T$, $m$, $r$; out $\mathcal{C}$)

1:$\textit{Cand}\leftarrow\overline{\texttt{TRUE}}$;
$\textit{Neighbor}\leftarrow\overline{\texttt{TRUE}}$

2:for all $T^{(i)}\in T$ do$\triangleright$ PARALLEL (block)

3: for all $\textit{Chunk}^{(j)}\in T^{(i)}$ where $i\leq j$
do$\triangleright$ PARALLEL (thread)

4: if $i=j$ then

5:
$\textit{QTrow}\leftarrow\textsc{CalcDotProducts}(T_{1,\,m}^{(i)},\textit{Chunk}^{(j)})$

6: continue

7:
$\textit{QTrow}\leftarrow\textsc{UpdateDotProducts}(\textit{QTrow},T_{1,\,m}^{(i)},\textit{Chunk}^{(j)})$

8:
$\textit{QTcol}\leftarrow\textsc{CalcDotProducts}(\textit{Chunk}_{1,\,m}^{(j)},T^{(i)})$

9:
$\textit{dist}\leftarrow\textsc{CalcDist}(\textit{Chunk}_{1,\,m}^{(j)},T^{(i)},\textit{QTcol},\bar{\mu},\bar{\sigma})$

10: if $\textit{dist}<r$ then

11:
$\textit{Cand}(i\cdotp\textit{segN}+tid)\leftarrow\texttt{FALSE};\textit{Neighbor}(j\cdotp\textit{segN}+1)\leftarrow\texttt{FALSE}$

12: else

13:
$\textit{nnDist}(j\cdot\textit{segN}+1)\leftarrow\min\bigl{(}\textit{dist},\textit{nnDist}(j\cdot\textit{segN}+1)\bigr{)}$

14: if not
$\bigvee_{k=i\cdotp\textit{segN}}^{(i+1)\cdotp\textit{segN}}\textit{Cand}(k)$
then

15: break

16: for all $\textit{Chunk}_{k,\,m}^{(j)}\in
S^{m}_{\textit{Chunk}^{(j)}}\smallsetminus\textit{Chunk}_{1,\,m}^{(j)}$
do$\triangleright$ PARALLEL (thread)

17:
$\textit{QTcol}\leftarrow\textsc{UpdateDotProducts}(\textit{QTcol},\textit{QTrow},\textit{Chunk}_{k,\,m}^{(j)},T^{(i)})$

18:
$\textit{dist}\leftarrow\textsc{CalcDist}(\textit{Chunk}_{k,\,m}^{(j)},T^{(i)},\textit{QTcol},\bar{\mu},\bar{\sigma})$

19: if $\textit{dist}<r$ then

20:
$\textit{Cand}(i\cdotp\textit{segN}+tid)\leftarrow\texttt{FALSE};\textit{Neighbor}(j\cdotp\textit{segN}+k)\leftarrow\texttt{FALSE}$

21: else

22:
$\textit{nnDist}(j\cdot\textit{segN}+1)\leftarrow\min\bigl{(}\textit{dist},\textit{nnDist}(j\cdot\textit{segN}+1)\bigr{)}$

23: if not
$\bigvee_{k=i\cdotp\textit{segN}}^{(i+1)\cdotp\textit{segN}}\textit{Cand}(k)$
then

24: break

25:$\mathcal{C}\leftarrow\bigl{\\{}\\{T_{i,\,m}\in
S^{m}_{T};\,\textit{nnDist}(i)\\}\mid 1\leq i\leq
n-m+1,\textit{Cand}(i)=\texttt{TRUE}\bigr{\\}}$

26:return $\mathcal{C}$

Alg. 3 depicts a pseudo code of our parallelization of the candidate selection
phase. The respective CUDA kernel forms a grid consisting of
$\lceil\tfrac{N}{\textit{segN}}\rceil$ blocks of segN threads in each block. A
thread block considers the segment subsequences as local candidates to
discords and performs chunk-wise processing of the subsequences that are
located to the _right_ of the segment and do not overlap with the candidates.
The processing of the subsequences is as follows. If the distance from the
candidate to the subsequence is less than the parameter $r$, then the
candidate and the subsequence are excluded from further processing as
obviously not discords (the corresponding flags in the bitmaps are set to
FALSE). If all the local candidates are discarded, the block terminates all
its threads ahead of schedule.

Figure 3: Computational kernel for the candidate selection phase

In Fig. 3, we show in detail how the block’s threads work. The thread block
loads its segment into the shared memory once before starting calculations,
and at each scanning step also loads there the current chunk located to the
right of the segment. This technique allows for increasing the algorithm
performance through reduction the number of reads the time series elements
from the global memory. Next, the block threads calculate scalar products,
storing the results in shared memory: firstly, products between the first
subsequence of the segment and all the subsequences of the current chunk, and
then these ones between the first subsequence of the current chunk and all the
subsequences of the segment (the vectors
$\textit{QTrow},\textit{QTcol}\in\mathds{R}^{\textit{segN}}$ in lines 4–7 and
8 in Alg. 3, respectively).

Further, based on the obtained vector QTcol and the pre-calculated vectors
$\bar{\mu},\bar{\sigma}$, we calculate the distances between the first
subsequence of the chunk and all the subsequences of the segment through the
Equation 6 (see line 9 in Alg. 3). Employing the calculated distances, we
discard unpromising candidates located in the segment and the current chunk
(see lines 10–11 in Alg. 3). If all the candidates in the segment are
discarded, then the block stops (see lines 14–15 in Alg. 3).

After that, the block threads perform similar actions over the remaining
subsequences of the current chunk, however, calculating scalar products more
efficiently (see lines 16–24 in Alg. 3). We calculate the scalar products
between the current subsequence of the chunk and all the subsequences of the
segment (i.e., the vector QTcol) based on the previously calculated vector
QTrow and the vector QTcol obtained at the previous iteration (see line 16 in
Alg. 3). To calculate the scalar product of the $k$-th
$(1<k\leq\textit{segN})$ subsequence in the $\textit{Chunk}^{(j)}$ and a
subsequence in the segment $T^{(i)}$, we employ the following formula:

$QTcol(tid)=\begin{cases}\textit{QTcol}(\textit{tid}-1)+T_{\textit{tid},\,m}^{(i)}\cdotp\textit{Chunk}_{k,\,m}^{(j)}(m)-\\\
\quad-
T_{\textit{tid}-1,\,m}^{(i)}(1)\cdotp\textit{Chunk}_{k-1,\,m}^{(j)}(1),&\text{$1<\textit{tid}\leq\textit{segN}$}\\\
\textit{QTrow}(k),&\text{$\textit{tid}=1$}\end{cases},$ (10)

where $tid$ denotes the number of a thread in the block. Since in Equation 10,
the first term is obtained from the previous iteration, we achieve the $O(1)$
complexity of the scalar product calculation instead of $O(m)$ as in the
straightforward case.

###  3.3 Parallelization of Discord Refinement

Algorithm 4 PD3refine (in $T$, $m$, $r$; out $\mathcal{D}$)

1:for all $T_{i,\,m}\in S^{m}_{T}$ do$\triangleright$ PARALLEL (thread)

2: $\textit{Cand}(i)\leftarrow\textit{Cand}(i)\wedge\textit{Neighbor}(i)$

3:for all $T^{(i)}\in T$ where
$\bigwedge_{k=i\cdotp\textit{segN}}^{(i+1)\cdotp\textit{segN}}\textit{Cand}(k)=\texttt{TRUE}$
do$\triangleright$ PARALLEL (block)

4: for all $Chunk^{(j)}\in T^{(i)}$ where $i\geq j$ do$\triangleright$
PARALLEL (thread)

5: if $i=j$ then

6:
$\textit{QTrow}\leftarrow\textsc{CalcDotProducts}(T_{1,\,m}^{(i)},\textit{Chunk}^{(j)})$

7: continue

8:
$\textit{QTrow}\leftarrow\textsc{UpdateDotProducts}(\textit{QTrow},T_{1,\,m}^{(i)},\textit{Chunk}^{(j)})$

9:
$\textit{QTcol}\leftarrow\textsc{CalcDotProducts}(\textit{Chunk}_{1,\,m}^{(j)},T^{(i)})$

10:
$dist\leftarrow\textsc{CalcDist}(\textit{Chunk}_{1,\,m}^{(j)},T^{(i)},\textit{QTcol},\bar{\mu},\bar{\sigma})$

11: if $\textit{dist}<r$ then

12: $\textit{Cand}(i\cdotp\textit{segN}+\textit{tid})\leftarrow\texttt{FALSE}$

13: else

14:
$\textit{nnDist}(j\cdot\textit{segN}+\textit{tid})\leftarrow\min\bigl{(}\textit{dist},\textit{nnDist}(j\cdot\textit{segN}+\textit{tid})\bigr{)}$

15: if not
$\bigvee_{k=i\cdotp\textit{segN}}^{(i+1)\cdotp\textit{segN}}\textit{Cand}(k)$
then

16: break

17: for all $\textit{Chunk}_{k,m}^{(j)}\in
S^{m}_{\textit{Chunk}^{(j)}}\smallsetminus\textit{Chunk}_{1,\,m}^{(j)}$
do$\triangleright$ PARALLEL (thread)

18:
$\textit{QTcol}\leftarrow\textsc{UpdateDotProducts}(\textit{QTcol},QTrow,Chunk_{k,\,m}^{(j)},T^{(i)})$

19:
$\textit{dist}\leftarrow\textsc{CalcDist}(\textit{Chunk}_{k,\,m}^{(j)},T^{(i)},\textit{QTcol},\bar{\mu},\bar{\sigma})$

20: if $dist<r$ then

21: $\textit{Cand}(i\cdotp\textit{segN}+\textit{tid})\leftarrow\texttt{FALSE}$

22: else

23: $\textit{nnDist}(j\cdot
segN+tid)\leftarrow\min\bigl{(}\textit{dist},\textit{nnDist}(j\cdot\textit{segN}+\textit{tid})\bigr{)}$

24: if not
$\bigvee_{k=i\cdotp\textit{segN}}^{(i+1)\cdotp\textit{segN}}\textit{Cand}(k)$
then

25: break

26:$\mathcal{D}\leftarrow\bigl{\\{}\\{T_{i,\,m}\in
S^{m}_{T};\,\textit{nnDist}(i)\\}\mid 1\leq i\leq
n-m+1,\textit{Cand}(i)=\texttt{TRUE}\bigr{\\}}$

27:return $\mathcal{D}$

The discord refinement phase (see Alg. 4) parallelization is implemented
through two CUDA kernels called one after the other. The first one trivially
refines discords obtained in the previous phase through the element-wise
conjunction of the $Cand$ and $Neighbor$ bitmap vectors, writing the result to
the former. This operation allows for pruning the nearest neighbors of the
subsequences discarded at the selection phase.

The second kernel performs non-trivial refinement and is parallelized similar
to the selection phase involving only those segments of the time series whose
set of local candidates is not empty (see line 3 in Alg. 4). The algorithm
scans and processes the subsequences that do not overlap with the candidates
and are located to the _left_ of the segment (see line 4 in Alg. 4). If the
distance between the candidate and the subsequence is less than the parameter
$r$, then the candidate is discarded as obvious false positive.

##  4 Experimental Evaluation

To evaluate the proposed algorithm, we carried out the experiments and study
the performance of PALMAD over various real-world and synthetic time series in
comparison with analogs, and investigated the algorithm’s scalability. We
designed the experiments to be easily reproducible with our repository [27]
that contains the algorithm’s source code and all the datasets used in this
work. Below, Section 4.1 describes hardware and time series employed in the
experiments, and Section 4.2 presents experimental results and discussion.

###  4.1 The Experimental Setup

Table 1: Time series employed in the experiments Time series |  | Length
---
$(n)$
| Discord length
---
$(\textit{minL}=\textit{maxL})$
Domain
Space shuttle | 50 000 | 150 |  | Measurements of a sensor on the NASA
---
spacecraft
ECG | 45 000 | 200 | Electrocardiogram of an adult patient
ECG-2 | 21 600 | 400
Koski-ECG | 100 000 | 458
Respiration | 24 125 | 250 | Human breathing by chest expansion
Power demand | 33 220 | 750 | Annual energy consumption of an office
RandomWalk1M | $10^{7}$ | 512 | Synthetic time series
RandomWalk2M | $2\cdot 10^{7}$ | 512

In our study, we employed the time series listed in Tab. 1 that are also used
in the experimental evaluation of HOTSAX [24], KBF_GPU [46], and Zhu et al.’s
[54] algorithm. The Space shuttle data [14] are solenoid current measurements
on a Marotta MPV-41 series valve as the valve is cycled on and off under
various test conditions in a laboratory where the valves are used to control
fuel flow on the NASA spacecraft. The ECG and ECG2 [17], and Koski-ECG [26]
time series are electrocardiograms of adult patients. The Respiration time
series [24] shows a patient’s breathing (measured by thorax extension), as
s/he wakes up. The Power demand time series reflects the energy consumption of
the research center in Netherlands for 1997 [49]. The RandomWalk1M and
RandomWalk2M are our generated time series through the random walk model [37].

Table 2: Hardware platform of the experiments Specifications | GPU-SUSU | GPU-MSU
---|---|---
Brand and product line | NVIDIA Tesla
Model | V100 | P100
# cores | 5 120 | 3 584
Core frequency, GHz | 1.3 | 1.19
Memory, Gb | 32 | 16
| Peak performance
---
(double precision), TFLOPS
7 | 4

Tab. 2 summarizes hardware platform of our experiments, where GPU-SUSU and
GPU-MSU denote graphics processors installed in the HPC centers of the South
Ural State University [12] and Moscow State University [47], respectively.

###  4.2 Results and Discussion

####  4.2.1 Comparison with Analogs

In the experiments, we compared PALMAD with two algorithms, namely, KBF_GPU
[46] and Zhu et al.’s [54] since our thorough review of related work (see
Section 1) did not reveal other GPU-oriented parallel competitors. Since the
authors of the above rivals do not provide their source codes, for a fair
comparison, in the experiments, we utilize time series and hardware identical
to those employed in [46] and [54], respectively, and compare our results with
ones reported in the original papers by the authors. The time series Koski-ECG
(see Tab. 1) was employed to compare PALMAD with KBF_GPU on the GPU-SUSU
system (see Tab. 2), and the rest time series were involved in comparison with
Zhu et al.’s algorithm on GPU-MSU. We omit the MERLIN performance results
since, as expected, the original serial algorithm is significantly inferior to
its parallel descendant, although in the experiments, we confirmed that PALMAD
produces exactly the same results as MERLIN. For each experiment, we ran
PALMAD 10 times and took the average value as the final running time.

Since the rival algorithms discover only the top-1 discord whereas our
algorithm finds all the discords of each length in a specified length range,
to provide a fair comparison, in the experiments, we employ two following
settings for PALMAD. First, we set $\textit{minL}=\textit{maxL}$ for the range
above. Second, we measure both running time of PALMAD and the number of
discords found to further show also the average time spent by our algorithm to
discover one discord.

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10^{0}$$10^{2}$Y1.0e0]pt5.3 (1)Y2.0e0]pt17
(1)Y3.0e0]pt36.6 (1)Y4.0e0]pt61 (1)Y5.0e0]pt90.8 (1)Y6.0e0]pt129
(1)Y7.0e0]pt172.9 (1)Y8.0e0]pt229.9 (1)Y9.0e0]pt280 (1)Y1.0e1]pt348.9
(1)Y1.0e0]pt0.4 (174)Y2.0e0]pt0.5 (63)Y3.0e0]pt1.18 (195)Y4.0e0]pt1.46
(221)Y5.0e0]pt1.73 (117)Y6.0e0]pt1.23 (286)Y7.0e0]pt0.68 (28)Y8.0e0]pt0.84
(28)Y9.0e0]pt1.11 (28)Y1.0e1]pt2.34 (644)Time series length, $\times 10^{4}$
Runtime, s (log scale) KBF_GPUPALMAD (# discords) (a) Overall running time

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10^{-2}$$10^{1}$$10^{4}$Y1.0e0]pt5.3Y2.0e0]pt17Y3.0e0]pt36.6Y4.0e0]pt61Y5.0e0]pt90.8Y6.0e0]pt129Y7.0e0]pt172.9Y8.0e0]pt229.9Y9.0e0]pt280Y1.0e1]pt348.9Y1.0e0]pt0.002Y2.0e0]pt0.008Y3.0e0]pt0.006Y4.0e0]pt0.007Y5.0e0]pt0.015Y6.0e0]pt0.004Y7.0e0]pt0.024Y8.0e0]pt0.03Y9.0e0]pt0.036Y1.0e1]pt0.004Time
series length, $\times 10^{4}$ Average runtime per discord, s (log scale)
KBF_GPUPALMAD (b) Average running time per one discord

Figure 4: Performance of PALMAD in comparison with KBF_GPU

In Fig. 4, we compare the performance of PALMAD and KBF_GPU. It can be seen
that our algorithm significantly outruns the rival in terms of both the
overall running time and the average running time to discover one discord.
Obviously, the reason is that KBF_GPU implements brute-force approach whereas
PALMAD avoids redundant calculations and exploits advanced data structures.

SpaceshuttleECGECG2PowerdemandRespirationRandomWalk1MRandomWalk2M$10^{1}$$10^{3}$$10^{5}$$10^{7}$Y0.0e0]pt7
(1)Y1.0e0]pt76 (1)Y2.0e0]pt26 (1)Y3.0e0]pt63 (1)Y4.0e0]pt26 (1)Y5.0e0]pt5 431
(1)Y6.0e0]pt15 190 (1)Y0.0e0]pt71 (963)Y1.0e0]pt1 056 (111)Y2.0e0]pt495 (1
618)Y3.0e0]pt353 (1 290)Y4.0e0]pt584 (7 915)Y5.0e0]pt394 838 (1 047
809)Y6.0e0]pt1 568 880 (2 096 385) Runtime, s (log scale) Zhu _et al._ PALMAD
(# discords) (a) Overall running time

SpaceshuttleECGECG2PowerdemandRespirationRandomWalk1MRandomWalk2M$10^{-1}$$10^{1}$$10^{3}$$10^{5}$$10^{7}$$10^{9}$6.575.5625.6163.3726.455
431.015 190.071.241 056.04495.45353.12584.4394 837.751 568
880.380.7495.143.062.740.743.777.480.338.061.231.090.291.512.990.1519.030.610.550.150.751.50.112.690.410.360.10.51.00.079.510.310.270.070.380.75
Average runtime per discord, ms (log scale) Zhu _et al._ PALMAD, topK =
1PALMAD, topK = 10% of $|\mathcal{D}|$PALMAD, topK = 25% of
$|\mathcal{D}|$PALMAD, $topK$ = 50% of $|\mathcal{D}|$PALMAD, topK = 75% of
$|\mathcal{D}|$PALMAD, topK = $|\mathcal{D}|$ (b) Average running time per one
discord

Figure 5: Performance of PALMAD in comparison with Zhu et al.’s algorithm

Figure 5 depicts the experimental results on the PALMAD performance comparing
with Zhu et al.’s algorithm. It can be seen that Zhu et al.’s algorithm
significantly outruns PALMAD: up to 20 times and up to two orders of magnitude
greater over real and synthetic time series, respectively. However, at the
same time, PALMAD discovers substantially more discords: at least to two and
seven orders of magnitude greater over real and synthetic time series,
respectively. Thus, comparing the average running time to discover one
discord, it can be seen that PALMAD significantly outruns the rival starting
at least from the moment when we set $topK$, the number of discords to be
discovered, as a quarter of actual number of discords found: at least two
times and three orders of magnitude over real and synthetic time series,
respectively.

####  4.2.2 Scalability of PALMAD

In addition to comparison of our algorithm with analogs, we also study the
scalability of PALMAD. First, we investigate the impact of the segment length
(the parameter $seglen$, see Section 3.2) on the PALMAD performance. Second,
we asses our algorithm’s performance depending on two input parameters that
directly affect the amount of calculations, namely the time series length and
discord range length.

SpaceshuttleRespirationECG2PowerdemandECG$200$$400$Y4.52e1]pt45.2Y2.372e2]pt237.2Y3.119e2]pt311.9Y3.96e2]pt396.0Y4.966e2]pt496.6Y3.69e1]pt36.9Y1.968e2]pt196.8Y2.33e2]pt233.0Y2.621e2]pt262.1Y4.197e2]pt419.7Y3.38e1]pt33.8Y1.799e2]pt179.9Y1.979e2]pt197.9Y1.986e2]pt198.6Y3.872e2]pt387.2
Average runtime per discord, ms seglen = 128seglen = 256seglen = 512

RandomWalk1MRandomWalk2M2$\cdot 10^{5}$3$\cdot 10^{5}$4$\cdot 10^{5}$6$\cdot
10^{5}$110 255.6430 711.7108 561.3428 368.0107 667.0425 177.3 Average runtime
per discord, ms (log scale)

Figure 6: Scalability of the PALMAD algorithm w.r.t. the segment length

In Fig. 6, we show experimental results regarding the impact of the segment
length on the PALMAD performance. It can be seen that the algorithm’s running
time is proportional to the segment length for both real-world and synthetic
time series, and the greater value of the segment length provides higher
performance. This can be explained by the fact that when the segment length
increases, the overhead of reading and writing segments in the GPU shared
memory decreases. Moreover, this was the reason that we took $seglen=512$ in
the above-described experiments.

$0$$5$$10$$25$$50$$75$$100$$100$$200$$300$$26.87$$52.7$$131.71$$210.82$$257.05$$324.24$Time
series length, $\times 10^{3}$ Runtime, s

$0$$5$$10$$25$$50$$75$$100$$0.2$$0.4$$0.6$$0.8$$0.06$$0.12$$0.29$$0.46$$0.56$$0.71$Time
series length, $\times 10^{3}$ Average runtime per discord, s

(a) Real dataset (Koski-ECG, discord range is 458..916)

$0$$5$$10$$25$$50$$75$$100$$500$$1\,000$$1\,500$$58.62$$139.9$$380.25$$731.8$$997.41$$1\,376.78$Time
series length, $\times 10^{4}$ Runtime, s

$0$$5$$10$$25$$50$$75$$100$$5$$10$$0.46$$1.09$$2.97$$5.72$$7.79$$10.76$Time
series length, $\times 10^{4}$ Average runtime per discord, s

(b) Synthetic dataset (RandomWalk1M, discord range is 128..256)

Figure 7: Scalability of the PALMAD algorithm w.r.t. the time series length

$2$$3$$4$$5$$6$$7$$8$$9$$10$$2\,000$$4\,000$$325.78$$1\,084.58$$1\,990.53$$3\,004.71$$4\,516.17$Discord
length range, $\times 458$ Runtime, s

$2$$3$$4$$5$$6$$7$$8$$9$$10$$0.8$$1$$0.71$$0.79$$0.87$$0.94$$1.1$Discord
length range, $\times 458$ Average runtime per discord, s

(a) Real dataset (Koski-ECG)

$2$$3$$4$$5$$6$$7$$8$$9$$10$$5\,000$$10\,000$$15\,000$$1\,376.89$$4\,365.74$$7\,446.32$$10\,652.79$$14\,499.2$Discord
length range, $\times 128$ Runtime, s

$2$$3$$4$$5$$6$$7$$8$$9$$10$$11$$11.5$$12$$12.5$$10.76$$11.37$$11.63$$11.89$$12.59$Discord
length range, $\times 128$ Average runtime per discord, s

(b) Synthetic dataset (RandomWalk1M)

Figure 8: Scalability of the PALMAD algorithm w.r.t. the discord length range

In Fig. 7 and Fig. 8, we depict the performance of our algorithm depending on
the time series length and on the discord length range, respectively, for the
cases of real-world and synthetic data. It can be observed that the
algorithm’s running time is proportional to the above-mentioned parameters for
both real-world and synthetic time series.

##  5 Case Study

In this section, we apply PALMAD to discover subsequence anomalies in a real-
world time series from a smart heating control system. The PolyTER system [64]
allows for intelligent monitoring and control of operating conditions of
utility systems through the analysis of the data from various IoT sensors
installed in the university campus buildings. We took a time series from a
temperature sensor installed in a lecture hall and discovered the anomalies in
a specified range. The sensor’s frequency is 4 times per hour, the time series
corresponds to annual measurements (i.e., time series length
$n=35\leavevmode\nobreak\ 040$), and we search for anomalies that range from
12 hours to 7 days (i.e., $\textit{minL}=48$ and $\textit{maxL}=672$,
respectively).

To visualize the results obtained, we propose the discord heatmap technique
that illustrates the anomaly score through the intensity of a color, and is
somewhat like the motif heatmap [33]. Formally speaking, we plot a one-color
heatmap as a matrix of the size $(\textit{maxL}-\textit{minL}+1)$ $\times$
$(n-\textit{minL})$, where the intensity of a pixel $(m,i)$ shows the anomaly
score of the discord $T_{i,\,m}\in D_{m}$, and the pixel’s intensity is
calculated as a normalization of the discord’s distance to its nearest
neighbor:

$\textit{heatmap}(m,\,i)=\frac{T_{i,\,m}.\textit{nnDist}}{2m},$ (11)

where we employ the normalizing divisor $2m$ according to Equation 6.

Despite the fact that we have proposed a visual tool to explore and discover
multiple length discords, there is an open question for a practitioner, how we
can rank discords of different lengths and extract the most interesting ones.
There is a number of a discord attributes that can be suggested to take into
account, e.g., its length and index, a distance to its nearest neighbor, a
number of its self-matches, etc. However, in this study, we employ a straight-
forward approach that considers a discord’s interest as a normalized distance
to its nearest neighbor comparing such distances among discords having the
same index. Thus, the most interesting discord among the ones of different
lengths is selected as below:

$\arg\max\limits_{1\leq i\leq N}\max\limits_{\textit{minL}\leq
m\leq\textit{maxL}}\textit{heatmap}(m,\,i).$ (12)

Clearly, through Equation 12, we can select top-$k$ interesting discords. In
our repository [27], the reader can find plots of discord heatmaps and top
discords for all the real-world time series listed in Tab. 1 and these ones
from other subject domains.

(a) Time series and its discord heatmap

(b) Zooming the most interesting intervals of the discord heatmap

Top-1 discord, $m=60$
(15 hr)

Top-2 discord, $m=248$
(2 days, 14 hr)

Top-3 discord, $m=90$
(22 hr 30 min)

Top-4 discord, $m=266$
(2 days, 18 hr 30 min)

Top-5 discord, $m=168$
(1 day, 18 hr)

Top-6 discord, $m=183$
(1 day, 21 hr 45 min)

(c) Top-6 discords of different lengths

Figure 9: Case of PolyTER

Figure 9 summarizes the results of the case study showing the time series and
its discord heatmap, zooming of the heatmap intervals with most interesting
discords, and top-6 discords according to Equation 12 (see sub-figures 9(a),
9(b), and 9(c), respectively). In addition, in top discord plots, we indicate
both indoor and outdoor temperature, where the latter is obtained from the
open weather archive [1]. Highly likely, each of top three discords
illustrates a long-term malfunction in a temperature sensor that outputs the
same measurements during a considerable time period. Next, top-4 and top-5
discords show a short-term failure of the sensor. Finally, top-6 discord may
indicate the fact that the system (or its operator) chose an inefficient
heating mode of the lecture hall in the considered period.

## Conclusions

In this article, we addressed the problem of accelerating the time series
subsequence anomaly discovery on a graphics processor. Such an anomaly refers
to successive points in time whose collective behavior is abnormal, although
each observation individually does not necessarily deviates. Currently,
discovering subsequence anomalies in time series remains one of the most
topical research problems.

Among a large number of approaches to discovering subsequence anomalies, the
discord concept [31] is considered one of the best. A time series discord is
intuitively defined as a subsequence that is maximally far away from its
nearest neighbor. However, application of discords is reduced by sensitivity
to a user’s choice of the subsequence length. A brute-force discovering
discords of all the possible lengths and then selecting the best discords with
respect to some measure is clearly computationally prohibitive. Recently
introduced the MERLIN algorithm [36] discovers time series discords of every
possible length in a specified range, being ahead of competitors in terms of
accuracy and performance. MERLIN employs repeated calls of the DRAG algorithm
[51] that discovers discords of a given length with a distance of at least $r$
to their nearest neighbors, and adaptive selection of the parameter $r$.
However, to the best of our knowledge, no research has addressed the
accelerating MERLIN with any parallel hardware architecture.

In the article, based on Keogh et al.’s works [51, 36], we proposed a novel
parallelization scheme PALMAD (Parallel Arbitrary Length MERLIN-based Anomaly
Discovery) for a graphics processor. When basically following the original
serial algorithm, PALMAD, however, employs our derived recurrent formulas to
calculate the mean values and standard deviations of subsequences of the time
series. Since that data are further involved in calculations of the normalized
Euclidean distances between subsequences, eventually, we significantly reduce
the amount of calculations. Furthermore, PALMAD repeatedly calls PD3 (Parallel
DRAG-based Discord Discovery) [28], our developed parallel version of the
original DRAG algorithm. Similar to its predecessor, PD3 performs in two
phases. To implement the candidate selection phase, we exploit the data
parallelism by division the time series into equal-length segments, where each
segment is processed separately by a block of GPU threads. The thread block
considers the segment subsequences as local candidates to discords and
processes the subsequences that are located to the right of the segment and do
not overlap with the candidates. Next, the thread block scans the subsequences
in chunks, the number of elements in which is equal to the segment length, and
the first chunk begins with the $m$-th element in the segment, where $m$ is
the discord length. Such a technique allows us for avoiding redundant checks
of the fact that candidates and subsequences in chunks overlap. In PD3, the
candidate refinement phase is parallelized similar to the selection phase.
Refinement involves only those segments of the time series whose set of local
candidates is not empty. The algorithm scans and processes the subsequences
that do not overlap with the candidates and are located to the left of the
segment.

We carried out an extensive experimental evaluation of PALMAD over real-world
and synthetic time series. In the experiments, we compared our development
with two algorithms, namely, KBF_GPU [46] and Zhu et al.’s [54] since our
thorough review of related work did not reveal other GPU-oriented parallel
competitors. Both rivals aim at discovery top-1 discord, where the former is a
parallelization of the brute-force approach while the latter employs
computational patterns to reduce amount of calculations. As expected, in the
experiments, our algorithm significantly outruns KBF_GPU. Next, being adapted
to discover the most important discord of a specified length, Zhu et al.’s
algorithm significantly outruns PALMAD that discovers discords of every
possible length in a specified range: up to 20 times and up to two orders of
magnitude greater over real and synthetic time series, respectively. However,
PALMAD discovers substantially more discords: at least to two and seven orders
of magnitude greater over real and synthetic time series, respectively. Thus,
PALMAD significantly outruns the rivals in terms of the average running time
to discover one discord. Finally, in the experiments, we also investigated the
scalability of PALMAD and found that the algorithm’s running time is
proportional to each of the following parameters for both real-world and
synthetic time series: the segment length, time series length, and discord
range length.

We also apply PALMAD to discover anomalies in a real-world time series from a
smart heating control system employing our proposed discord heatmap technique
to illustrate the results.

Our further studies might elaborate on the following topics: (a) discords
discovery in a large time series that cannot be entirely placed in RAM with a
high-performance cluster based on GPU nodes, and (b) application of PALMAD in
a deep learning-based online time series anomaly detection.

## Acknowledgment

This work was financially supported by the Russian Science Foundation (grant
No. 23-21-00465). The research is carried out using the equipment of the
shared research facilities of HPC computing resources at Lomonosov Moscow
State University and supercomputer resources of the South Ural State
University.

## References

  * [1] Weather for 243 countries of the world: Chelyabinsk. https://rp5.ru/Weather_in_Chelyabinsk/ (2022), accessed: 2022-11-17
  * [2] Ahmad, S., Lavin, A., Purdy, S., Agha, Z.: Unsupervised real-time anomaly detection for streaming data. Neurocomputing 262, 134–147 (2017). 10.1016/j.neucom.2017.04.070
  * [3] Avogadro, P., Dominoni, M.A.: A fast algorithm for complex discord searches in time series: HOT SAX Time. Applied Intelligence (2022). 10.1007/s10489-021-02897-z
  * [4] Blázquez-García, A., Conde, A., Mori, U., Lozano, J.A.: A review on outlier/anomaly detection in time series data. ACM Comput. Surv. 54(3), 56:1–56:33 (2021). 10.1145/3444690
  * [5] Bu, Y., Leung, O.T., Fu, A.W., Keogh, E.J., Pei, J., Meshkin, S.: WAT: finding top-_k_ discords in time series database. In: Proceedings of the 7th SIAM International Conference on Data Mining, April 26-28, 2007, Minneapolis, Minnesota, USA. pp. 449–454. SIAM (2007). 10.1137/1.9781611972771.43
  * [6] Buu, H.T.Q., Anh, D.T.: Time series discord discovery based on _i_ SAX symbolic representation. In: 3rd International Conference on Knowledge and Systems Engineering, KSE 2011, Hanoi, Vietnam, October 14-17, 2011. pp. 11–18. IEEE Computer Society (2011). 10.1109/KSE.2011.11
  * [7] Chandola, V., Banerjee, A., Kumar, V.: Anomaly detection: A survey. ACM Comput. Surv. 41(3), 15:1–15:58 (2009). 10.1145/1541880.1541882
  * [8] Chandola, V., Cheboli, D., Kumar, V.: Detecting anomalies in a time series database. Retrieved from the University of Minnesota Digital Conservancy. https://hdl.handle.net/11299/215791 (2009), accessed: 2022-04-12
  * [9] Chau, P.M., Duc, B.M., Anh, D.T.: Discord detection in streaming time series with the support of R-tree. In: 2018 International Conference on Advanced Computing and Applications (ACOMP), 27-29 November 2018, Ho Chi Minh City, Vietnam. pp. 96–103 (2018). 10.1109/ACOMP.2018.00023
  * [10] Choi, K., Yi, J., Park, C., Yoon, S.: Deep learning for anomaly detection in time-series data: Review, analysis, and guidelines. IEEE Access 9, 120043–120065 (2021). 10.1109/ACCESS.2021.3107975
  * [11] Dean, J., Ghemawat, S.: MapReduce: simplified data processing on large clusters. In: Brewer, E.A., Chen, P. (eds.) 6th Symposium on Operating System Design and Implementation (OSDI 2004), San Francisco, California, USA, December 6-8, 2004. pp. 137–150. USENIX Association (2004), http://www.usenix.org/events/osdi04/tech/dean.html
  * [12] Dolganina, N., Ivanova, E., Bilenko, R., Rekachinsky, A.: HPC resources of South Ural State University. In: Sokolinsky, L., Zymbler, M. (eds.) 16th International Conference on Parallel Computational Technologies, PCT 2022, Dubna, Russia, March 29-41, 2022, Revised Selected Papers. Communications in Computer and Information Science. vol. 1618, pp. 43–55. Springer (2022). 10.1007/978-3-031-11623-0_4
  * [13] Farahani, I.V., Chien, A., King, R.E., Kay, M.G., Klenz, B.: Time series anomaly detection from a Markov chain perspective. In: Wani, M.A., Khoshgoftaar, T.M., Wang, D., Wang, H., Seliya, N. (eds.) 18th IEEE International Conference On Machine Learning and Applications, ICMLA 2019, Boca Raton, FL, USA, December 16-19, 2019. pp. 1000–1007. IEEE (2019). 10.1109/ICMLA.2019.00170
  * [14] Ferrell, B., Santuro, S.: NASA shuttle valve data. http://www.cs.fit.edu/~pkc/nasa/data/ (2005), accessed: 2022-03-19
  * [15] Fredkin, E.: Trie memory. Commun. ACM 3(9), 490–499 (1960). 10.1145/367390.367400
  * [16] Fu, A.W., Leung, O.T., Keogh, E.J., Lin, J.: Finding time series discords based on Haar transform. In: Li, X., Zaïane, O.R., Li, Z. (eds.) Advanced Data Mining and Applications, Second International Conference, ADMA 2006, Xi’an, China, August 14-16, 2006, Proceedings. Lecture Notes in Computer Science, vol. 4093, pp. 31–41. Springer (2006). 10.1007/11811305_3
  * [17] Goldberger, A.L., Amaral, L.A.N., Glass, L., Hausdorff, J.M., Ivanov, P.C., Mark, R.G., Mietus, J.E., Moody, G.B., Peng, C.K., Stanley, H.E.: PhysioBank, PhysioToolkit, and PhysioNet components of a new research resource for complex physiologic signals. Circulation 101(23), 215–220 (2000). 10.1161/01.CIR.101.23.e215
  * [18] Huang, T., Zhu, Y., Mao, Y., Li, X., Liu, M., Wu, Y., Ha, Y., Dobbie, G.: Parallel discord discovery. In: Bailey, J., Khan, L., Washio, T., Dobbie, G., Huang, J.Z., Wang, R. (eds.) Advances in Knowledge Discovery and Data Mining \- 20th Pacific-Asia Conference, PAKDD 2016, Auckland, New Zealand, April 19-22, 2016, Proceedings, Part II. Lecture Notes in Computer Science, vol. 9652, pp. 233–244. Springer (2016). 10.1007/978-3-319-31750-2_19
  * [19] Huang, T., Zhu, Y., Wu, Y., Shi, W.: J-distance discord: An improved time series discord definition and discovery method. In: 2015 IEEE International Conference on Data Mining Workshop (ICDMW). pp. 303–310 (2015). 10.1109/ICDMW.2015.120
  * [20] Hundman, K., Constantinou, V., Laporte, C., Colwell, I., Söderström, T.: Detecting spacecraft anomalies using LSTMs and nonparametric dynamic thresholding. In: Guo, Y., Farooq, F. (eds.) Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD 2018, London, UK, August 19-23, 2018. pp. 387–395. ACM (2018). 10.1145/3219819.3219845
  * [21] Imani, S., Keogh, E.J.: Matrix profile XIX: time series semantic motifs: A new primitive for finding higher-level structure in time series. In: Wang, J., Shim, K., Wu, X. (eds.) 2019 IEEE International Conference on Data Mining, ICDM 2019, Beijing, China, November 8-11, 2019. pp. 329–338. IEEE (2019). 10.1109/ICDM.2019.00043
  * [22] Imani, S., Madrid, F., Ding, W., Crouter, S.E., Keogh, E.J.: Introducing time series snippets: a new primitive for summarizing long time series. Data Min. Knowl. Discov. 34(6), 1713–1743 (2020). 10.1007/s10618-020-00702-y
  * [23] Keogh, E.J., Chakrabarti, K., Pazzani, M.J., Mehrotra, S.: Dimensionality reduction for fast similarity search in large time series databases. Knowl. Inf. Syst. 3(3), 263–286 (2001). 10.1007/PL00011669
  * [24] Keogh, E.J., Lin, J., Fu, A.W.: HOT SAX: efficiently finding the most unusual time series subsequence. In: Proceedings of the 5th IEEE International Conference on Data Mining (ICDM 2005), 27-30 November 2005, Houston, Texas, USA. pp. 226–233. IEEE Computer Society (2005). 10.1109/ICDM.2005.79
  * [25] Kirk, D.B.: NVIDIA CUDA software and GPU parallel computing architecture. In: Morrisett, G., Sagiv, M. (eds.) Proceedings of the 6th International Symposium on Memory Management, ISMM 2007, Montreal, Quebec, Canada, October 21-22, 2007. pp. 103–104. ACM (2007). 10.1145/1296907.1296909
  * [26] Koski, A.: Primitive coding of structural ECG features. Pattern Recognit. Lett. 17(11), 1215–1222 (1996). 10.1016/0167-8655(96)00079-7
  * [27] Kraeva, Y., Zymbler, M.: PALMAD: Parallel MERLIN-based Anomaly Discovery algorithm for GPU. https://github.com/kraevaya/PALMAD (2022), accessed: 2022-12-01
  * [28] Kraeva, Y., Zymbler, M.: Parallel algorithm for time series discord discovery on a graphics processor. Pattern Recognition and Image Analysis 33(2), accepted for publication (2023)
  * [29] Kraeva, Y., Zymbler, M.L.: Scalable algorithm for subsequence similarity search in very large time series data on cluster of Phi KNL. In: Manolopoulos, Y., Stupnikov, S.A. (eds.) Data Analytics and Management in Data Intensive Domains - 20th International Conference, DAMDID/RCDL 2018, Moscow, Russia, October 9-12, 2018, Revised Selected Papers. Communications in Computer and Information Science, vol. 1003, pp. 149–164. Springer (2018). 10.1007/978-3-030-23584-0_9
  * [30] Li, G., Bräysy, O., Jiang, L., Wu, Z., Wang, Y.: Finding time series discord based on bit representation clustering. Knowl.-Based Syst. 54, 243–254 (2013). 10.1016/j.knosys.2013.09.015
  * [31] Lin, J., Keogh, E.J., Fu, A.W., Herle, H.V.: Approximations to magic: Finding unusual medical time series. In: 18th IEEE Symposium on Computer-Based Medical Systems (CBMS 2005), 23-24 June 2005, Dublin, Ireland. pp. 329–334. IEEE Computer Society (2005). 10.1109/CBMS.2005.34
  * [32] Lin, J., Keogh, E.J., Lonardi, S., Chiu, B.Y.: A symbolic representation of time series, with implications for streaming algorithms. In: Zaki, M.J., Aggarwal, C.C. (eds.) Proceedings of the 8th ACM SIGMOD workshop on Research issues in data mining and knowledge discovery, DMKD 2003, San Diego, California, USA, June 13, 2003. pp. 2–11. ACM (2003). 10.1145/882082.882086
  * [33] Madrid, F., Imani, S., Mercer, R., Zimmerman, Z., Senobari, N.S., Keogh, E.J.: Matrix profile XX: Finding and visualizing time series motifs of all lengths using the matrix profile. In: Gao, Y., Möller, R., Wu, X., Kotagiri, R. (eds.) 2019 IEEE International Conference on Big Knowledge, ICBK 2019, Beijing, China, November 10-11, 2019. pp. 175–182. IEEE (2019), https://doi.org/10.1109/ICBK.2019.00031
  * [34] Mueen, A., Nath, S., Liu, J.: Fast approximate correlation for massive time-series data. In: Elmagarmid, A.K., Agrawal, D. (eds.) Proceedings of the ACM SIGMOD International Conference on Management of Data, SIGMOD 2010, Indianapolis, Indiana, USA, June 6-10, 2010. pp. 171–182. ACM (2010). 10.1145/1807167.1807188
  * [35] Munir, M., Siddiqui, S.A., Dengel, A., Ahmed, S.: DeepAnT: A deep learning approach for unsupervised anomaly detection in time series. IEEE Access 7, 1991–2005 (2019). 10.1109/ACCESS.2018.2886457
  * [36] Nakamura, T., Imamura, M., Mercer, R., Keogh, E.J.: MERLIN: parameter-free discovery of arbitrary length anomalies in massive time series archives. In: Plant, C., Wang, H., Cuzzocrea, A., Zaniolo, C., Wu, X. (eds.) 20th IEEE International Conference on Data Mining, ICDM 2020, Sorrento, Italy, November 17-20, 2020. pp. 1190–1195. IEEE (2020). 10.1109/ICDM50108.2020.00147
  * [37] Pearson, K.: The problem of the random walk. Nature 72(1865), 294 (1905). 10.1038/072342a0
  * [38] Pfeilschifter, G.: Time series analysis with matrix profile on HPC systems. Master thesis, Department of Informatics, Technical University of Munich, Germany (2019), http://mediatum.ub.tum.de/doc/1471292/1471292.pdf, accessed: 2022-03-17
  * [39] Reyes, R., López-Rodríguez, I., Fumero, J.J., de Sande, F.: A preliminary evaluation of OpenACC implementations. J. Supercomput. 65(3), 1063–1075 (2013). 10.1007/s11227-012-0853-z
  * [40] Senin, P., Lin, J., Wang, X., Oates, T., Gandhi, S., Boedihardjo, A.P., Chen, C., Frankenstein, S.: Time series anomaly discovery with grammar-based compression. In: Alonso, G., Geerts, F., Popa, L., Barceló, P., Teubner, J., Ugarte, M., den Bussche, J.V., Paredaens, J. (eds.) Proceedings of the 18th International Conference on Extending Database Technology, EDBT 2015, Brussels, Belgium, March 23-27, 2015. pp. 481–492. OpenProceedings.org (2015). 10.5441/002/edbt.2015.42
  * [41] Shieh, J., Keogh, E.J.: _i_ SAX: indexing and mining terabyte sized time series. In: Li, Y., Liu, B., Sarawagi, S. (eds.) Proceedings of the 14th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, Las Vegas, Nevada, USA, August 24-27, 2008. pp. 623–631. ACM (2008). 10.1145/1401890.1401966
  * [42] Snir, M.: Technical perspective: The future of MPI. Commun. ACM 61(10), 105 (2018). 10.1145/3264415
  * [43] Son, N.T.: An improvement of disk aware discord discovery algorithm for discovering time series discord. In: 2020 5th International Conference on Green Technology and Sustainable Development (GTSD), November 27-28, 2020, Ho Chi Minh City, Vietnam. pp. 19–23 (2020). 10.1109/GTSD50082.2020.9303111
  * [44] de Supinski, B.R., Scogland, T.R.W., Duran, A., Klemm, M., Bellido, S.M., Olivier, S.L., Terboven, C., Mattson, T.G.: The ongoing evolution of OpenMP. Proc. IEEE 106(11), 2004–2019 (2018). 10.1109/JPROC.2018.2853600
  * [45] Thuy, H.T.T., Anh, D.T., Chau, T.N.V.: An effective and efficient hash-based algorithm for time series discord discovery. In: 2016 3rd National Foundation for Science and Technology Development Conference on Information and Computer Science (NICS), 14-16 September 2016, Danang, Vietnam. pp. 85–90 (2016). 10.1109/NICS.2016.7725673
  * [46] Thuy, T.T.H., Anh, T.D., Chau, T.N.V.: A new discord definition and an efficient time series discord detection method using GPUs. In: ICSED 2021: 2021 3rd International Conference on Software Engineering and Development, November 19-21, 2021, Xiamen, China. pp. 63–70 (2021). 10.1145/3507473.3507483
  * [47] Voevodin, V.V., Antonov, A.S., Nikitenko, D.A., Shvets, P.A., Sobolev, S.I., Sidorov, I.Y., Stefanov, K.S., Voevodin, V.V., Zhumatiy, S.A.: Supercomputer Lomonosov-2: large scale, deep monitoring and fine analytics for the user community. Supercomput. Front. Innov. 6(2), 4–11 (2019). 10.14529/jsfi190201
  * [48] Wei, L., Keogh, E., Xi, X.: SAXually explicit images: Finding unusual shapes. In: 6th International Conference on Data Mining (ICDM’06). pp. 711–720 (2006). 10.1109/ICDM.2006.138
  * [49] van Wijk, J.J., van Selow, E.R.: Cluster and calendar based visualization of time series data. In: IEEE Symposium on Information Visualization 1999 (INFOVIS’99), San Francisco, California, USA, October 24-29, 1999. pp. 4–9. IEEE Computer Society (1999). 10.1109/INFVIS.1999.801851
  * [50] Wu, Y., Zhu, Y., Huang, T., Li, X., Liu, X., Liu, M.: Distributed discord discovery: Spark based anomaly detection in time series. In: 17th IEEE International Conference on High Performance Computing and Communications, HPCC 2015, 7th IEEE International Symposium on Cyberspace Safety and Security, CSS 2015, and 12th IEEE International Conference on Embedded Software and Systems, ICESS 2015, New York, NY, USA, August 24-26, 2015. pp. 154–159. IEEE (2015). 10.1109/HPCC-CSS-ICESS.2015.228
  * [51] Yankov, D., Keogh, E.J., Rebbapragada, U.: Disk aware discord discovery: Finding unusual time series in terabyte sized datasets. In: Proceedings of the 7th IEEE International Conference on Data Mining (ICDM 2007), October 28-31, 2007, Omaha, Nebraska, USA. pp. 381–390. IEEE Computer Society (2007). 10.1109/ICDM.2007.61
  * [52] Yankov, D., Keogh, E.J., Rebbapragada, U.: Disk aware discord discovery: Finding unusual time series in terabyte sized datasets. Knowl. Inf. Syst. 17(2), 241–262 (2008). 10.1007/s10115-008-0131-9
  * [53] Yeh, C.M., Zhu, Y., Ulanova, L., Begum, N., Ding, Y., Dau, H.A., Zimmerman, Z., Silva, D.F., Mueen, A., Keogh, E.J.: Time series joins, motifs, discords and shapelets: A unifying view that exploits the matrix profile. Data Min. Knowl. Discov. 32(1), 83–123 (2018). 10.1007/s10618-017-0519-9
  * [54] Zhu, B., Jiang, Y., Gu, M., Deng, Y.: A GPU acceleration framework for motif and discord based pattern mining. IEEE Trans. Parallel Distributed Syst. 32(8), 1987–2004 (2021). 10.1109/TPDS.2021.3055765
  * [55] Zhu, Y., Imamura, M., Nikovski, D., Keogh, E.J.: Introducing time series chains: a new primitive for time series data mining. Knowl. Inf. Syst. 60(2), 1135–1161 (2019). 10.1007/s10115-018-1224-8
  * [56] Zhu, Y., Yeh, C.M., Zimmerman, Z., Kamgar, K., Keogh, E.J.: Matrix profile XI: SCRIMP++: time series motif discovery at interactive speeds. In: IEEE International Conference on Data Mining, ICDM 2018, Singapore, November 17-20, 2018. pp. 837–846. IEEE Computer Society (2018). 10.1109/ICDM.2018.00099
  * [57] Zimmerman, Z., Kamgar, K., Senobari, N.S., Crites, B., Funning, G.J., Brisk, P., Keogh, E.J.: Matrix profile XIV: scaling time series motif discovery with GPUs to break a quintillion pairwise comparisons a day and beyond. In: Proceedings of the ACM Symposium on Cloud Computing, SoCC 2019, Santa Cruz, CA, USA, November 20-23, 2019. pp. 74–86. ACM (2019). 10.1145/3357223.3362721
  * [58] Zymbler, M.: A parallel discord discovery algorithm for time series on many-core accelerators. Numerical Methods and Programming 20(3), 211–223 (2019). 0.26089/NumMet.v20r320, (in Russian)
  * [59] Zymbler, M., Goglachev, A.: Fast summarization of long time series with graphics processor. Mathematics 10(10), 1781 (2022). 10.3390/math10101781
  * [60] Zymbler, M., Grents, A., Kraeva, Y., Kumar, S.: A parallel approach to discords discovery in massive time series data. Computers, Materials & Continua 66(2), 1867–1878 (2021). 10.32604/cmc.2020.014232
  * [61] Zymbler, M., Ivanova, E.: Matrix profile-based approach to industrial sensor data analysis inside RDBMS. Mathematics 9(7), 2146 (2021). 10.3390/math9172146
  * [62] Zymbler, M., Kraeva, Y.: Discovery of time series motifs on Intel many-core systems. Lobachevskii J. Math. 40(12), 2124–2132 (2019). 10.1134/S199508021912014X
  * [63] Zymbler, M., Kraeva, Y.: Parallel algorithm for time series motif discovery on graphics processor. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering 9(3), 17–34 (2020). 10.14529/cmse200302, (in Russian)
  * [64] Zymbler, M., Kraeva, Y., Latypova, E., Kumar, S., Shnayder, D., Basalaev, A.: Cleaning sensor data in smart heating control system. In: Proceedings of 2020 Global Smart Industry Conference, GloSIC 2020, Chelyabinsk, Russia, November 17-19, 2020. pp. 375–381 (2020). 10.1109/GloSIC50886.2020.9267813
  * [65] Zymbler, M., Polyakov, A., Kipnis, M.: Time series discord discovery on Intel many-core systems. In: Sokolinsky, L., Zymbler, M. (eds.) 13th International Conference, PCT 2019, Kaliningrad, Russia, April 2-4, 2019, Revised Selected Papers. Communications in Computer and Information Science. vol. 1063, pp. 168–182. Springer, Cham (2019). 10.1007/978-3-030-28163-2_12

## Appendix

######  Lemma 5.1.

Let us have the time series $T$, $|T|=n$, and two its $m$-length subsequences,
$T_{i,\,m}$ and $T_{i,\,m+1}$, where $1\leq i\leq n-m$ and $3\leq m\ll n$.
Then the following holds:

$\mu_{T_{i,\,m+1}}=\dfrac{1}{m+1}\bigl{(}m\mu_{T_{i,\,m}}+t_{i+m}\bigr{)},$
$\sigma^{2}_{T_{i,\,m+1}}=\dfrac{m}{m+1}\Bigl{(}\sigma_{T_{i,\,m}}^{2}+\dfrac{1}{m+1}\bigl{(}\mu_{T_{i,\,m}}-t_{i+m}\bigr{)}^{2}\Bigr{)}.$

######  Proof 5.2.

First, let us prove the equation on the mean value. According to the
definition

$\mu_{T_{i,\,m}}=\dfrac{1}{m}\sum_{k=0}^{m}t_{i+k}.$

Then

$\sum_{k=0}^{m}t_{i+k}=m\mu_{T_{i,\,m}}.$

Next, let us consider $\mu_{T_{i,\,m+1}}$:

$\mu_{T_{i,\,m+1}}=\dfrac{1}{m+1}\sum_{k=0}^{m}t_{i+k}=\dfrac{1}{m+1}(m\mu_{T_{i,\,m}}+t_{i+m}),$

so, the first equation is proved. Further, let us prove the equation on the
standard deviation. According to the definition

$\sigma^{2}_{T_{i,\,m}}=\dfrac{1}{m}\sum_{k=0}^{m-1}t^{2}_{i+k}-\mu^{2}_{T_{i,\,m}}.$

Then

$\sum_{k=0}^{m-1}t^{2}_{i+k}=m(\sigma^{2}_{T_{i,\,m}}+\mu^{2}_{T_{i,\,m}}).$

Next, let us consider $\sigma^{2}_{T_{i,\,m+1}}$:

$\sigma^{2}_{T_{i,\,m+1}}=\dfrac{1}{m+1}\sum_{k=0}^{m}t^{2}_{i+k}-\mu^{2}_{T_{i,\,m+1}}=\dfrac{1}{m+1}\bigl{(}m(\sigma^{2}_{T_{i,\,m}}+\mu^{2}_{T_{i,\,m}})+t^{2}_{i+m}\bigr{)}-\mu^{2}_{T_{i,\,m}}.$

Employing the above-proved proposition on the mean value, we obtain

$\sigma^{2}_{T_{i,\,m+1}}=\dfrac{1}{m+1}\bigl{(}m(\sigma^{2}_{T_{i,\,m}}+\mu^{2}_{T_{i,\,m}})+t^{2}_{i+m}\bigr{)}-\bigl{(}\dfrac{1}{m+1}(m\mu_{T_{i,\,m}}+t_{i+m})\bigr{)}^{2}.$

Performing operations in parentheses and further collecting terms, we obtain

$\sigma^{2}_{T_{i,\,m+1}}=\dfrac{1}{m+1}\bigl{(}m(\sigma^{2}_{T_{i,\,m}}+\mu^{2}_{T_{i,\,m}})+t^{2}_{i+m}-\dfrac{1}{m+1}(m\mu_{T_{i,\,m}}+t_{i+m})^{2}\bigr{)}=$
$=\dfrac{1}{m+1}\Bigl{(}m\sigma^{2}_{T_{i,\,m}}+\dfrac{1}{m+1}\bigl{(}m(m+1)\mu^{2}_{T_{i,\,m}}-m^{2}\mu_{T_{i,\,m}}+(m+1)t^{2}_{i+m}-t^{2}_{i+m}-2m\mu_{T_{i,\,m}}t_{i+m}\bigr{)}\Bigr{)}=$
$=\dfrac{1}{m+1}\Bigl{(}m\sigma^{2}_{T_{i,\,m}}+\dfrac{m}{m+1}\bigl{(}\mu^{2}_{T_{i,\,m}}-2\mu_{T_{i,\,m}}t_{i+m}+t^{2}_{i+m}\bigr{)}\Bigr{)}=$
$=\dfrac{m}{m+1}\Bigl{(}\sigma^{2}_{T_{i,\,m}}+\dfrac{1}{m+1}\bigl{(}\mu_{T_{i,\,m}}-t_{i+m}\bigr{)}^{2}\Bigr{)}.$

This concludes our proof.