Admissible Time Series Motif Discovery with Missing Data
Yan Zhu, ‡Abdullah Mueen, Eamonn Keogh
University of California, Riverside, ‡University of New Mexico
yzhu015@ucr.edu, eamonn@cs.ucr.edu, ‡mueen@cs.unm.edu

Abstract
The discovery of time series motifs has emerged as one of the most
useful primitives in time series data mining. Researchers have
shown its utility for exploratory data mining, summarization,
visualization, segmentation, classification, clustering, and rule
discovery. Although there has been more than a decade of extensive
research, there is still no technique to allow the discovery of time
series motifs in the presence of missing data, despite the welldocumented ubiquity of missing data in scientific, industrial, and
medical datasets. In this work, we introduce a technique for motif
discovery in the presence of missing data. We formally prove that
our method is admissible, producing no false negatives. We also
show that our method can “piggy-back” off the fastest known motif
discovery method with a small constant factor time/space overhead.
We will demonstrate our approach on diverse datasets with varying
amounts of missing data.

1 INTRODUCTION
Time series motifs are short approximately repeated patterns within
a longer time series dataset. The fact that such patterns are
conserved often suggests underlying structure and regularities that
can be exploited in many ways. Some examples include: rule
discovery [16][13], forecasting [10], or building better classifiers
[2]. However, despite over a decade of active research, there is no
known method to allow the discovery of the motif in the presence
of missing data.
Paradoxically, in spite of improvements in sensor technology,
missing data is becoming more prevalent. This is because sensors
are now so cheap that we are willing to deploy them in hostile
environments with intermittent and unreliable bandwidth [17]. In
many cases, sensors have become “throwable” and disposable [7].
Figure 1 shows an example of a motif in a music processing
domain [12]. Note that both occurrences of the motif contain
sections of missing data.
At 34 seconds

At 75 seconds

0

100

200

Figure 1: A four-second long motif that appears in the pitch contour
time series of a Cypriot folk song, Kotsini Trantafillia (Red Rose-tree).
Note that both occurrences have multiple instances of missing data
[12].

Here the data could be missing for one of two reasons. It might be
meaningful, i.e., the vocalist may not have produced a sound at
these times. It is also possible that the missing data reflects poor
audio quality, loose wires, etc. We are agnostic to such issues in
this work. We simply note that in either case, there is no motif
discovery algorithm defined in the literature for such data.
More generally, in time series data mining, missing data is often
handled by filling in the missing values with some interpolation or

imputation method [8] and running the analysis unaltered.
However, as we shall demonstrate, no matter what interpolation or
imputation method is used, and no matter what current motif
discovery algorithm is used, this can result in false negatives.
To solve this problem, we introduce a novel algorithm that does not
allow false negatives. The method may allow occasional false
positives, but since the discovered motifs are typically examined by
the human eye [11][13][25], or some subsequent analysis, the false
positives (if any) can be filtered out at a later stage.
We call our algorithm Motif Discovery with Missing Data, MDMS.
Our MDMS algorithm is built on top of the recently introduced
Matrix Profile data structure [22]. Matrix Profile-based motif
discovery has been shown to have following advantageous features,
all of which are inherited by our proposed MDMS:
•
It is simple and parameter-free. In contrast, the more general
motif discovery algorithms require building and tuning spatial
access methods and/or hash functions.
•
It requires an inconsequential space overhead, just O(n) with
a small constant factor.
•
While the exact algorithm is scalable, for extremely large
datasets we can compute the results in a timely fashion,
allowing ultra-fast approximate solutions.
•
Having computed the motifs for a dataset, we can
incrementally update them efficiently. Thus for reasonable
arrival rates (< 20Hz), we can effectively maintain exact
motifs on streaming data forever.
•
The algorithm is easy to parallelize, both on multicore
processors and in distributed systems.
The rest of the paper is organized as follows. Section 2 reviews
related work and introduces the necessary background materials
and definitions. In Section 3, we introduce our algorithms. Section
4 gives a detailed empirical evaluation. Finally, in Section 5 we
offer conclusions and directions for future work.

2 RELATED WORK AND BACKGROUND
Time series motif discovery has been an active area of research for
over a decade [3][11][25]. However, to the best of our knowledge,
there is no existing algorithm that can find motifs in the presence
of missing data. In principle, any motif discovery algorithm could
be used with missing data, if we use some imputation algorithm to
fill in the missing values. There are hundreds of imputation
algorithms in the literature (see [21] and the references therein) to
choose from, but no matter which one we use, we may obtain false
negatives with respect to the oracle data (the true underlying data,
without missing values).
To see this let us consider a simple example. Suppose we have a
dataset that is composed of just three (sub)sequences:
A={0,2,0,2}, B={0,2,0,2}, C={0,-1,0,2}
Clearly the pair A|B is a perfect motif. Now suppose that the second
value in A is missing. The most obvious imputation technique is
interpolation from the two neighbors of the missing data point,
giving us Amiss = {0,0,0,2}. As we can see in Figure 2, this one
change means that we no longer discover the pair A|B as a perfect

motif, but instead we are lead to believe that the pair A|C is the best
motif in the dataset.
C

B

B

C

A

Amissing

Figure 2: left) A contrived dataset in which the pair A|B is a perfect
motif. right) If A had its second value replaced by the most common
imputation algorithm, we would fail to discover A|B as the motif.

While this is a trivial example, it is easy to see that no matter what
imputation algorithm is used, using a simple adversarial augment
we can always construct datasets for which the classic motif
algorithms would produce false negatives. This problem is even
more severe with complex datasets that contain a lot of highfrequency patterns and noise, for example seismology datasets.
In [26] the authors showed two repeating earthquake patterns that
appeared approximately 14 years apart. We have found that if there
were just a handful of missing data points in one of the earthquake
samples, we would be unable to detect a match between them with
any common imputation method, as such high-frequency and noisy
data defies the assumptions that most imputation techniques
assume.
Figure 2 showed that imputation methods can produce possible
false negatives even if we have a random single value missing in
the data. Moreover, another disadvantage of current imputation
methods is that they cannot predict block-missing data well [24]. In
some circumstances, due to malfunction of the sensor or other
anomaly factors, we may lose reading from a sensor at consecutive
timestamps (instead of sparse, single missing timestamps at random
locations of the data). This is often called block-missing data. In [24]
the authors proposed a state-of-the-art spatial-temporal imputation
method to predict the block-missing data by learning from not only
the real-valued reading of the same sensor, but also from the
reading of several geographically nearby sensors. The information
from other sensors greatly improved their imputation accuracy.
However, this method does not apply when we only have access to
one single sensor, or when all the sensors contain block-missing
data at the same timestamps due to regional power outage or
communication errors.
Finally, note that our problem is not artificial or contrived in
anyway. The literature is replete with examples of data analysts
frustrated by the inability to perform motif discovery in the
presence of missing data. For example, a recent paper studies
recurrent water consumption behavior by Australian consumers
[19]. The authors observe “A small proportion of all hourly
readings are missing…, probably due to server failure or
maintenance”. The authors observe “A small proportion of all
hourly readings are missing…, probably due to server failure or
maintenance.” The authors realize that any imputation method used
here has a risk of producing false negatives by noting that “…hourly
water consumption is highly unpredictable, we ignored the points
of missing hours for the routine discovery, rather than
approximating missing readings.” However, their solution of
ignoring some data runs the risk of missing interesting patterns.

2.1 Dismissing Apparent Solutions
In this section, we continue the discussion of related work, while
explicitly dismissing some proposed solutions to the task at hand.

