In this work, we are dealing with time series. A time series X = ( x 1 , x 2 , … , x N ) is an ordered sequence of real-valued numbers, where N = | X | is the length of X. A subsequence S, S ′ or X [ i , j ] = ( x i , x i + 1 , … , x j ) ( 1 ⩽ i ⩽ j ⩽ n ) denotes the continuous sequence of length j − i + 1 starting from the i-th position in X. Note that the subsequence is itself a time series.

Given a time series X, a query sequence Q, and a distance function D, the problem of subsequence matching is to find out the top-K subsequences from X, denoted as R = { S 1 , S 2 , … , S K }, which are most similar to Q.

The two representative distance measures are Euclidean Distance (ED) and Dynamic Time Warping (DTW). Formally, given two length-L sequences, S and S ′ , their ED and DTW distance can be computed as the following:

In this paper, instead of processing one single query, we attempt to find out subsequences similar to multiple queries. That is, given a set of queries, Q = { Q 1 , Q 2 , … , Q N }, our objective is to find out the top-K subsequences similar to queries in Q, denoted as R.

Since each query sequence varies in length, we do not impose a constraint on the length of subsequences in R. In this way, we find out variable-length subsequences answering multiple queries, which is worthy of wide use in real time series applications.

In this paper, instead of processing the queries in Q independently, we first represent them by a unified formation, which is a multi-dimensional probability distribution. Then we find target subsequences from X based on the representation.

In many real world applications, the meaningful query sequence can be approximately represented as a sequence of line segments. Recall the EOG pattern in Fig. 1. The query sequence Q can be approximated with 4 line segments. These line segments capture the most representative characteristics in Q.

In response, we propose a two-step approach to represent the bundle of queries together. In the first step, we represent each single query Q i in Q individually by a traditional segmentation in which each segment is described by some features. Then, in the second step, we represent each feature as a Gaussian distribution over the values from multiple queries.

Given the fact that the query representation consists of 4 m number of Gaussian distributions rather than a sequence, the existing distance measures, like ED and DTW, are inapplicable. In this paper, we propose a novel distance function D ( S , Q ) based on the probability distribution.

Formally, to define the distance between subsequence S and the query representation P Q , we first approximate S with m line segments. We indicate the segmentation as s e g = ( S 1 , S 2 , … , S m ), where S j ( 1 ⩽ j ⩽ m) denotes the j-th segment of S. Note that the segment here infers subsequence rather than line segment. We extract four features, the length l j , the slope θ j , the value of the starting point v j , and the MSE error ε j from the linear representation of S j . For ease of presentation, we represent all the j-th segments in Q as Q j . That is, Q j = ( Q 1 j , Q 2 j , … , Q N j ). Then the distance between S j and Q j is (2) dist ( S j , Q j ) = − log ( p l j ( l j ) p θ j ( θ j ) p v j ( v j ) p ε j ( ε j ) ) , dist ( S j , Q j ) is the negative logarithm of the probability. The smaller the value, the more similar S j and Q j are. Accordingly, under segmentation s e g, the distance between S and the query set Q can be computed as (3) D ( S , Q , seg ) = ∑ j = 1 m dist ( S j , Q j ) . Since S can be segmented by different segmentations, the value of D ( S , Q , seg ) may be different. In this paper, we define the distance between S and Q as the minimal one among all possible segmentations, that is, (4) D ( S , Q ) = min seg D ( S , Q , seg ) .

We first introduce our search strategy to generate the candidate set CS with size not exceeding the parameter n c . We utilize an iterative approach, and generate the candidates segment-by-segment. In the first round, we generate at most n c number of candidates with only one segment. The candidate set is denoted as CS 1 = { c s 1 , c s 2 , … , c s n c }. Each candidate, c s i , is a triple ⟨ s i , e i , d i ⟩, in which s i is its starting point and e i is its ending point. So c s i corresponds to the subsequence X [ s i , e i ]. The third element d i is distance dist ( c s 1 , Q 1 ). All candidates in CS 1 are ordered based on the values of d i ascendingly. In other words, CS 1 contains n c number of subsequences with smallest distance with Q 1 . We discuss how to select top- n c candidates in the next section.

