Time-series alignment is an important question in many different applications, including bioinformatics [1], computer vision [2], speech recognition [3], speech synthesis [4,5] and human action recognition [6]. The most commonly used method to perform this alignment is dynamic time warping (DTW), which uses the Bellman's dynamic programming technique to search for the optimal time warping [7,8]. However, it is difficult to apply DTW over large sets of time-series because the computational complexity 1 of DTW grows exponentially with the number of time-series to be aligned [9,10].
In this paper, we introduce a new time warping algorithm, called trainable time warping (TTW) with linear computational complexity in both N and T , where N and T are the number and the length of time-series, respectively 2 . TTW first aligns the input sequences and then takes an average over the synchronized sequences to calculate the centroid signal. We propose a sinc convolutional kernel that is able to warp signals in a nonlinear (elastic) manner. TTW uses this convolutional kernel along with a gradient-based optimization technique to synchronize multiple time-series.
DTW extensions for multiple time-series include a progressive DTW-based method, called NonLinear Alignment and Averaging Filters (NLAAF) [11]. In each iteration of NLAAF, time-series are randomly grouped into a set of pairs and then each pair is fused using the DTW averaging. This procedure is repeated until the centroid time-series is achieved. However, although NLAAF is practical for averaging multiple sequences, it has two problems: (1) orderingthe centroid is sensitive to the order of pairing time-series and (2) error propagation -there is no feedback in the algorithm to reduce the errors that occur at early steps.
Many papers presented methods to address these two weaknesses. Niennattrakul et al. [12], addressed the first problem, introducing a method that uses hierarchical clustering to automatically create a reasonable order for pairing time-series. Petitjean et al. [13,14], addressed both, introducing an averaging method, DTW Barycenter Averaging (DBA), which iterates over two steps:
(1) aligning: DBA computes a DTW mapping between the current centroid and each input time-series and (2) updating: DBA uses the mappings obtained from the previous step to align the time-series and update the centroid sequence. However, DBA is sensitive to the quality of the initial average signal. When the initial average is considerably different from the best solution, DBA is highly likely to converge to a weak local minimum [15]. Extensions have been proposed to improve the performance of the DBA algorithm (e.g., soft-DTW [15], stochastic subgradient (SSG) [8] and constraint DTW [7]). However, all have the computational complexity of O(N T 2 ). The complexity is quadratic in T , which is problematic in many applications that must process long-term time-series [16].
To the best of our knowledge, generalized time warping (GTW) [17,18] is the only algorithm that is able to align multiple signals with linear complexity in T . GTW approximates the optimal temporal warping by linearly combining a fixed set of monotonic basis functions. They introduced a Gauss-Newton-based procedure to learn the weights of the basis functions. However, in the cases where the temporal relationship between the time-series is complex, GTW requires a large number of complex basis functions to be effective; defining these basis functions is very difficult.
We evaluate the efficacy of TTW on 85 datasets of the UCR archive [19] for two applications: DTW averaging and time-series classification. We implement GTW as the baseline system, because GTW has similar computational complexity. The results show that TTW is able to outperform GTW on most of the datasets.
In this section, we first describe how the shifted sinc kernel can be used to introduce a learnable uniform shift to a time-series. We then extend to a more general non-linear time warping case. Finally, we introduce the TTW algorithm.
We first consider a simple case where the time alignment between the signals can be performed by applying a uniform delay τn (i.e., xn[t] = xn[t -τn]). We can implement this alignment by convolving the input signal xn with a shifted dirac-delta, δ[t -τn], kernel:
where * is the convolution operation.
In order to perform temporal alignment, we must estimate the best value for the shift parameter τn, which can be done through a search algorithm [22]. However, search algorithms are inefficient when there are multiple input signals, because the search space grows exponentially with the number of the signals (similar to the standard DTW). Gradient-based optimization techniques are commonly used to deal with very large search spaces. However, τn is not learnable in this way because: (1) τn is an integer, while gradientbased optimization strategies have been designed for a continuous search space and (2) δ(t) is not differentiable whenever t = 0.
As a solution, we propose to transfer the time-series to the continuous-time domain and perform the temporal alignment in that domain. In this approach, τn is applied to a continuous time-series and can be any real number. To implement this idea, three main steps are needed: (1) Interpolation -In the first step, we use sinc interpolation to transfer to the continuous-time domain. Ideal sinc interpolation preserves all frequency components of the input signals [23]. The relationship between discrete and continuous time signals in the sinc interpolation technique is 3 :
where xn[m] is the m-th sample of the n-th signal, xn(t) is an interpolation of the n-th signal at time t, Ts is the sampling rate and the sinc function is defined as:
Shifting -In the second step, we temporally align the continuoustime signals obtained from the previous step. In the simplest case, the alignment is performed by applying a fixed shift, i.e.,
In this equation, the size of the shift is τn samples or τnTs seconds. We can simplify Equation (5) using properties of the sinc function:
3 In this paper, we use [.] for descrete-time signals and (.) for continuoustime signals. For example, f [m] is the m-th sample of a descrete-time signal f , and f (t) is the value of a continuous-time signal, f , at time t.
This equation shows that we can use a shifted low-pass (sinc) filter to introduce a fixed shift to the interpolated signals.
Sampling -Finally, in the third step, we transfer the shifted signals to the original discrete-time domain using uniform sampling. Uniform sampling is performed by replacing t with tTs in Equation (6) which results in:
This equation is the final expression for the shifted sinc kernel. In the case of integer τn, the kernel (Equation ( 7)) is equal to the dirac-delta function in the discrete-time domain (Equation ( 2)), and therefore the sinc kernel just introduces a fixed shift. In other cases, the kernel interpolates the signal to continuous-time domain, introduces a shift there and returns the signal to the discrete-time domain.
Equation (7) shows that introducing a shift in the continuoustime domain is equivalent to convolving the signals with a shifted sinc kernel in the discrete-time domain. Therefore, to align the signals with a simple shift, we just need to convolve them with the shifted sinc kernels (Equation ( 7)) and learn the shift parameters τn. The sinc kernel gives us the ability to introduce real-valued shifts to discrete-time signals. It is also a differentiable function and can therefore be used in gradient-based optimization techniques.
The shifted sinc kernel is able to apply a fixed shift; however, a fixed shift is generally not sufficient for aligning time-series. In DTWbased algorithms, alignment is performed by applying an elastic (non-linear) transform to the time axis of the signals. This transform is called time warping or time alignment function. Inspired by these algorithms, we extend our linear transform (t -→ t -τn), to a more general nonlinear transform (t -→ τn[t]) by replacing the term t -τn with τn[t] in Equation (7):
We represent the set of all warping functions with T = {τ1, . . . , τN }, where τn = [τn [1] 8) is computational expensive. For a signal of length T , it needs O(T 2 ) arithmetic operations, which is problematic when the algorithm processes long-term signals. In order to overcome this problem, we approximate the ideal sinc with a windowed sinc function. We know that the sinc function is a sine wave that decays by increasing the time; more than 99% of its power lies in the range of [-10, 10]. Therefore, we assume sinc(t) = 0 whenever |t| > 10.
Using this assumption we can approximate Equation (8) by:
where . denotes the floor operator. In other words, we use a truncated sinc (windowed sinc with a rectangular window) to reduce the computational complexity of Equation ( 9) from O(T 2 ) to O(T ).
To accomplish DTW averaging, we first synchronize the given timeseries through an optimal set of time warping functions, T ; we then take average of the synchronized sequences. We estimate the optimal set of time warping functions by numerically solving the following optimization problem using the Adam algorithm [24]:
where D( X) is the within-group mean square error (MSE) of the synchronized signals, X, and y[t] is the t-th sample of the average signal. According to this Equation, we find the optimal alignments that satisfy the DTW conditions and minimizes the within-group distance of the synchronized signals. The DTW conditions guarantee that the synchronized sequences preserve the shape of the original time-series. The next section will explain these conditions in detail.
In this section, we ensure that our time warping algorithm satisfies the DTW constraints including continuity, monotonicity and boundary conditions [17]. This mitigates concerns regarding unreasonable choices of T , which could otherwise change the shape of the input time-series in an unsatisfactory manner. Continuity: Sudden changes in the warping function distort the overall shape of the input time-series [17,7]. In order to avoid sudden changes, we enforce a smoothness over τn [t]. There are different methods for generating a smooth sequence. In this paper, we propose to estimate the first K coefficients of the discrete sine transform (DST) of the warping functions instead of estimating all samples [25]. This guarantees the smoothness of the alignment by eliminating high frequency components. In other words, we assume that the warping functions can be written as:
where an = [a n 1 , . . . , a n K ] is the DST coefficients of τn[t]. In the alignment process, we find these coefficients instead of τn[t]. The parameter K (number of sine components) controls the smoothness of the warping function such that with a smaller value of K we generate a smoother warping function. By increasing the value of K, we can capture higher-frequency variations of the input signals. However, our experiments show that TTW with a higher value of K is more likely to converge to a weak local minimum. In addition, The DST expressed by Equation (11) significantly reduces the number of the learning parameters, which enables us to leverage more complex (e.g., quasi-Newton) optimization algorithms.
Monotonicity: Time warping algorithms must preserve the order of the samples. This constraint guarantees that the time warping algorithm provides consistent temporal changes between input and output sequences. The warping function must be monotonically increasing to satisfy this constraint [17] (i.e., for all t ∈ {2, . . . , T }, ∆τn[t] ≥ 0). In the TTW algorithm, we enforce this constraint in our optimization process; after each iteration (update) of the Adam algorithm, we traverse our warping functions from beginning to end; if a sample does not follow the monotonicity condition, we set the sample to its previous sample; more precisely, for all t ∈ {2, . . . , T } and ∆τn[t] ≥ 0, we set τn[t] to τn[t -1].
Boundary conditions: τn[1] = 1 and τn[T ] = T . This constraint (along with the other constraints) guarantees that the aligned signal contains all parts of the original signal [26]. Our algorithm is guaranteed to satisfy this constraint because the DST transform expressed by Equation ( 11) generates 1 and T for the starting (τn [1]) and the ending (τn[T ]) points of the warping functions.
The goal is to synchronize a set of N time-series of length T , X = {x1, . . . , xN } where xn ∈ R T , and fuse them into a single centroid time-series y ∈ R T . To do so, we define K learning parameters, an = {a n 1 , . . . , a n K }, for each input signal xn and initialize them with zero. TTW finds the optimal set of learning parameters through iterating over the following forward, backward and update steps.
Forward: takes the given signals X and current values of the learning parameters, A = {a1, . . . , aN } as the input to compute the current average sequence. In the forward step, TTW generates T using Equation (11); fixes the monotonicity problems of T using the algorithm explained in Section 3.4 (Monotonicity); generates X using Equation ( 9); and calculates the within-group distance D and the current mean signal y using Equation (10).
Backward: takes all intermediate signals of the forward step to generate the derivatives of the within-group distance D( X) with respect to the learning parameters A, (i.e., ∂D ∂A ). ∂D ∂A is an N -by-K matrix that can be calculated through the following expression:
where × denotes matrix multiplication and ∂D ∂T is an N -by-T matrix that contains the derivatives of the within-group distance D( X) with respect to the warping functions; the following equation calculates the (n, t)-th element of this matrix:
where sinc is the derivative of the sinc function. In addition, ∂T ∂A denotes the derivatives of the warping functions with respect to the learning parameters which is a T -by-K matrix with the following (t, k)-th element:
Update: We use the Adam algorithm [24] with the step size of 0.01 to update the learning parameters. The Adam algorithm requires ∂D ∂A obtained in the backward step and it performs O(N K) arithmetic operations to update the parameters, where N K is the number of learning parameters.
Computational Complexity of TTW -Calculating Equation (12) is the most computationally expensive step in the TTW algorithm. Time and space complexities needed to calculate this step is O(KN T ). We run this step in a loop for I iterations (I = 100 in our experiments); therefore, time and space complexities of TTW are O(IKN T ) and O(KN T ) respectively. According to our experiments, a good choice for K is a number between 8 and 16 for the datasets of the UCR archive.
We use the UCR (University of California, Riverside) time-series classification archive to conduct experiments of this paper [19]. This archive contains 85 datasets from a wide variety of applications (e.g., Table 1: Results of the DTW averaging experiment. AVG refers to the conventional sample-by-sample averaging technique. We report two numbers in each cell (α-β)%: α is the percentage of the datasets on which the TTW system is significantly better; β is the percentage that the baseline is significantly better. K denotes the number of DST coefficients in Equation (11). Optimal K refers to the experiment in which we validate K for each dataset. TTW K = 4 K = 8 K = 16 K = 32 Optimal K AVG (65-13)% (66-15)% (72-19)% (62-26)% (82-4)% GTW (49-31)% (53-31)% (53-40)% (38-55)% (57-25)% medical imaging, geology and astronomy). Each dataset is split into train and test sets and both sets contain class labels up to 60 classes. All time-series in a dataset of UCR have the same length; for different datasets this length varies from 24 to 2,709 samples.
Baseline System: We implement the GTW algorithm, proposed by Zhou et al. [17], as the baseline method for comparison. The basis functions and other parameters of the implemented GTW are consistent with Section 5.4 of their paper [17].
The goal of DTW averaging task is to find a time-series that have the minimum DTW distance with the given signals [27,28]. In this section, we compare three time-series averaging techniques (i.e., TTW, GTW and AVG) for the DTW averaging task. AVG is the straightforward sample-by-sample averaging technique. We selected GTW and AVG as the baseline systems because their computational complexity is similar or better than TTW.
For class of a dataset, we randomly select ten sets of timeseries, each set contains ten signals. We estimate the DTW average signal of each set using TTW, GTW and AVG methods. We calculate the DTW distance (Equation (1)) for each generated average signal. We finally compare our proposed method with baseline systems (i.e., GTW and AVG) in each dataset by applying pair-wise t-test on the DTW distances. For each dataset, we identify if the systems are comparable (two-sided p-value > 0.05) or one of them (the one with lower mean DTW distance) is significantly better.
Table 1 shows the percentage of the datasets on which TTW vs. the other baseline system is significantly better. In this table, we report the results of TTW with different numbers of DST coefficients (i.e., K). According to Table 1, TTW outperforms both GTW and AVG for all values of K except 32. GTW approximates the temporal warping by combining a set of basis functions. The results show that the GTW algorithm is not flexible enough to find optimal warping functions in many datasets; however, TTW offers a more flexible time warping that is able to learn complex temporal mappings.
The performance of the TTW improves by increasing K up to 8; however, we observe diminishing gains in performance for K > 16.
It is because TTW with a smaller K, estimates a smoother warping function that is able to smooth out local optima and provide a better optimization landscape. Figure 1(a) shows an example of aligning ten time-series using GTW, TTW(K=4) and TTW(K=16). In this example TTW(K=16) converges to a local minimum and cannot find its best solution; however, TTW(K=4) provides a better alignment and considerably outperforms the other methods.
It is important to run TTW with a reasonable value of K: TTW with a small K is not flexible enough to estimate complex warping functions and TTW with a large K is more likely to converge to a local minimum. A straightforward method to define K is to tune it on a subset of the dataset. We use the training sets to tune the value of K and we compare different methods on the test sets. We increase K exponentially from 2 0 to 2 4 and select the one that minimizes the DTW distance of the training set of each dataset. We found that by tuning the value of K for each dataset, our system outperforms GTW on 67% of the datasets (significantly for 57% of the datasets). Table 1 (Optimal K column) summarizes the results of this experiment; the results confirm the importance of tuning K for each dataset.
In this section we investigate the performance of the TTW algorithm in a nearest centroid classification application [29]. We divide the training part of each UCR dataset into two equal parts randomly. We use the first part to find the centroid time-series using TTW and GTW methods. We validate the parameter K on the second part. We select K from five log scale values between 2 0 and 2 4 . We evaluate our centroids using the UCR test partition. In the test phase, conventional DTW distance assigns each sample to a centroid.
Figure 1(b) summarizes the results of this experiment. There are 85 points in this figure; each point compares the classification accuracy of TTW and GTW on a specific dataset. The green rectangular points represent a database for which our proposed algorithm outperforms the GTW algorithm. The results show that the TTW algorithm is able to improve the GTW classification in 61.2% of the datasets.
This paper offers a new solution to the problem of aligning multiple time-series. We first introduce a convolutional sinc kernel that is able to apply non-linear temporal warping functions to time-series. We then use this kernel along with a gradient-based optimization technique to develop a new multiple sequence alignment algorithm: trainable time warping (TTW). The computational complexity of TTW is linear with the number and the length of the input timeseries. Our experiments show that TTW provides an effective time alignment; it outperforms GTW (the previous time-warping algorithm that has similar computational complexity) on DTW averaging and nearest centroid classification tasks.
In this paper, we used a sinc filter to approximate the dirac-delta function. Dirac-delta can also be approximated by other functions such as Gaussian. Future work will explore the effect of using different types of kernels.
This work was partially supported by the National Science Foundation (NSF CAREER 1651740), NIMH R34MH100404, the HC Prechter Bipolar Research Program, and Richard Tam Foundation.
Computational complexity refers to both time and space complexities throughout this paper.
For the sake of simplicity in notations, we assume time-series have the same length. Extending TTW to variable-length signals is straightforward. One method is to interpolate the signals to the maximum length of them.