The last decade has seen several distance measures for handling
uncertain time series, including PROUD [23], DUST [15], PBRQ,
MUNICH etc. One might believe that these measures could be used
to replace the Euclidean Distance subroutine in an existing motif
discovery algorithm. However, we believe this is not possible for
the following reasons:
•
These methods assume not missing data, but uncertain data.
For example, they can address the situation where an
observation is not known precisely but comes from some
known distribution. In [15] for example, they explicitly model
the normal, exponential, and uniform error distributions.
However, for generality, we do not wish to make such strong
assumptions.
•
Even if we assume that we could somehow avail of an existing
uncertain distance measure, none of them are metrics (only
measures). However, all speedup techniques for exact motif
discovery that we are aware of require and exploit metric
properties [3][11]. This suggests we must resort to a O(n2m)
brute force search (n is the length of the time series, m the
motif length). However, our proposed method is O(n2). As m
may be in the thousands (see Section 4), this suggests a three
order of magnitude time difference.
•
Finally, we want to be able to guarantee that our search
produces no false negatives in the face of missing data. To the
best of our knowledge, no existing uncertain time series
similarity measures can support this requirement.
The reader may not appreciate why our task-at-hand is hard,
because the analogue problem with strings is trivial. Suppose we
are asked to compare the following text strings “Norwegian blue”
and “Norwegian wood” under the Hamming Distance, and we
encounter missing values, represented here by “*”.
Norwe*ian_blu*
N*rwe*ian_wood
We can easily compute both the lower and upper bound of the
Hamming Distance. In the former case, we would assume that all
the missing values in one word are the same as their counterparts
in the other word. Given that the only letters we can be sure are
differing are “woo” vs “blu,” and we have a lower bound of three.
In the latter case, we would assume that all the missing values in
one word are different than their counterparts in the other word.
These three pessimistic differences combined with the three
observed differences give us an upper bound of six.
However, consider the time series version of this problem. Suppose
we have the following two time series:
[ 0.5, 0.1, ***, ***, -0.6, -0.7, 0.0 ]
[ 0.3, 0.1, 1.1, ***, -0.6, ***, 0.1 ]

One might consider applying similar logic here. For example,
accumulating 0.2 error (i.e. |0.5 - 0.3|) from the first pair of
numbers, then 0.0 error (i.e. |0.1 - 0.1|) from the second pair, etc.
However, the critical difference is that the time series must be
normalized before comparison. This is because, aside from the rare
and well understood exceptions [5], it is meaningless to compare
time series without normalizing them first. This presents a problem
as z-normalization (the most common normalization technique [4])
requires us to know the exact mean and standard deviation of the
data, which are undefined when we have a single missing data point.
Thus, for any pair of corresponding known points (for example the
0.5 and the 0.3 in the above), it is possible that the true (had we
known the mean and standard deviation of the data to allow the
correct normalization) difference between them could stay the same,

increase, or decrease by arbitrary amounts. This suggests that
producing tight upper or lower bounds will be nontrivial.

2.2 Definitions and Notations
In this section, we introduce all of the necessary notations and
definitions needed to explain our solution. Some of the text is
redundant with [26], but is included for completeness.
We begin by defining the data type of interest, time series:
Definition 1: A time series T is a sequence of real-valued numbers
ti: T = t1, t2, ..., tn, where n is the length of T.
For motif discovery, we are not interested in the global properties
of time series, but in the local regions, known as subsequences:
Definition 2: A subsequence Ti,m of a time series T is a continuous
subset of the values from T of length m starting from position i.
Formally, Ti,m = ti, ti+1,…, ti+m-1, where 1 ≤ i ≤ n-m+1.
We can generalize both definitions to allow for the possibly that at
least one value is missing. For clarity, we differentiate such time
series (including subsequences) with a “bar.” To keep both in a
common notation and for implementation purposes (i.e. Matlab,
etc.), we use NaNs as a placeholder for missing values.
Definition 3: A missing value time series 𝑇 is a sequence of values
that are either real-valued numbers or NaNs, 𝑡̅i: 𝑇 = 𝑡̅1, 𝑡̅2, ..., 𝑡̅n,
where n is the length of 𝑇.
In the rest of the paper, we assume that T is the actual time series
of 𝑇 before the missing values were created by some process. That
is to say, we would have obtained T instead of 𝑇 if the sensors were
functioning properly. We do not know T precisely, but we have 𝑡i̅ =
ti (1 ≤ i ≤ n) if 𝑡̅i ≠NaN.
The following restates definitions introduced with the STOMP
framework [26], but are reviewed for clarity.
We can take any subsequence from a time series and compute its
distance to all sequences. We call an ordered vector of such
distances a distance profile:
Definition 4: A distance profile Di of time series T is a vector of
the Euclidean distances between a given query subsequence Ti,m
and each subsequence in the time series T. Formally, Di = [di,1,
di,2,…, di,n-m+1], where di,j (1 ≤ i, j ≤ n-m+1) is the distance between
Ti,m and Tj,m.
We assume that the distance is measured by Euclidean distance
between z-normalized subsequences [4]. Unless otherwise stated,
in the rest of the paper, we always use di,j to represent the znormalized Euclidean distance between Ti,m and Tj,m.
We are interested in finding the nearest neighbors of all
subsequences in T, as the closest pairs are the classic definition of
time series motifs [3][11]. Note that by definition, the ith location
of distance profile Di is zero and very close to zero just before and
after this location. In the literature, such matches are called trivial
matches [11]. We avoid such matches by ignoring an “exclusion
zone” of length m/4 before and after the location of the query. In
practice, we set di,j to infinity (i-m/4 ≤ j ≤ i+m/4) while evaluating
Di.
We use a vector called a matrix profile to represent the distances
between all subsequences and their nearest neighbors:
Definition 5: A matrix profile P of time series T is a vector of the
Euclidean distances between each subsequence Ti,m and its nearest
neighbor (closest match) in time series T. Formally, P = [min(D1),
min(D2),…, min(Dn-m+1)], where Di (1 ≤ i ≤ n-m+1) is the distance
profile Di of time series T.
This vector is called the matrix profile, because one way to compute
it would be to compute the full distance matrix of all pairs of

subsequences in time series T, and then evaluate the minimum
value of each column. Figure 3 illustrates both a distance profile
and a matrix profile created on the same raw time series T.
Q, query of
length m

T, synthetic data

D, a distance profile

Note that |D| = |T|-|Q|+1
m/4 m/4

0

2,500

T, synthetic data
P, a matrix
profile

Note that |P| = |T|-|Q|+1

0

2,500

Figure 3: top) One distance profile (Definition 4) created from a
random subsequence Q of T. If we created distance profiles for all
possible subsequences of T, the element-wise minimum of this set
would be the matrix profile (Definition 5) shown at (bottom). Note that
the two lowest values in P are at the location of the 1st motif [3][11].

The ith element in the matrix profile P tells us the Euclidean
distance from the subsequence Ti,m to its nearest neighbor in the
time series T. However, it does not tell us where that neighbor is
located. This information is recorded in a companion data structure
called the matrix profile index.
Definition 6: A matrix profile index I of time series T is a vector of
integers: I=[I1, I2, … In-m+1], where Ii=j if di,j = min(Di).
By storing the neighboring information this way, we can efficiently
retrieve the nearest neighbor of query Ti,m by accessing the ith
element in the matrix profile index.
For clarity of presentation, we have confined this work to the single
dimensional case; however, nothing about our work intrinsically
precludes generalizations to multidimensional data.
To briefly summarize this section: we can create two meta time
series, the matrix profile and the matrix profile index, to annotate a
time series T with the distance and location of all its subsequences’
nearest neighbors within itself. For a time series T of length n, the
recently introduced STOMP algorithm [26] is able to compute the
two meta files with a mere O(n2) time complexity and O(n) space
complexity, which enables a fast motif discovery in a massive time
series. However, the STOMP algorithm is not applicable to any
time series 𝑇 with missing values (NaNs). Here we claim our
MDMS algorithm can solve the missing data problem with the
same time and space complexity. We leave the detailed discussion
of the algorithm to Section 3.