In the second round, we obtain candidate set CS 2 by extending the candidates in CS 1 with the second segment. Specifically, given any candidate subsequence c s = ⟨ s , e , d ⟩ in CS 1 , if we want to extend c s to c s ′ with a length-L segment, the new candidate c s ′ = ⟨ s , e ′ , d ′ ⟩ contains two segments: one corresponds to X [ s , e ] and the other to X [ e + 1 , e + L ]. Note that the new starting point s keeps unchanged, and the new ending point e ′ changes to e + L. Also, the new distance d ′ is updated to d + dist ( X [ e + 1 , e ′ ] , Q 2 ). Since from each candidate c s in CS 1 , we can extend it to multiple candidates by concatenating X [ s , e ] with variable length segments, we generate all possible candidates, compute the distance and add top- n c candidates into CS 2 .

After m rounds, we obtain the candidate set CS m , which is the final candidate set CS. Now each candidate in CS consists of m segments.

Now we introduce our approach to generate possible candidates in each round.

In the first round, we enumerate all subsequences X [ s , e ] from all possible starting points, that is, we try all s’s within [ 1 , n ]. To avoid the low-quality candidates, we only select subsequences with length satisfying the 3 σ standard. Formally, given any starting point s, the ending point e must satisfy e ∈ [ μ l 1 − 3 σ l 1 , μ l 1 + 3 σ l 1 ]. For each candidate subsequence X [ s , e ], we compute the optimal linear approximation, y = θ · x + b, as well as the MSE error ε as follows, (5) θ = 12 ∑ i x − 6 ( l + 1 ) ∑ x l ( l + 1 ) ( l − 1 ) , b = 6 ∑ i x − 2 ( 2 l + 1 ) ∑ x l ( 1 − l ) , ε = ∑ x 2 + θ 2 ∑ i 2 + l b 2 − 2 θ ∑ i x − 2 b ∑ x + 2 θ b ∑ i , where l = e − s + 1 is the length of X [ s , e ]. After that, we compute dist ( X [ s , e ] , Q 1 ). Also, if only any other feature of a candidate (θ, v or ε) violates the 3 σ standard, we ignore this candidate. During enumerating different candidates, we maintain CS 1 as a priority queue to keep top- n c candidates.

In the j-th round ( 2 ⩽ j ⩽ m), we generate candidates by extending previous ones in CS j − 1 . For each candidate c s = ⟨ s , e , d ⟩ in CS j − 1 , we try all possible segments next to X [ s , e ] whose length falls within [ μ l j − 3 σ l j , μ l j + 3 σ l j ]. Similar with the first round, we also dismiss the candidates violating the 3 σ standard. When extending candidate c s = ⟨ s , e , d ⟩ to c s ′ = ⟨ s , e ′ , d ′ ⟩, except the new ending point e ′ , we also update the distance as d ′ = d + dist ( X [ e + 1 , e ′ ] , Q j ).

Note that the subsequences in C S are not approximated optimally. Here, an additional refinement step has to be performed, where subsequences are verified and re-ordered using the optimal segmentation. Specifically, we fetch the subsequences in CS, approximate each subsequence c s with m line segments via a dynamic programming algorithm, and thus get the actual distance D ( c s , Q ) under the new segmentation.

