Subsequence Matching in Data Streams

2.1 Related work

Similarity search over time-series data has been studied for many years. Most methods focus on similarity queries for static datasets with a query sequence. Specifically, given a query sequence and a distance threshold, a similarity query finds all the sequences or subsequences that are similar to the query sequence. The basic idea is to transform the sequence into the frequency domain using a discrete Fourier transform (DFT) and then to extract the few features from the resultant frequency-domain sequence [1]. Euclidean distance is used as a similarity measure. Other feature extraction functions include discrete wavelet transform (DWT) [2], singular value decomposition (SVD) [3], piecewise aggregate approximation (PAA) [4], and adaptive piecewise constant approximation (APCA) [5]. These functions have been proposed to reduce the number of dimensions of the time-series.

Since the Euclidean distance treats sequence elements independently, it cannot be used to calculate the distance between sequences whose lengths are different. DTW has been adopted to overcome these problems [6]; it is a widely used similarity measure that allows time-axis scaling. The DTW distance can be computed using dynamic programming techniques. Therefore, DTW is typically a much costlier approach than the Euclidean distance. To address the cost issue, several methods have been proposed that use a lower bound to refine the results and envelope techniques to constrain the computation cells [7]−[10].

Similarity search over data streams has recently attracted more research interest. In contrast to a setting where the data sequences are fixed, data streams arrive irregularly and are frequently updated. Therefore, incremental computation techniques that use the previous feature when computing the new feature are required. Several discriminative methods have been proposed for subsequence matching with a query sequence. These include subsequence matching based on prediction [11], approximate subsequence matching with data segmentation and piecewise line representation [12], subsequence matching supporting shifting and scaling in the time and amplitude domains [13], and subsequence matching based on DTW with batch filtering [14].

Another important area is pattern discovery between data streams. Mueen et al. [15] presented the first online motif discovery algorithm to accurately monitor and maintain motifs, which represent repeated subsequences in time-series, in real time. Papadimitriou et al. [16] proposed an algorithm for discovering optimal local patterns, which concisely describe the main trends in data streams.

In summary, there have been many previous studies on similarity search over data streams; however, none of them have addressed effective subsequence matching based on DTW. We formally define two problems for subsequence matching and present efficient and effective algorithms to solve them.

2.2 DTW

DTW is a transformation that allows sequences to be stretched along the time axis to minimize the distance between them. The DTW distance of two sequences is the sum of the tick-to-tick distances after the two sequences have been optimally warped to match each other. An illustration of DTW is shown in Fig. 1. The left figure is the alignment by DTW for measuring the DTW distance. To align two sequences, we construct a time warping matrix as shown in the right figure (the warping scope w will be described in section 4.1). The warping path is a set of grid cells in the time warping matrix; this set of grid cells represents the optimal alignment between the sequences. Consider two sequences, X = (x₁, x₂, ..., x_n) of length n and Y = (y₁, y₂, ..., y_m) of length m. Their DTW distance D(X, Y ) is defined as

 

 

Fig. 1. Illustration of DTW.

Note that ||x_i − y_j|| = (x_i - y_j)² is the distance between two numerical values in cell (i, j) of the time warping matrix. Note that other choices (e.g., absolute difference ||x_i − y_j|| = |x_i - y_j|) can also be used; the algorithms we present in this article are completely independent of the choice made. Specifically, DTW requires O(nm) time because the time warping matrix consists of nm elements. Note that the space complexity is O(m) because the algorithm needs only two columns (i.e., the current and previous columns) of the time warping matrix to compute the DTW distance.

In the problem of subsequence matching with a query sequence, given an evolving sequence X = (x₁, x₂, ..., x_n) and a fixed-length query sequence Y = (y₁, y₂, ..., y_m), we want to find the subsequences of X that are similar to Y in the sense of the DTW distance. By contrast, in the problem of subsequence matching between data streams, given two evolving sequences X = (x₁, x₂, ..., x_n) and Y = (y₁, y₂, ..., y_m), we want to find the subsequence pairs, namely the common local patterns over data streams. We will give the exact definitions for these problems later (in sections 3.1 and 4.1). In both settings, naïve ways of subsequence matching require unfeasible computation time and memory consumption in data stream processing. We show that our approaches offer considerable improvement without loss of accuracy.

3.1 Problem definition