2.3 Pseudo Missing Data
Before introducing our solution to the missing data problem, we
take the time to point out that the problem is more general than one
may assume when we consider the generalization of the Pseudo
Missing Data (PMD). Informally, we define PMD as data that
technically is not missing, but effectively is. Figure 4 illustrates
three kinds of PMD frequently encountered.
(A) Spike

100
50
0
-50
-100

(B) Noise

(C) Single Byte
Maximum Value

EOGright
0

400

800

1200

Figure 4: A snippet of an Electrooculogram (EOG) exhibits three
kinds of pseudo missing data

In Figure 4.A, we see a “spike.” Given what we know about this
domain, it is inconceivable that the human eye could move fast

enough to produce this data, so it is clearly an artifact. Likewise, in
Figure 4.B, the dramatic increase in variance suggests that this
section of data is not likely to faithfully represent the underlying
physical process. Finally, in Figure 4.C the perfectly flat plateau is
not reflective of reality, but is simply a region where the physical
process exceeds the 8-bit precision available to record it.
In all three cases, the best thing to do would be to treat the data as
missing. Note that this decision is domain-dependent; there are
clearly domains where a short spike represents some physical event,
or where a perfectly flat plateau represents a physical limitation,
which is not a quirk of the hardware/software use.

3 ALGORITHMS
3.1 An Intuitive Preview
We begin by previewing our solution. As shown in [22][26], if T is
used to compute a matrix profile then finding the motifs is trivial.
The location of the smallest pair of values in the matrix profile is
also the location in T of the optimal motif pair. Moreover, other
definitions of motifs, such as the top-K motifs or range motifs
[11][22], can also be easily extracted from the matrix profile. Given
this, our solution to the missing data problem is to build a special
matrix profile using 𝑇. This special matrix profile will be very
similar to the true matrix profile, and be a (in general, very tight)
lower bound for it. If we use the existing motif extraction
algorithms [22][26] to pull out motifs from this matrix profile, we
may have false positives, but we will have no false negatives [5].
Thus most of our contribution outlined below is to show how we
can build this special matrix profile.

3.2 Lower Bound Matrix Profile
To create the special matrix profile data structure, our MDMS
algorithm evaluates the z-normalized Euclidean distance between
every pair of subsequences within a missing value time series 𝑇.
Depending on whether or not the subsequences have missing values,
we may encounter three different situations. Assume that the pair
of subsequences under consideration is 𝑇𝑖,𝑚 and 𝑇𝑗,𝑚 .
•

•
•

Case 1: Neither 𝑇𝑖,𝑚 nor 𝑇𝑗,𝑚 contain any missing value
(NaN). Normally, this applies to most subsequence pairs
within time series 𝑇 if 𝑇 contains more real-valued numbers
than NaNs. The traditional exact z-normalized Euclidean
distance between 𝑇𝑖,𝑚 or 𝑇𝑗,𝑚 can be evaluated in this case.
Case 2: 𝑇𝑖,𝑚 contains missing values (NaNs) while 𝑇𝑗,𝑚 does
not, or vice versa.
Case 3: Both 𝑇𝑖,𝑚 and 𝑇𝑗,𝑚 have missing values (NaNs).

In cases 2 and 3, the exact z-normalized Euclidean distance 𝑑𝑖,𝑗
between 𝑇𝑖,𝑚 or 𝑇𝑗,𝑚 cannot be evaluated. However, we can
𝐿𝐵
evaluate a lower bound of the distance between 𝑇𝑖,𝑚 and 𝑇𝑗,𝑚 , 𝑑𝑖,𝑗
,
𝐿𝐵
such that 𝑑𝑖,𝑗
≤ min⁡(𝑑𝑖,𝑗 ). That is to say, no matter what the
missing values in 𝑇𝑖,𝑚 or 𝑇𝑗,𝑚 are, we guarantee that 𝑑𝑖,𝑗 , the actual
𝐿𝐵
distance between 𝑇𝑖,𝑚 and 𝑇𝑗,𝑚 , is no less than 𝑑𝑖,𝑗
.
Now we are ready to redefine Definition 4 in the context of missing
values. We keep those 𝑑𝑖,𝑗 values corresponding to Case 1
𝐿𝐵
unchanged, and use 𝑑𝑖,𝑗
to replace the 𝑑𝑖,𝑗 values corresponding to
Cases 2 and 3. We will then obtain a lower bound distance profile:
Definition 7: A lower bound distance profile 𝐷𝑖𝐿𝐵 of a missing
value time series 𝑇 is a vector 𝐷𝑖𝐿𝐵 =[ ⁡𝑑̅𝑖,1 , 𝑑̅𝑖,2 ,…, 𝑑̅𝑖,𝑛−𝑚+1 ,],
where 𝑑̅𝑖,𝑗 =𝑑𝑖,𝑗 (1 ≤ i, j ≤ n-m+1) if neither 𝑇𝑖,𝑚 nor 𝑇𝑗,𝑚 contains
𝐿𝐵
NaNs (Case 1), and 𝑑̅𝑖,𝑗 =𝑑𝑖,𝑗
(1 ≤ i, j ≤ n-m+1) otherwise.

Similarly, we can redefine Definition 5 and Definition 6 in the
context of missing values.
Definition 8: A lower bound matrix profile 𝑃 𝐿𝐵 of a missing value
time series 𝑇̅ is a vector of the lower bound Euclidean distances
between each subsequence 𝑇̅𝑖,𝑚 and its nearest possible neighbor
(closest possible match) in 𝑇̅ . Formally, 𝑃 𝐿𝐵 = [min( 𝐷1𝐿𝐵 ),
𝐿𝐵
min(𝐷2𝐿𝐵 ),…, min(𝐷𝑛−𝑚+1
)], where 𝐷𝑖𝐿𝐵 (1 ≤ i ≤ n-m+1) is the
lower bound distance profile 𝐷𝑖𝐿𝐵 of the time series T.
Definition 9: A lower bound matrix profile index 𝐼𝐿𝐵 of 𝑇̅ is a
𝐿𝐵
vector of integers: 𝐼𝐿𝐵 =[𝐼1𝐿𝐵 , 𝐼2𝐿𝐵 , … 𝐼𝑛−𝑚+1
], where 𝐼𝑖𝐿𝐵 =j if 𝑑̅𝑖,𝑗 =
𝐿𝐵
min(𝐷𝑖 ).
The lower bound matrix profile gives us an optimistic
approximation of the distance between every subsequence and its
nearest possible neighbor in 𝑇̅. This approximation may produce
false positives, but it will not allow false negatives.
We believe this is the best approach for the task at hand: as [26]
shows, the cost to filter out false positive motifs is very low once
we have the matrix profile, but we cannot afford the occurrence of
any false negatives, since they may include the most important
patterns in the time series. In the following sections, we will
introduce the lower bounds corresponding to Cases 2 and 3
respectively; then we will introduce our MDMS algorithm, which
evaluates a lower bound of the matrix profile.

3.3 The Lower Bound Euclidean Distance
3.3.1 Case 2
Let us first consider Case 2, where 𝑇𝑖,𝑚 contains missing values
(NaNs) while 𝑇𝑗,𝑚 does not, or vice versa. Without loss of
generality, for now we assume 𝑇𝑖,𝑚 is the subsequence that
contains missing values, and 𝑇𝑗,𝑚 is the subsequence without
missing values. Figure 5 shows a visual example of this case.
𝑇𝑖,𝑚
𝑇𝑗,𝑚
0