The objective of the segmentation is to minimize the distance between c s and Q. We search the optimal segmentation from left to right sequentially on c s. We define E ( i , j ) ( 1 ⩽ i ⩽ m and 1 ⩽ j ⩽ | c s |) to be the minimal distance between the prefix of c s (i.e. c s = X [ 1 , j ]) and the prefix of Q with i segments (i.e. [ Q 1 , Q 2 , … , Q i ]). We begin by initializing E ( 1 , j ) to be dist ( X [ 1 , j ] , Q 1 ). When computing E ( i , j ) ( 2 ⩽ i ⩽ m ), we consider all the possible segmentations of X [ 1 , k ] ( i ⩽ k ⩽ j ) with i − 1 segments, compare the sum of E ( i − 1 , k ) and dist ( X [ k , j ] , Q k ), and define E ( i , j ) to be the minimal one. Formally, the dynamic programming equation is presented as the following (6) E ( i , j ) = dist ( X [ 1 , j ] , Q 1 ) , i = 1 , + ∞ , i > j , min i ⩽ k ⩽ j ( E ( i − 1 , k ) + dist ( X [ k , j ] , Q k ) ) , otherwise . We re-compute the distance between each subsequence in CS and Q. Before re-ordering, we have to dismiss the trivial match subsequences since candidates can overlap with each other. Formally, given two candidates c s = X [ s , e ] and c s ′ = X [ s ′ , e ′ ], if their overlapping ratio, c s ∩ c s ′ c s ∪ c s ′ , exceeds min ( 0.8 , min ( | Q | ) max ( | Q | ) ), where the latter indicates the ratio between the minimum length and the maximum length of queries Q, we dismiss the subsequence with the larger distance.

Afterwards, we simply sort the remaining subsequences in CS according to D ( c s , Q ) and select the K smallest as the final results R.

In this section, we propose two optimization strategies to further accelerate the search process.

The overall process of our approach can be divided into two phases: query representation and query processing.

Before the two phases, we have to scan the whole time series X once to store the basic aggregates. Its time cost is O ( n ), where n is the length of X.

Then, in the first phase, we scan all the queries and perform traditional piecewise linear approximations. Thanks to the basic aggregates, the time cost of conducting linear regression is negligible. This phase is then O ( N | Q | 2 ) in time complexity, where N is the number of queries. For ease of presentation, although multiple queries can vary in length, we denote | Q | as the length of the queries.

In the second phase, we first assume w i = 6 σ l i + 1, where σ l i is the standard deviation of the length values of ( H 1 i , H 2 i , … , H N i ). Suppose we generate candidates from left to right sequentially; when considering the first segments, it requires O ( w 1 n log ( n c ) ) time to maintain the candidate set. For the following i-th segments, the time complexity is O ( w i n c log ( n c ) ). In the post-processing process, we have to verify and compute the actual distance between Q and every subsequence c s in the candidate set. It takes O ( n c m | c s | 2 ) time to conduct the dynamic programming algorithm for the n c candidates, where m is the number of line segments. Moreover, the re-ordering process is at most O ( n c log n c ) in time complexity.

Since w i is a constant, n c is proportional to n, and | c s | and | Q | is much smaller than n, we can infer that the time complexity of our approach is O ( n log ( n c ) ) theoretically.

Note that it is difficult to find reasonable baselines to compare our approach to because the existing methods are only devised for the case of a single query. Since no approach can deal with multiple queries as a whole, we choose two representative subsequence matching algorithms, UCR Suite (Rakthanmanon et al., 2012) and SpADe (Chen et al., 2007), and then enable them to handle the problem of multiple queries. UCR Suite finds the best normalized subsequences and supports both Euclidean distance and Dynamic Time Warping (UCR-ED and UCR-DTW for short). SpADe finds the shortest path within several local patterns, and is able to handle shifting and scaling both in temporal and amplitude dimensions.

To make UCR-Suite (or SpADe) support multiple queries, we first find out top-K similar subsequences for each query based on UCR-Suite (or SpADe). We then sort these K · N number of subsequences in the ascending order of their distances (normalized by the subsequence length) and pick out the top-K ones excluding any trivial results as final. For UCR-Suite, we utilize both ED and DTW as the distance, denoted as UCR-ED and UCR-DTW, respectively.

Let N be the number of queries in Q, we denote our approach and the other three competitors as MQ-N, UCR-ED-N, UCR-DTW-N and SpADe-N, respectively. Note that the three rivals have to scan the time series for N times. For fairness, we do not count I/O time when comparing efficiency.