Data stream X is a discrete, semi-infinite sequence of numbers x₁, x₂, …, x_n, …, where x_n is the most recent value. Note that n increases with every new time-tick. Let X [i_s : i_e] be the subsequence of X that starts from time-tick i_s and ends at i_e. The subsequence matching problem is to find the subsequence X [i_s : i_e] that is highly similar to a fixed-length query sequence Y (i.e., the subsequence with a small value of D(X [i_s : i_e], Y). However, a subtle point should be noted: whenever the query Y matches a subsequence of X (e.g., X [i_s : i_e]), we expect that there will be several other matches with subsequences that heavily overlap the local minimum best match. Overlaps provide the user with redundant information and would slow down the algorithm since all useless solutions are tracked and reported. In our solution, we discard all these extra matches. Specifically, the problem we want to solve is as follows:

Problem 1 Given a stream X (that is, an evolving data sequence, which at the time of interest has length n), a query sequence Y of fixed-length m, and a distance threshold e, report all subsequences X [i_s : i_e] such that

Hereafter we use the term optimal subsequences to refer to subsequences that satisfy both conditions.

3.2 SPRING

The most straightforward (and slowest) solution to this problem is to consider all possible subsequences X [i_s : i_e] (1 ≤ i_s ≤ i_e ≤ n) and apply the standard DTW dynamic programming algorithm, which requires O(n²) matrices. The time complexity would be O(n³m) (or O(n²m) per time-tick). This method is extremely expensive. Moreover, it cannot be extended to the streaming case. To solve the problem, we therefore proposed SPRING [17], which is an efficient algorithm. SPRING uses only a single matrix and detects optimal subsequences in stream processing.

3.2.1 Basic ideas

Our solution to the problem is based on two ideas. The first idea is star-padding, in which sequence Y is prefixed with a special value ("*") that always gives zero distance. This value stands for the don't care interval, that is, the interval (−∞ : +∞).

Definition 1 (Star-padding) Given a query sequence Y, star-padding of Y is defined as follows.

We use Y' to compute the DTW distances of Y and subsequences of X, instead of operating on the original sequence of Y.

Given a sequence Y = (y₁, y₂, ..., y_m), we have the star-padding of Y. Let X be a sequence of length n. We can then derive the minimum distance D(X [i_s : i_e], Y) from the matrix of X and Y'.

Star-padding dramatically reduces both time and space since we need to update only O(m) numbers per time-tick to derive the minimum distance.

Star-padding is a good first step, and it can tell us (a) where the subsequence match ends and (b) what the resulting distance is. However, such applications often also need the starting time-tick of the match. Our second idea is to use a subsequence time warping matrix (STWM): we augment the time warping matrix and record the starting position of each candidate subsequence.

Definition 2 (STWM) The STWM contains the value d(i, j), which is the best distance to match the prefix of length i from X with the prefix of length j from Y, and the starting position s(i, j) corresponding to d(i, j). The starting position is computed as

We obtain the starting position of D(X [i_s : i_e], Y) as

We update the starting position accompanied by the distance value as well as the distance value itself. By using the matrix, we can identify which subsequence gave the minimum distance during stream processing.

3.2.2 Algorithm

The way SPRING detects optimal subsequences is shown in Fig. 2. SPRING is carefully designed to (a) guarantee no false dismissals for the second condition of Problem 1 and (b) report each match as early as possible (detailed proofs are given in [17]). For each incoming data point x_n, we first incrementally update the m distance values d(n, j) and determine the m starting positions s(n, j). The algorithm reports the subsequence after confirming that the current optimal subsequence cannot be replaced by the upcoming subsequences. That is, we report the subsequence that gives the minimum distance d_min when the d(n, j) and s(n, j) arrays satisfy

which means that the captured optimal subsequence cannot be replaced by the upcoming subsequences. Otherwise, the upcoming candidate subsequences would not overlap the captured optimal subsequence. We initialize d_min and the d(n, j) arrays after the output.

Fig. 2. Illustration of SPRING algorithm, which uses only a single matrix to capture all qualifying sub­­se­quences.

Example 1 Assume that ε = 15, X = (5, 12, 6, 10, 6, 5, 13), and Y = (11, 6, 9, 4) in Fig. 3. The cell (i, j) of the matrix contains d(i, j) and s(i, j). At i = 3, we found candidate subsequence X[2 : 3] whose distance d(3, 4) = 14 below ε. At i = 4, although the distance d(4, 4) = 38 is larger than ε, we do not report X [2 : 3] since d(4, 3) = 2, which means X [2 : 3] can be replaced by the upcoming subsequences. We then capture the optimal subsequence X [2 : 5] at i = 5. X [2 : 5] is reported at i = 7 since we now know that none of the upcoming subsequences will be the optimal subsequence. Finally, because subsequences starting from i = 7 may be candidates for the next group, we do not initialize d(7, 1).

Fig. 3. Illustration of SPRING. The upper numbers show the distance in each element of the matrix. The numbers in parentheses show the starting position.

4.1 Problem definition

We have focused on finding subsequences similar to a query sequence. In this section, we address subsequence matching between data streams. That is, both sequences X and Y are co-evolving data streams, and we want to identify common local patterns between them.

Like data stream X, data stream Y is a discrete, semi-infinite sequence of numbers y₁, y₂, …, y_m, …, where y_m is the most recent value. Note that m increases with every new time-tick. Let Y [j_s : j_e] be the subsequence of Y that starts from time-tick j_s and ends at j_e. The lengths of X [i_s : i_e] and Y [j_s : j_e] are l_x = i_e − i_s + 1 and l_y = j_e − j_s + 1, respectively. The goal is to find the common local patterns of sequences by data stream processing based on DTW. That is, we want to detect subsequence pairs that satisfy

D(X [i_s : i_e],Y [j_s : j_e]) ≤ εL(l_x, l_y), (7)

where L is a function that sets the length of the subsequence. The DTW distance increases as the subsequence length increases since it is the sum of the distances between elements. Therefore, unlike the problem where sequence Y is fixed, the distance threshold should be proportional to the subsequence length to detect the subsequence pairs without depending on the subsequence length. Accordingly, we set the distance threshold at εL(l_x, l_y). In this article, the algorithm uses L(l_x, l_y) = (l_x+ l_y)/2, which is the average length of the two subsequences, but the user can make any other choice (e.g., L(l_x, l_y) = max(l_x, l_y) or L(l_x, l_y) = min(l_x, l_y)).

Equation (7) allows us to detect subsequence pairs without regard to the subsequence length. In practice, however, we might detect shorter and meaningless matching pairs owing to the influence of noise. We introduce the concept of subsequence match length to enable us to discard such meaningless pairs. We formally define the cross-similarity between X and Y, which indicates common local patterns.

Definition 3 (Cross-similarity) Given two sequences X and Y, a distance threshold ε, and a threshold of subsequence length l_min, X [i_s : i_e] and Y [j_s : j_e] have the property of cross-similarity if this sequence pair satisfies the condition

D(X [i_s : i_e],Y [j_s : j_e]) ≤ ε (L(l_x, l_y)− l_min). (8)

The minimum length l_min of subsequence matches should be given by the users. The subsequences that satisfy this equation are guaranteed to have lengths exceeding l_min.

In addition to subsequence matching with a query sequence, we need to consider several other matches that strongly overlap the local minimum best match. Specifically, in this setting, an overlap is simply where two subsequence pairs have a common alignment, which is defined as follows:

Definition 4 (Overlap) Given two warping paths for subsequence pairs of X and Y, their overlap is defined as the condition where the paths share at least one element.

Our solution for subsequence matching between data streams is to detect the local best subsequences from the set of overlapping subsequences. Thus, our goal is to find the best match of cross-similarity.

Problem 2 Given two sequences X and Y, thresholds ε, and l_min, report all subsequence pairs X [i_s : i_e] and Y [j_s : j_e] that satisfy the following conditions.

Hereafter we use the term qualifying subsequence pairs to refer to pairs that satisfy the first condition, and we use optimal subsequence pairs to refer to pairs that satisfy both conditions.

Typically, new elements in data streams are usually more significant than those in the distant past. To limit the number of cells in the matrix and focus on recent elements, we utilize a global constraint for DTW, namely the Sakoe-Chiba band [18] that restricts the warping path to the |i − j| ≤ w range (i.e., the gray cells in Fig. 1), where w is called the warping scope. In data stream processing, we compute the cells from a recent element to an earlier element of the warping scope w. If m = n, the warping scope is exactly equal to the Sakoe-Chiba band.

4.2 CrossMatch

To solve the problem of cross-similarity, the most straightforward solution is to consider all possible subsequences of X [i_s : i_e] (1 ≤ i_s < i_e ≤ n), and all possible subsequences of Y [j_s : j_e] (1 ≤ j_s < j_e ≤ m) in the warping scope and apply the standard DTW dynamic programming algorithm. Let w be the warping scope (the gray cells in the figure). The solution requires O(nw²+mw²) time (per update) and space because it has to handle a total of O(nw+mw) matrices to compute the DTW distance. If we prune dissimilar subsequence pairs and reduce the number of matrices, the distance computations become much more efficient. Our method, CrossMatch [19], finds good matches in a single matrix efficiently by pruning the subse-quences (see Fig. 4).

Fig. 4. Illustration of CrossMatch. Black cells indicate the warping paths of the optimal subsequence pairs, and gray cells indicate the warping scope.

4.2.1 Basic ideas

To identify the dissimilar subsequences early, we propose computing the DTW distance indirectly by using a scoring function. The scoring function computes the maximum cumulative score corresponding to the DTW distance with a score matrix. The score is determined by accumulating the difference between the threshold and the distance between the elements in the score matrix. Thus, we can recognize a dissimilar subsequence pair since the score has a negative value if the subsequence pair does not satisfy the first condition of Problem 2. The scoring function initializes the negative score to zero and then restarts the computation from the cell. This operation allows us to discard unqualifying, non-optimal subsequence pairs.

Definition 5 (Score matrix) Given two sequences, X = (x₁, x₂, ..., x_n) and Y = (y₁, y₂, ..., y_m), and the warping scope w, score V(X [i_s : i_e], Y [j_s : j_e]) of X [i_s : i_e] and Y [j_s : j_e] defined as follows:

Symbols b_v, b_h, and b_d in Eq. (9) indicate a weight function for each direction, which makes transformation between the score and the DTW distance possible. These values are determined as subsequence length L. For example, for L(l_x, l_y) = (l_x+l_y)/2, we obtain b_v = b_h =1/2 and b_d =1, respectively. The scoring function is designed so that the sum of the weights on the warping path is equal to subsequence length L. Therefore, it guarantees transformation between the DTW distance and the score, and it finds the qualifying subsequence pairs without any omissions (detailed proofs are given in [19]). The DTW distance of a subsequence pair is computed from the score and the subsequence length as

Equation (10) holds for the time warping and the score matrices, which have the same starting position (i_s, j_e).

The scoring function tells us the subsequence match ends and the subsequence score. However, we lose the information about the starting position of the subsequence. This is the motivation behind our second idea, a position matrix: we store the starting position to keep track of qualifying subsequence pairs in a streaming fashion.

Definition 6 (Position matrix) The position matrix stores the starting position of each subsequence pair. The starting position p(i, j) corresponding to score v(i, j) is computed as:

The starting position is described as a coordinate value; p(i_e, j_e) indicates the starting position (i_s, j_s) of the subsequence pair X [i_s : i_e] and Y [j_s : j_e]. We can identify the optimal subsequence that gives the maximum score during stream processing by using the score and position matrices. Moreover, the starting position of the shared cell is maintained through the subsequent alignments because we repeat the operation, which maintains the starting position of the selected previous cell. Thus, we know the overlapping subsequence pairs from the fact that the starting positions match.

4.2.2 Algorithm

We now have all the pieces needed to answer the question: how do we find the optimal subsequence pairs? Every time x_n or y_m is received at time-tick n or m, our CrossMatch algorithm incrementally updates the score and starting position and retains the end position. We use a candidate array to find the optimal subsequence pair, and we store the best pair in a set of overlapping subsequence pairs. CrossMatch reports the optimal subsequence pair after confirming that it cannot be replaced by the upcoming subsequence pairs. The upcoming candidate subsequence pairs do not overlap the captured optimal subsequence pair if the starting positions in the position matrix satisfy the condition:

CrossMatch requires only O(w) (i.e., constant) time (per update) and space. This is a great reduction because the straightforward solution requires O(nw²+mw²) time (per update) and space.

Example 2 Assume that we have two sequences X = (5, 12, 6, 10, 3, 18) and Y = (11, 9, 4, 2, 9, 13) as well as ε = 14, l_min = 2, and w = 3 in Fig. 5. To simplify our example with no loss of generality, we assume that x_i and y_j arrive alternately. At each time-tick, the algorithm updates the scores and the starting positions. At i = 4, we update the cells from (4, 1) to (4, 3) and identify a candidate subsequence, X[2 : 4] and Y[1 : 3], starting at (2, 1), whose score v(4, 3) = 36 is greater than εl_min. At j = 4, we update the cells from (1, 4) to (4, 4). Although we detect no subsequences satisfying the condition, we do not report the subsequence X[2 : 4] & Y[1 : 3] since this pair might be replaced by upcoming subsequences. We then capture the optimal subsequence pair X[2 : 5] & Y[1 : 4] at i = 5. Finally, we report the subsequence as the optimal subsequence at j = 6 since we can confirm that none of the upcoming subsequences can be optimal.

Fig. 5. Example of cross-similarity detection. The lightly shaded cells signify cross-similarity, and the dark cell in each matrix shows the best match.