40

𝑀

80

𝑅

Figure 5: top) A subsequence with missing values. bottom) A
subsequence without missing values.

Here 𝑀 = {𝑘⁡|⁡1 ≤ 𝑘 ≤ 𝑚⁡and⁡𝑡𝑖+𝑘−1 = 𝑁𝑎𝑁⁡} is a set of the
locations of missing values within 𝑇𝑖,𝑚 , and 𝑅 = {𝑘⁡|⁡1 ≤ 𝑘 ≤
𝑚⁡and⁡𝑡𝑖+𝑘−1 ≠ 𝑁𝑎𝑁⁡and⁡𝑡𝑗+𝑘−1 ≠ 𝑁𝑎𝑁⁡} is the intersection of
the real-valued locations within 𝑇𝑖,𝑚 and those within 𝑇𝑗,𝑚 . To
evaluate the lower bound distance of the two subsequences, first we
need to z-normalize them [4][11][22][26].
We assume that 𝜇𝑖 and⁡𝜎𝑖 are the mean and standard deviation of
𝑇𝑖,𝑚 , 𝜇𝑗 and⁡𝜎𝑗 are the mean and standard deviation of 𝑇𝑗,𝑚 . Note
that because 𝑇𝑖,𝑚 has missing values, we cannot evaluate 𝜇𝑖 ⁡and 𝜎𝑖 ;
however, we can treat them as variables. Assume 𝑑𝑖,𝑗 is the
distance between 𝑇𝑗,𝑚 and the oracle subsequence of 𝑇𝑖,𝑚 , we can
easily obtain a lower bound distance of 𝑑𝑖,𝑗 by ignoring the missing
part 𝑀 in Figure 5:
2

𝑡𝑖+𝑘−1 − 𝜇𝑖 𝑡𝑗+𝑘−1 − 𝜇𝑗
𝐿𝐵
𝑑𝑖,𝑗
= √min ∑ (
−
)
𝜇𝑖 ,𝜎𝑖
𝜎𝑖
𝜎𝑗
𝑘∈𝑅

2

𝐿𝐵
Assume 𝑓1 = 𝑑𝑖,𝑗
. We can linearly transform 𝑓1 as:

(1)

2

𝜎𝑅

2

𝑡𝑖+𝑘−1 −𝜇

𝐿𝐵
𝑓1 = ⁡𝑑𝑖,𝑗
= ( 𝑗 ) 𝑚𝑖𝑛 ∑𝑘∈𝑅(
𝜎𝑗

𝜎

𝜇,𝜎

−

𝑡𝑗+𝑘−1 −𝜇𝑗𝑅 2
)
𝜎𝑗𝑅

(2)

𝑅
𝑅 is shown in Figure 5. We assume 𝑇𝑖,𝑚
is the real-valued part of
𝑅
𝑇𝑖,𝑚 , 𝑇𝑗,𝑚
is the subset of 𝑇𝑗,𝑚 corresponding to 𝑅, 𝜇𝑖𝑅 and 𝜎𝑖𝑅 are
𝑅
the mean and standard deviation of 𝑇𝑖,𝑚
, 𝜇𝑗𝑅 and 𝜎𝑗𝑅 are the mean
𝑅
and standard deviation of 𝑇𝑗,𝑚 .
𝑅
Figure 6 visualizes (2). Note that 𝑇𝑗,𝑚
is z-normalized in (2), so its
𝐿𝐵
offset and scale are fixed in Figure 6; to obtain 𝑑𝑖,𝑗
, we would like
𝑅
to adjust 𝜇 (corresponding to the offset of 𝑇𝑖,𝑚
) and 𝜎
𝑅
(corresponding to the scale of 𝑇𝑖,𝑚 ), such that the Euclidean
𝑅
𝑅
distance between 𝑇𝑖,𝑚
and 𝑇𝑗,𝑚
is minimized.

scale

𝑅
𝑇𝑖,𝑚

𝐿𝐵
We can see from (6) that ⁡𝑑𝑖,𝑗
is dependent on 𝜎𝑗 (controlled by the
𝐿𝐵
missing part of 𝑇𝑗,𝑚 in Figure 7): the larger 𝜎𝑗 , the smaller ⁡𝑑𝑖,𝑗
.
𝐿𝐵
Note that 𝜎𝑗 can be as large as +∞; in that case ⁡𝑑𝑖,𝑗 becomes zero.
This is a very undesirable lower bound, as any pair of missing value
subsequences can be reported as a motif, even if they look very
different from each other. Figure 8 shows an example of this.

offset

𝑅
𝑇𝑗,𝑚

0