In the first experiment, we compare our approach MQ with UCR-ED, UCR-DTW, and SpADe on synthetic datasets. Both accuracy and efficiency are tested. To compare these approaches extensively, we vary the parameter r for all three variations, length, amplitude and shape. Specifically, for the length variation, we set the maximal scaling factor r from 0.25 to 0.7 with step 0.05. For the amplitude variation, we varied r in a range of 0.2 to 2 with step 0.2. For the shape variation, we changed r in a range of 0.05 to 0.5 with step 0.05.

We set the number of queries, N, as 5 and 10, respectively. The experimental results are then indicated by MQ-5, MQ-10, UCR-DTW-5, UCR-DTW-10, SpADe-5, and SpADe-10.1¹

UCR-ED is omitted for its consistent inferiority to UCR-DTW.

In each set of experiments, we attempt to find out the top-40 subsequences in the time series. The accuracy is the ratio between the number of correct subsequences and 40. Subsequence S is correct, if the overlapping ratio between S and certain planted instance exceeds the tolerance parameter, ϵ = min ( 0.8 , min ( | Q | ) max ( | Q | ) ).

The results are shown in Fig. 4 and Fig. 5, respectively. It can be seen that MQ outperforms UCR-DTW and SpADe in both accuracy and efficiency under all variations. The reason is that MQ summarizes common characteristics in multiple queries while the other approaches are only able to find out the subsequences that are similar to certain given query. As a consequence, when the number of query sequences N increases, UCR-DTW and SpADe yield better results. Nevertheless, UCR-DTW-10 (or SpADe-10) is twice as slow as UCR-DTW-5 (or SpADe-5) on average, which means with more query sequences provided, UCR-DTW and SpADe can find out more satisfying subsequences, but at the cost of efficiency. Instead, MQ captures latent semantics and thus demonstrates its superiority in both accuracy and efficiency under different variations. The only exception is in Fig. 4(c), where the accuracy of MQ sharply decreases when r = 0.4, which is mainly because the undue scaling in shape will result in ambiguity of what users really want.

In this experiment, we compare our approach with other ones as the number of query sequences N varied on real dataset. We pick out queries from the query set of size 15 so that their lengths distribute uniformly. Note that no matter the value of N, we find out top-100 subsequences from the real time series. The results are shown in Fig. 6.

It can be seen that in Fig. 6(a), the accuracy of MQ, UCR-ED, and UCR-DTW increases when N increases while the accuracy of SpADe is consistently low. The reason is that SpADe extracts local patterns by using a fixed size of sliding window, and consequently, it fails to capture the query intuition. Moreover, it is noteworthy that when N = 6, the accuracy of MQ has already exceeded 0.9, which means that MQ is able to find out what users really want with only a small query set.

Figure 6(b) compares MQ and other approaches on efficiency. Obviously, MQ outperforms all other approaches, and its running time is insensitive to N since MQ summarizes all the queries and searches for results in the time series only once. Due to the specific pattern shape and the lack of abundant pruning strategies, UCR-ED is inferior to UCR-DTW in terms of efficiency in this dataset.

In this experiment, we investigate the influence of the parameter n c , the size of the candidate set. On the real dataset, we varied n c from 0.01 to 0.1, and the number of queries, N was set to 3, 6, and 9 respectively.

Results are shown in Fig. 7. It can be seen that as n c gets larger, the accuracy of MQ increases while the efficiency decreases. It is because once the query set is fixed, we can maintain more candidates by enlarging the candidate set, i.e. increasing the value of n c at the cost of more running time. Also, we can find that MQ has already achieved satisfying results when n c = 0.05, so we set the default value of n c to 0.05 in previous experiments. Generally, as shown in Fig. 7(a), the accuracy of MQ increases as the number of queries N increases. The only exception is when n c = 0.01. The reason is that the standard deviation of the query feature values experiences an increase as N increases, resulting in a larger search space. The small candidate set then fails to maintain enough candidates, and thus achieves poor results. Meanwhile, the changes in the search space cause the running time to increase proportionally, as shown in Fig. 7(b).