Multi-Querying: A Subsequence Matching Approach to Support Multiple Queries

3.1Definition

In this work, we are dealing with time series. A time series X=(x1,x2,…,xN) is an ordered sequence of real-valued numbers, where N=|X| is the length of X. A subsequence S, S′ or X[i,j]=(xi,xi+1,…,xj) (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={S1,S2,…,SK}, 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={Q1,Q2,…,QN}, 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.

4.1Query Representation

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 Qi 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.

4.2Distance Definition

Given the fact that the query representation consists of 4m 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 PQ, we first approximate S with m line segments. We indicate the segmentation as seg=(S1,S2,…,Sm), where Sj (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 lj, the slope θj, the value of the starting point vj, and the MSE error εj from the linear representation of Sj. For ease of presentation, we represent all the j-th segments in Q as Qj. That is, Qj=(Q1j,Q2j,…,QNj). Then the distance between Sj and Qj is

5.1BFS-Based Search Process

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

In the second round, we obtain candidate set CS2 by extending the candidates in CS1 with the second segment. Specifically, given any candidate subsequence cs=⟨s,e,d⟩ in CS1, if we want to extend cs to cs′ with a length-L segment, the new candidate cs′=⟨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′],Q2). Since from each candidate cs in CS1, 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-nc candidates into CS2.

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

5.2Candidate Generation

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∈[μl1−3σl1,μl1+3σl1]. For each candidate subsequence X[s,e], we compute the optimal linear approximation, y=θ·x+b, as well as the MSE error ε as follows,

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

5.3Post-Processing

Note that the subsequences in CS 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 cs with m line segments via a dynamic programming algorithm, and thus get the actual distance D(cs,Q) under the new segmentation.

The objective of the segmentation is to minimize the distance between cs and Q. We search the optimal segmentation from left to right sequentially on cs. We define E(i,j) (1⩽i⩽m and 1⩽j⩽|cs|) to be the minimal distance between the prefix of cs (i.e. cs=X[1,j]) and the prefix of Q with i segments (i.e. [Q1,Q2,…,Qi]). We begin by initializing E(1,j) to be dist(X[1,j],Q1). 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],Qk), and define E(i,j) to be the minimal one. Formally, the dynamic programming equation is presented as the following

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

5.4Optimization

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

5.5Complexity Analysis

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 wi=6σli+1, where σli is the standard deviation of the length values of (H1i,H2i,…,HNi). Suppose we generate candidates from left to right sequentially; when considering the first segments, it requires O(w1nlog(nc)) time to maintain the candidate set. For the following i-th segments, the time complexity is O(winclog(nc)). In the post-processing process, we have to verify and compute the actual distance between Q and every subsequence cs in the candidate set. It takes O(ncm|cs|2) time to conduct the dynamic programming algorithm for the nc candidates, where m is the number of line segments. Moreover, the re-ordering process is at most O(nclognc) in time complexity.

Since wi is a constant, nc is proportional to n, and |cs| and |Q| is much smaller than n, we can infer that the time complexity of our approach is O(nlog(nc)) theoretically.

6.1Datasets

6.1.1Synthetic Datasets

We generate the synthetic sequence as follows. First, we generate a long random walk sequence T. Then we embed in T some meaningful pattern instances, in which some are taken as queries, and others as target results.

Fig. 3

Two fundamental patterns.

(a) W-wave

(b) Backward wave

We generate two types of patterns, W-wave from the Adiac dataset in UCR archive (Dau et al., 2019), and backward wave, common in stock series. For each pattern, we first generate a seed instance and add some noise following a Gaussian distribution N(0,0.05), as shown in Fig. 3. Then we modify them to generate more instances. Specifically, we adopt three types of variations, length, amplitude, and shape as follows.

• ∙ Length: It is obvious that the W-wave in Fig. 3(a) can be split into four segments. For the middle two segments, we change the segment length with a scaling factor λ by inserting or deleting some data points. In other words, for a length-l segment, the length of the new segment is l·λ.

• ∙ Amplitude: Similar to the case of length scaling, we increase or decrease the amplitude of the middle two segments in the W-wave. We still use a factor λ to control the extent. Specifically, every value in new segment changes to v·λ, where v is the original value.

• ∙ Shape: For the backward wave shown in Fig. 3(b), we change the global shape of the pattern by modifying the last two segments on both length and amplitude. To be more specific, we change the segment length from l to l·λ, and the data values from v to v·λ.

Note that the three types of variations will influence different features in the line segments. We list them respectively in Table 1.

Table 1

Influence of variations.

	l	θ	v	ε
Length	√			√
Amplitude		√	√	√
Shape	√		√	√

For each type of variation, we test our approach under different extents. Take the length variation as an example. To generate a dataset, we first set a parameter r, which determines the length variation range. Specifically, given r, we can only set the length scaling factor λ within the range [1/(1+r),1+r]. Obviously, the larger the value of r, the larger the length scaling extent.

For each variation, given the fixed r, we pick out 50 values from the range [1/(1+r),1+r] as the value of λ. Then we generate 50 corresponding pattern instances, 40 of which are randomly planted into the random walk long time series. The rest 10 instances form the query set Q. A summary of the synthetic datasets is presented in Table 2.

Table 2

A summary of the synthetic datasets.

Dataset	Pattern	Variation	Length
1	W-wave	Length	0.7 million
2	W-wave	Amplitude	0.7 million
3	Backward wave	Shape	0.4 million

6.2Counterpart Approaches

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.

Fig. 4

Accuracy comparisons under different variations.

(a) Length scaling

(b) Amplitude scaling

(c) Shape scaling

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.

6.3Results on Synthetic Datasets

Fig. 5

Efficiency comparisons under different variations.

(a) Length scaling

(b) Amplitude scaling

(c) Shape scaling

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.11 For MQ, we set the only parameter, the size of the candidate set nc, to be 0.05·n, where n is the length of the time series X.

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.

6.4Results on Real Datasets

Fig. 6

Comparison on real dataset.

(a) Accuracy vs. N

(b) Efficiency vs. N

Fig. 7

Influence of parameter nc.

(a) Accuracy vs. nc

(b) Efficiency vs. nc

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.

6.5Influence of Parameter nc

In this experiment, we investigate the influence of the parameter nc, the size of the candidate set. On the real dataset, we varied nc 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 nc 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 nc at the cost of more running time. Also, we can find that MQ has already achieved satisfying results when nc=0.05, so we set the default value of nc 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 nc=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).