Figure 7). A diligent reader may have noticed that the lower bound
𝐿𝐵
expression 𝑑𝑖,𝑗
in (5) subsumes (1), the lower bound expression in
Case 2. We can visualize this in Figure 7: if we remove 𝑀𝑗 , the
problem is transformed to Case 2. Therefore, we can directly
substitute (3), the result of Case 2, into (5):
1
𝜎𝑗𝑅 √|𝑅| min
𝑖𝑓⁡𝑞𝑖,𝑗 ≤ 0
𝜎𝑗 𝜎𝑗
𝐿𝐵
(6)
⁡𝑑𝑖,𝑗
=
1
2
𝜎𝑗𝑅 √|𝑅|(1 − 𝑞𝑖,𝑗
) min
𝑖𝑓⁡𝑞𝑖,𝑗 > 0
𝜎𝑗 𝜎𝑗
{

40

80

𝑅

Figure 6: Different setting of 𝝁⁡and 𝝈 changes the offset and the scale
of 𝑻𝑹𝒊,𝒎 . Note that the offset and scale of 𝑻𝑹𝒋,𝒎 are fixed.

By solving

𝜕𝑓1
𝜕𝜇

= 0 and

𝜕𝑓1
𝜕𝜎

= 0, and substituting 𝜎 and 𝜇 back into

(2), we have:
𝜎𝑗𝑅
𝜎𝑗

𝐿𝐵
⁡𝑑𝑖,𝑗
=

𝜎𝑗𝑅
{ 𝜎𝑗

√|𝑅|

𝑖𝑓⁡𝑞𝑖,𝑗 ≤ 0
(3)

2
) 𝑖𝑓⁡𝑞𝑖,𝑗 > 0
√|𝑅|(1 − 𝑞𝑖,𝑗

Here⁡𝑅 is the intersection of the real-valued locations within 𝑇𝑖,𝑚
and those within 𝑇𝑗,𝑚 , 𝑞𝑖,𝑗 is the Pearson Correlation Coefficient
𝑅
𝑅
between 𝑇𝑖,𝑚
and 𝑇𝑗,𝑚
:
𝑞𝑖,𝑗 =

∑𝑘∈𝑅 𝑡𝑗+𝑘−1 𝑡𝑖+𝑘−1 − |𝑅|𝜇𝑖𝑅 𝜇𝑗𝑅

(4)

|𝑅|𝜎𝑖𝑅 𝜎𝑗𝑅

The analysis of Case 2 is now complete. Let us turn to Case 3,
where both 𝑇𝑖,𝑚 and 𝑇𝑗,𝑚 contain missing values.

3.3.2

Case 3

Figure 7 shows an example of Case 3.
𝑇𝑖,𝑚

0

40

Figure 8: Two subsequences with missing values. The real-valued parts
of the subsequences look very different from each other, but if we fill
the missing parts with infinitely large numbers, the z-normalized
Euclidean distance of the two subsequences will become zero.

Fortunately, sensor readings normally have physical limits. The
accelerometer values on an iPhone 7 are limited to ± 8g (± 78.48
m/s2); virtually all medical sensors come with carefully specified
limits to meet regulations (i.e., EU directive 93/42/EEC mandates
that a pediatric lung ventilator monitor produces values in the range
of 0 to 125cmH20), etc. Therefore, we can assume that the missing
values in 𝑇𝑗,𝑚 are bounded by [𝑉𝑚𝑖𝑛 , 𝑉𝑚𝑎𝑥 ]. With this bound, we
can derive the following inequality for 𝜎𝑗 2 (we refer the interested
reader to [18] for the complete derivation):
𝑅𝑗

𝜎𝑗 2 ≤

2
𝐶 2 ∑𝑘∈𝑅𝑗 𝑡𝑗+𝑘−1 + |𝑅𝑗 | (𝐵 − 𝜇𝑗 𝐴)
+
4
𝑚

40

𝑅

Figure 7), 𝜇𝑗 𝑗 is the mean of the real-valued part of 𝑇𝑗,𝑚 ,𝐶 =
𝑉𝑚𝑎𝑥 − 𝑉𝑚𝑖𝑛 ,⁡𝐵 = 𝑉𝑚𝑎𝑥 𝑉𝑚𝑖𝑛 , 𝐴 = 𝑉𝑚𝑎𝑥 + 𝑉𝑚𝑖𝑛 . In practice, we
set 𝑉𝑚𝑖𝑛 and 𝑉𝑚𝑎𝑥 as the minimum and maximum value of the realvalued part of 𝑇𝑗,𝑚 .
𝐿𝐵
We can now evaluate⁡𝑑𝑖,𝑗
by substituting (7) back into (6):

𝑀𝑗

𝐿𝐵
𝑑𝑖,𝑗
={

80

𝑅
𝑅𝑗

2

Here 𝑓𝐿𝐵 (𝑗) = 𝜎𝑅𝑗 ⁄[

Figure 7: Two subsequences with missing values.

As both 𝑇𝑖,𝑚 and 𝑇𝑗,𝑚 have missing values, their means (𝜇𝑖 , 𝜇𝑗 )
and standard deviations (𝜎𝑖 ,⁡𝜎𝑗 ) can be arbitrary values. We have:
𝑡𝑖+𝑘−1 −𝜇𝑖

𝑑𝑖,𝑗 ≥ 𝑚𝑖𝑛 √𝑚𝑖𝑛 ∑𝑘∈𝑅 (
𝜇𝑗 ,𝜎𝑗

𝜇𝑖 ,𝜎𝑖

𝜎𝑖

−

𝑡𝑗+𝑘−1 −𝜇𝑗
𝜎𝑗

2

𝐿𝐵
) = 𝑑𝑖,𝑗

(5)

Here 𝑅 = {𝑘⁡|⁡1 ≤ 𝑘 ≤ 𝑚, 𝑡𝑖+𝑘−1 ≠ 𝑁𝑎𝑁, 𝑡𝑗+𝑘−1 ≠ 𝑁𝑎𝑁} is the
intersection of the real-valued locations within 𝑇𝑖,𝑚 and 𝑇𝑗,𝑚 (see

(7)

Here 𝑅𝑗 is a set of the locations of all the real values within 𝑇𝑗,𝑚 (see

𝑇𝑗,𝑚
0

80

√|𝑅|𝑓𝐿𝐵 (𝑗)

𝑖𝑓⁡𝑞𝑖,𝑗 ≤ 0

2
)𝑓𝐿𝐵 (𝑗) 𝑖𝑓⁡𝑞𝑖,𝑗 > 0
√|𝑅|(1 − 𝑞𝑖,𝑗

𝐶2
4

(8)

𝑅𝑗

+

∑𝑘∈𝑅𝑗 𝑡𝑗+𝑘−1 2 +|𝑅𝑗 |(𝐵−𝜇𝑗 𝐴)

]. Note that 𝑓𝐿𝐵 (𝑗)

𝑚

is based on 𝑅𝑗 , the real-valued part of 𝑇𝑗,𝑚 (recall Figure 7).
However, as in Case 3 both 𝑇𝑖,𝑚 and 𝑇𝑗,𝑚 have missing values, we
can analogously derive a similar lower bound expression as (8)
𝐿𝐵
based on 𝑇𝑖,𝑚 , and set the larger expression as 𝑑𝑖,𝑗
.
2

Formally, let 𝑓𝐿𝐵 (𝑝) = 𝜎𝑅𝑝 ⁄[

𝐶2
4

𝑅

+

∑𝑘∈𝑅𝑝 𝑡𝑝+𝑘−1 2 +|𝑅𝑝|(𝐵−𝜇𝑝 𝑝 𝐴)

the lower bound distance for Case 3 is:

𝑚

] , then

√|𝑅|⁡𝑚𝑎𝑥⁡(𝑓𝐿𝐵 (𝑖), 𝑓𝐿𝐵 (𝑗))
𝐿𝐵
⁡𝑑𝑖,𝑗
⁡= {
2
)𝑚𝑎𝑥⁡(𝑓𝐿𝐵 (𝑖), 𝑓𝐿𝐵 (𝑗))
√|𝑅|⁡(1 − 𝑞𝑖,𝑗

•

𝑖𝑓⁡𝑞𝑖,𝑗 ≤ 0
𝑖𝑓⁡𝑞𝑖,𝑗 > 0

(9)

The analysis of Case 3 is now complete. Finally, for completeness,
let us briefly discuss Case 1.

3.3.3

Case 1

As neither 𝑇𝑖,𝑚 nor 𝑇𝑗,𝑚 contain any missing value in Case 1, we
𝐿𝐵
set 𝑑𝑖,𝑗
as the exact z-normalized Euclidean distance between 𝑇𝑖,𝑚
and 𝑇𝑗,𝑚 , using the following equation ([26]):
𝐿𝐵
𝑑𝑖,𝑗
= √2𝑚(1 − 𝑞𝑖,𝑗 )

(10)

Note that all three cases use the same expression of 𝑞𝑖,𝑗 in (4).
Now that we have the lower bound distance for any subsequence
pair in 𝑇, we can also evaluate the lower bound matrix profile.

3.4 The MDMS Algorithm
The STOMP algorithm [26] can obtain the matrix profile of a time
series that is free of missing values, in O(n2) time with only O(n)
space; as we will now show, in the face of missing data, our MDMS
algorithm can obtain the lower-bound matrix profile with the same
time and space complexity. TABLE I summarizes the algorithm.
TABLE I. MDMS ALGORITHM
Algorithm MDMS(T, m)
Input: A missing value time series T, subsequence length m
Output: Lower bound matrix profile P and the associated lower bound
matrix profile index I of T
1
n ← Length(T), len ← n-m+1
2
vmax ← SlidingMax(T), vmin ← SlidingMin(T)
Z ← PadZero(T), B ← OneZero(T), X← ElementWiseSquare(Z)
3
μz, σz ← ComputeMeanStd(Z, m)
// see [14]
4
μb, σb ← ComputeMeanStd(B, m)
// see [14]
5
QZ ← SlidingDotProduct(Z[1:m], Z), QZ_first ← QZ //see [26]
6
QB ← SlidingDotProduct(B[1:m], B), QB_first ← QB //see [26]
7
BZ ← SlidingDotProduct(B[1:m], Z), BZ_first ← BZ //see [26]
8
ZB ← SlidingDotProduct(Z[1:m], B), ZB_first ← ZB//see [26]
9
10 BX ← SlidingDotProduct(B[1:m], X), BX_first ← BX //see [26]
11 XB ← SlidingDotProduct(X[1:m], B), XB_first ← XB //see [26]
12 P ← CalculateLBDistance (n, m, vmax, vmin, QZ, QB, BZ, ZB,
BX, XB, μz, σz, μb, σb, i), I← ones
// initialization
// in-order evaluation
13 for i = 2 to len
for j= len downto 2
// update dot product, see [26]
14
QZ[j] ← QZ[j-1]-Z[i-1]×Z[j-1]+Z[i+m-1]×Z[j+m-1]
15
QB[j] ← QB[j-1]-B[i-1]×B[j-1]+B[i+m-1]×B[j+m-1]
16
BZ[j] ← BZ[j-1]-B[i-1]×Z[j-1]+B[i+m-1]×Z[j+m-1]
17
ZB[j] ← ZB[j-1]-Z[i-1]×B[j-1]+Z[i+m-1]×B[j+m-1]
18
BX[j] ← BX[j-1]-B[i-1]×X[j-1]+B[i+m-1]×X[j+m-1]
19
XB[j] ← XB[j-1]-X[i-1]×B[j-1]+X[i+m-1]×B[j+m-1]
20
end for
21
QZ[1]←QZ_first[i], QB[1]←QB_first[i], BZ[1]←ZB_first[i]
22
ZB[1]←BZ_first[i], BX[1]←XB_first[i], XB[1]←BX_first[i]
23
D ← CalculateLBDistance (n, m, vmax, vmin, QZ, QB, BZ, ZB,
24
BX, XB, μz, σz, μb, σb, i)
P, I ← ElementWiseMin(P, I, D, i)
25
26 end for
27 return P, I

Before discussing the algorithm in detail, we first need to introduce
three important auxiliary time series (shown in line 3), 𝑍, 𝑋 and 𝐵.
•
We define 𝑍 = 𝑧1 , 𝑧2 , … 𝑧𝑛 , such that 𝑧𝑖 = 𝑡𝑖 if 𝑡𝑖 ≠ 𝑁𝑎𝑁 ,
and 𝑧𝑖 = 0 if 𝑡𝑖 = 𝑁𝑎𝑁. We can simply obtain 𝑍 by filling
zeros in the locations of 𝑇 where the data is missing.
•
We define 𝑋 = 𝑥1 , 𝑥2 , … 𝑥𝑛 = 𝑧1 2 , 𝑧2 2 , … 𝑧𝑛 2 .

We define 𝐵 = 𝑏1 , 𝑏2 , … 𝑏𝑛 , such that 𝑏𝑖 = 1 if 𝑡𝑖 ≠ 𝑁𝑎𝑁 ,
and 𝑏𝑖 = 0 if 𝑡𝑖 = 𝑁𝑎𝑁 . We can see that 𝐵 indicates the
locations of the real-valued numbers and missing values in 𝑇.
With these three auxiliary time series and the techniques introduced
in [14][26], we can evaluate any lower bound distances introduced
in the last section in O(1) time with O(n) space.
The MDMS algorithm is very similar to the STOMP framework
[26]. In line 2 we evaluate the maximum and minimum values of
the real-value part of every subsequence in T. Lines 4-5 evaluate
the mean and standard deviation of every subsequence in 𝑍 and 𝐵,
lines 6-11 initialize all the dot product vectors. We initialize the
lower bound matrix profile P and matrix profile index I in line 12.
Lines 13-26 iteratively evaluate the lower bound distance profile D,
and update P and I if necessary. The CalculateLBDistance
algorithm in lines 12 and 24 is shown in TABLE II.
TABLE II. ALGORITHM TO CALCULATE LOWERBOUND DISTANCE PROFILE
Algorithm CalculateLBDistance (n, m, vmax, vmin, QZ, QB, BZ, ZB,
BX, XB, μz, σz, μb, σb, i)
Input: the length n of time series T, the subsequence length m, the
maximum/minimum possible value vector vmax/vmin, dot product
vectors QZ, QB, BZ, ZB, BX, XB, means and standard deviations μz, σz,
μb, σb of time series Z and B, current subsequence index i
Output: Lower bound distance profile D
1
for j= 1 to n-m+1
2
ui ← ZB[j] / QB[j], uj ← BZ[j] / QB[j]
//𝜇𝑖𝑅 and 𝜇𝑗𝑅
3
vi← XB[j] / QB[j] - ui ^2, vj← BX[j] / QB[j] - uj ^2 //𝜎𝑖𝑅 and 𝜎𝑗𝑅
q ← (QZ[j] / QB[j]- ui × uj) / sqrt(vi × vj)
// (4)
4
if QB[j]==m then
// Case 1, |𝑅| = 𝑚
5
6
D[j] ←2×m×(1-q)
// (10)
else
7
if max(μb[i],μb[j])==1 then
//Case 2
8
if μb[i]>μb[j] then vo←vi, v←σz[i]^2
9
else vo←vj, v←σz[j]^2
10
end if
11
if q<=0 then D[j] ← QB[j] × vo / v // (3)
12
else D[j] ← QB[j] × vo / v × (1 - q^2)
13
end if
14
else
//Case 3
15
v1←vmax[i], v2←vmin[i], C←v1–v2, B←v1×v2, A←v1+v2
16
ur←μz[i]/μb[i], vr←(σz[i]^2+μz[i]^2)/μb[i]-ur^2
17
f1=vi/(ub[i]×(vr+B+ur×(ur–A))+C^2/4)
18
v1←vmax[j], v2←vmin[j], C←v1–v2, B←v1×v2, A←v1+v2
19
ur←μz[j]/μb[j], vr←(σz[j]^2+μz[j]^2)/μb[j]-ur^2
20
f2=vj/(ub[j]×(vr+B+ur×(ur–A))+C^2/4)
21
if q<=0 then D[j] ← QB[j] × max(f1, f2)
// (9)
22
else D[j] ← QB[j] × max(f1, f2) × (1 - q^2)
23
end if
24
end if
25
end if
26
27 end for
28 return D

The CalculateLBDistance algorithm evaluates all n-m+1 lower
bound distance values in D with equations (3), (9) and (10). Line 4
evaluates (4). Case 1 is handled in lines 5-6. Lines 8-14 handle Case
2, and Case 3 is evaluated in lines 15-24.
We can see that each loop of the CalculateLBDistance algorithm in
TABLE II (lines 2-26) can be evaluated in O(1) time, so the time
complexity of CalculateLBDistance is O(n). The space needed to
store all the vectors in TABLE I is O(n), and each loop in lines 1425 of TABLE I takes O(n) time. Therefore, the time complexity of
MDMS algorithm is O(n2) and the space complexity is O(n), the
same as STOMP [26].

Furthermore, we can see that unlike most imputation algorithms
[21], our MDMS algorithm is extremely model-free and parameterfree. The only inputs to the algorithm are the time series and the
subsequence length. In the next section, we will use two case
studies to demonstrate the robustness and efficacy of our ultra-fast,
parameter-free motif discovery algorithm in the face of missing
data.

4 EXPERIMENTAL EVALUATION
We begin by noting that all the code and data used in this work are
archived in perpetuity at [18]. In the following two case studies, we
consider both random-missing data (data with sparsely located,
random missing timestamps) and block-missing data (data with
consecutive missing timestamps).
For each case study, we compare our method with the commonly
used strawman in the literature to handle missing data: linear
imputation, as shown in Figure 9.

besides the two deep valleys). This shows the robustness of our
algorithm in the face of random missing data.
Next, we consider the “larger gaps” (or the block-missing data) case.
Here, instead of removing individual data points, we removed two
consecutive blocks of length 25. As shown in Figure 10 (orange
curve), the shape of the lower bound matrix profile still looks very
similar to that of the oracle matrix profile (red curve). We see only
two deep valleys at the vicinity of the motifs, so no false positive
patterns are discovered.
The result demonstrates that our lower bound matrix profile is
robust against producing false positives.
To test the robustness of our algorithm against false negatives, we
removed a block of missing data at the center of the second
earthquake pattern. The length of the missing block is 400, which
is 20% of the subsequence length. In Figure 11, we compared our
lower-bound matrix profile result with the matrix profile generated
by linear imputation.
We removed a block of
400 data points here

imputation

Raw data

Figure 9: To evaluate our method, we compare our result with that of
linear imputation.

Oracle Matrix Profile
Matrix Profile (MDMS)

4.1 Case Study: Seismological Data
Repeated pattern (i.e. motif) discovery is a fundamental tool in
seismology, which allows for the discovery of foreshocks, triggered
earthquakes, swarms, volcanic activity, and induced seismicity [25].
However, this domain is replete with missing data.
For example, a classic paper notes “A frequent dilemma in spectral
analysis (in seismology) is the incompleteness of the data record,
either in the form of occasional missing data or as larger gaps”
(our emphasis) [1]. In this experiment, we demonstrate that we can
handle both cases.
We consider a dataset for which we know the answer from external
sources. On April 30th, 1996, there was an earthquake of magnitude
2.12 in Sonoma County, California. Then, on December 29 th, 2009,
about 13.6 years later, there was another earthquake with a similar
magnitude. To allow the results to be visualized in a single plot, we
edited this data such that the two earthquakes happen just 15
seconds apart. We set the subsequence length as 2,000, which
corresponds to one second of data. As shown in Figure 10 (red
curve), when there is no missing data the matrix profile correctly
discovers the locations of the two earthquakes.
Raw data

Oracle Matrix Profile
Matrix Profile
(Missing points)
50 data points missing
Matrix Profile
(Missing blocks)
50 data points missing
0

25s

Figure 10: A raw seismograph contrived such that two earthquakes
from the same region happen 15 seconds apart. The matrix profile
computed with no missing data (red curve) finds the true event, as does
MDMS even in the presence of missing points (green curve) or missing
blocks (orange curve).

To test our algorithm for the “occasional missing data” case, we
randomly deleted 50 data points. As Figure 10 (green curve) shows,
the matrix profile is still minimized at the correct location, and there
are no false positives (no other small values in the matrix profile

Matrix Profile (imputation)
Oracle Matrix Profile
0

25s

Figure 11: We removed 400 consecutive data points at the center of the
second earthquake pattern. The oracle matrix profile computed with
no missing data (black curve) finds the true event, as does MDMS (red
curve) even in the presence of a large missing block. The Matrix Profile
generated after linear imputation (green curve) fails to capture the
minimum points within the oracle matrix profile.

We can see that the lower-bound matrix profile generated by our
MDMS algorithm (red curve) agrees closely with the oracle
matrix profile (black curve) at the vicinity of the two earthquake
patterns, while the matrix profile generated after linear imputation
(green curve) shows a high value at these locations.

0

2000

0

2000

Figure 12: The first motif found by the MDMS algorithm (right) in the
presences of a large missing block is identical with the first motif found
in the oracle data (left).

As a result, the MDMS algorithm successfully captures the 1 st
motif (as shown in Figure 12) even in the presence of a large
missing block within one earthquake pattern, while the imputation
method misses the 1st motif within the oracle data. This illustrates
two major strength of our algorithm over imputation methods.
Firstly, our algorithm does not allow false negatives. Secondly,
our algorithm is more robust to large missing blocks as it does not
change the data, while imputation method can change the data a
lot. In the next case study, we will further demonstrate the
robustness of MDMS in the presence of missing blocks.

4.2 Case Study: Activity Data from Video
Time series extracted from video often has missing data reflecting
“frame drops” due to bandwidth congestion [6]. To test our
algorithm in this context, we examine the activity dataset of [20].
This dataset consists of a 13.3 minute 10-fps video sequence of an
actor performing various activities. From this data, the original

authors extracted 721 channels of the optical flow time series, and
the length of each time series is 8,000. We consider the time series
corresponding to the 533th channel, which is suggestive of the
structure in places but is noisy. The data is shown in Figure 13.top,
the subsequence length is 120. From the oracle matrix profile in
Figure 13 (black curve), we can extract the 1st motif in the oracle
data. To test the performance of the MDMS algorithm, here we
remove 12 consecutive data points in the center of one of the 1st
motif patterns. In Figure 13, we compare our lower-bound matrix
profile result (red curve) with the matrix profile generated by linear
imputation (green curve).

We removed two blocks of length p
Data
Oracle
Matrix Profile

p=100
1st motif

p=400
1st motif

p=700

1st motif

p=800
1st motif

We removed a block of 12 data points here

Data
Oracle Matrix Profile
1st motif

Matrix Profile (MDMS)
1st motif

Matrix Profile (Imputation)
The deep valley disappears
0

8000

Figure 13: A raw activity time series. We removed 12 consecutive data
points in one of the 1st motif patterns in the time series. The oracle
matrix profile computed with no missing data (black curve) finds the
true motif starting at the 540th and the 622nd data points. With the
presence of 12 missing data points, the MDMS algorithm finds the
same motif as the oracle result (red curve), starting at 520th and 602nd
data points. The Matrix Profile generated after linear imputation
(green curve) fails to capture the two deep valleys within the oracle
matrix profile and thus misses the 1st motif.

We can see that the oracle matrix profile (black curve) shows two
apparent valleys at the locations of the 1st motif, as does the lower
bound matrix profile generated by MDMS (red curve). The 1st motif
discovered by the MDMS algorithm (shown in Figure 14.right) is
identical to the oracle motif (shown in Figure 14.left). In contrast,
the matrix profile generated by imputation does not have these
valleys, and thus misses the 1st motif of the oracle data.

Figure 15: Lower bound matrix profiles corresponding to different
missing block lengths. We removed 2 blocks of length p from the
seismograph. The oracle matrix profile (black curve) finds the true
motif. For p=100, p=400 and p=700, MDMS is able to find the true
event as the 1st motif. When p=800, MDMS finds a false positive as the
1st motif.

Note that the removed blocks are not within the two repeated
earthquake patterns. As a result, the lower bound matrix profiles
are the same as the oracle matrix profile at the vicinity of the two
repeated patterns, while lower at other locations. In other words, it
is easier to detect false positives with such setting.
We can see that when p=100 (5% of the subsequence length m), the
side valleys in the oracle Matrix Profile become deeper, and two
more side valleys show up at the vicinity of the removed blocks. As
p increases, all the side valleys become deeper and deeper. For
p=100, p=400 and p=700, we are able to find the true event as the
1st motif with MDMS. However, when p=800 (40% of m), the 1st
motif (corresponding to the minimum point of the lower bound
matrix profile) is no longer the true event. We show this false
positive motif pair in Figure 16.

0

2000

Figure 16: The 1st motif found by the MDMS algorithm when p=800.

0

120

0

120

Figure 14: The first motif found by the MDMS algorithm (right) is
identical to the first motif within the oracle data (left), despite a small
phase shift.

We can see from Figure 14 that though a large portion of the blue
pattern is missing, our MDMS algorithm still finds it very similar
with the red pattern. This example further demonstrates that our
algorithm is robust against missing true motif patterns.

4.3 Quantifying the Robustness of MDMS
As MDMS evaluates the lower bound matrix profile, it naturally
does not allow false negatives, but it can produce false positives.
Here we perform two “stress tests” to evaluate the robustness and
limitations of our MDMS algorithm against false positives.
We use the seismograph data in Figure 10 again for the stress test.
The subsequence length in this dataset is m=2,000.
We first test the sensitivity of MDMS over the length of missing
blocks. Here we remove two missing blocks of length p, located at
7.5s and 15s respectively, from the data. In Figure 15, we show
how the lower bound matrix profile varies as we increase p.

The two subsequences are both at the vicinity of the second missing
block in Figure 15. With a close inspection we can see why this
pair is reported as the 1st motif by MDMS. If we fill in the missing
part of each subsequence with their counterpart in the other
subsequence (shown by the arrows in Figure 16), the two can be
very similar to each other. Since MDMS does not allow false
negatives, it will capture and report this possible matching pair.
Furthermore, Figure 16 implies that p cannot be larger than 50%
of the subsequence length, otherwise MDMS will be able to find a
perfect match (with one subsequence missing the first 50% and the
other missing the second 50%). When p > 700 (35% of the
subsequence length), we are already very close to this limit.
Therefore, in Figure 15 we can see very deep side valleys at the
vicinity of the missing blocks, and we are prone to detect false
positives. When p ≤ 400 (20% of the subsequence length), the two
main valleys corresponding to the true events dominate, so we will
not detect false positives.
We have demonstrated that MDMS is robust against discovering
false positives when there are two missing blocks, and when the
length of the two blocks are within a reasonable range. Next, we
“stress test” the sensitivity of MDMS over the total number of
missing values in the data. We again use the seismograph dataset

in Figure 10, which consist of 50,000 data points. The lower bound
matrix profile results are shown in Figure 17.
10%

10%

20%

20%

30%

30%

40%

40%

Figure 17: Lower bound matrix profiles corresponding to various
percentage of data missing. left) random-missing data right) blockmissing data.

In the first run, we randomly selected 5,000, 10,000, 15,000 and
20,000 points to remove from the data. From Figure 17.left, we can
see that the scale of the matrix profile decreased as more points are
missing. However, even when 40% of the data is missing, the two
valleys corresponding to the true events (recall Figure 10) still
dominate. In the second run, we removed blocks of length 200 from
the data. The missing blocks were uniformly distributed, and the
number of missing blocks increased from 25 to 100. Figure
17.right shows that even when 30% of the data is missing, the two
main valleys still dominate. When 40% of the data is missing, the
two main valleys are no longer apparent, but we can still find the
true events as the 1st motif. The experiment demonstrates that
MDMS is robust against detecting false positives even if a large
percentage of data is missing.

5 CONCLUSIONS
We introduced what we believe to be the first time series motif
discovery algorithm that can find motifs in the presence of missing
data. The algorithm has the same time and space complexity as the
fastest known algorithm for motif discovery [26]. We formally
proved the admissibility of our algorithm, it does not produce any
false negatives. Experimental results show that our algorithm is
also robust against false positives even when a large portion of the
data is missing. Because our algorithm is based on creating a
special version of the matrix profile [22][26], our work may have
implications for other algorithms that can exploit the matrix profile,
including discord discovery and time series joins. The lower
bounds introduced can also be used to accelerate various length
motif discovery. We leave such considerations to future work.

6 REFERENCES
[1] Baisch, S. and Bokelmann, G.H. Spectral analysis with
incomplete time series: an example from seismology.
Computers and Geosciences, 25.7(1999): 739-750.
[2] Buza, K., and Schmidt-Thieme, L. Motif-based classification of
time series with bayesian networks and svms. In Advances in
Data Analysis, Data Handling and Business Intelligence 2009:
105-114.
[3] Chiu, B., Keogh, E., and Lonardi, S. Probabilistic discovery of
time series motifs. ACM SIGKDD 2003: 493-498.
[4] Ding, H., Trajcevski, G., Scheuermann, P., Wang, X., and
Keogh, E. J. Querying and Mining of Time Series Data:
Experimental Comparison of Representations and Distance
Measures. PVLDB 1.2 (2008): 1542-1552.
[5] Faloutsos, C., Ranganathan, M., Manolopoulos, Y. Fast
subsequence matching in time-series databases. SIGMOD
1994: 419-429.
[6] Gangadharan, D. et al. Video Quality Driven Buffer Sizing via
Frame Drops. RTCSA. 1 (2011): 319-328.

[7] Gao, H., Bi, S., Zhang, R. and Tang, S., The design of a
throwable two-wheeled reconnaissance robot. In 2012 IEEE
Robotics and Biomimetics (ROBIO), pp 2150-2155.
[8] Honaker, J. and King, G. What to do about missing values in
time‐series cross‐section data. American Journal of Political
Science, 54.2 (2010): 561-581.
[9] Kelly, G. Body temperature variability: a review of the history
of body temperature and its variability due to site selection,
biological rhythms, fitness, and aging. Med Rev. 11.4 (2006):
278.
[10] McGovern, A., Rosendahl, D.H., Brown, R.A. and
Droegemeier, K.K. Identifying predictive multi-dimensional
time series motifs: an application to severe weather
prediction. Data Mining and Knowledge Discovery, 22.1
(2011): 232-258.
[11] Mueen, A., Keogh, E., Zhu, Q., Cash, S. and Westover, B.
Exact Discovery of Time Series Motifs. SDM 2009: 473-484.
[12] Neocleous, A., Petkov, N., Schizas, C.: Finding repeating
Stanzas in monophonic folk songs of Cyprus. 6th Cyprus
Workshop on Signal Processing and Informatics 2013. pp 26.
[13] Patnaik, D., Manish, M., Sharma, R. K., and Ramakrishnan, N.
Sustainable Operation and Management of Data Center
Chillers Using Temporal Data Mining. KDD 2009: 1305-1314.
[14] Rakthanmanon, T., et al. Searching and Mining Trillions of
Time Series Subsequences under Dynamic Time Warping.
KDD 2012: 262-270.
[15] Sarangi, S. R. and Murthy, K. DUST: a generalized notion of
similarity between uncertain time series. KDD 2010: 383-392.
[16] Shokoohi-Yekta, M., Chen, Y., Campana, B., Hu, B., Zakaria,
J., and Keogh, E. Discovery of meaningful rules in time series.
KDD 2015: 1085-1094.
[17] Song, W.Z., Huang, R., Xu, M., Ma, A., Shirazi, B. and
LaHusen, R. Air-dropped sensor network for real-time highfidelity volcano monitoring. In Proceedings of the 7th
international conference on Mobile systems, applications, and
services, ACM, 2009, 305-318.
[18] Supporting Webpage:
https://sites.google.com/site/motifmissingdata/
[19] Wang, J., Cardell-Oliver, R. and Liu, W. An Incremental
Algorithm for Discovering Routine Behaviours from Smart
Meter Data. Knowledge-Based Systems, 113 (2016): 61-74.
[20] Veeraraghavan, A., Chellappa R., and Roy-Chowdhury, A.
The function space of an activity. Computer Vision and Pattern
Recognition. 1 (2006): 959-968.
[21] Weber, M. and Denk, M. Imputation of Cross-Country Time
Series: Techniques and Evaluation. In European Conference
on Quality in Official Statistics, 2010.
[22] Yeh, C.C.M. et al. Matrix Profile I: All Pairs Similarity Joins
for Time Series: A Unifying View that Includes Motifs,
Discords and Shapelets. ICDM 2016: 1317-1322.
[23] Yeh, M.Y., Wu, K. L., Yu, P. S. and Chen, M. S. PROUD: a
probabilistic approach to processing similarity queries over
uncertain data streams. EDBT 2009: 684-695
[24] Yi, X., Zheng, Y., Zhang, J. and Li, T. ST-MVL: Filling
Missing Values in Geo-sensory Time Series Data. IJCAI 2016:
2704-2710
[25] Yoon, C., O’Reilly, O., Bergen, K., and Beroza, G.
Earthquake Detection through Computationally Efficient
Similarity Search. Sci. Adv. 11.1(2015): e1501057.
[26] Zhu, Y. et al. Matrix Profile II: Exploiting a Novel Algorithm
and GPUs to Break the One Hundred Million Barrier for Time
Series Motifs and Joins, ICDM 2016: 739-748.

