Soft-DTW: a Differentiable Loss Function for Time-Series

Marco Cuturi, Mathieu Blondel

Introduction

The goal of supervised learning is to learn a mapping that links an input to an output objects, using examples of such pairs. This task is noticeably more difficult when the output objects have a structure, i.e. when they are not vectors (Bakir et al. 2007). We study here the case where each output object is a time series, namely a family of observations indexed by time. While it is tempting to treat time as yet another feature, and handle time series of vectors as the concatenation of all these vectors, several practical issues arise when taking this simplistic approach: Time-indexed phenomena can often be stretched in some areas along the time axis (a word uttered in a slightly slower pace than usual) with no impact on their characteristics; varying sampling conditions may mean they have different lengths; time series may not synchronized.

The DTW paradigm. Generative models for time series are usually built having the invariances above in mind: Such properties are typically handled through latent variables and/or Markovian assumptions (Lütkepohl 2005, Part I,§18). A simpler approach, motivated by geometry, lies in the direct definition of a discrepancy between time series that encodes these invariances, such as the Dynamic Time Warping (DTW) score (Sakoe & Chiba 1971; Sakoe & Chiba 1978). DTW computes the best possible alignment between two time series (the optimal alignment itself can also be of interest, see e.g. Garreau et al. 2014) of respective length $n$ and $m$ by computing first the $n\times m$ pairwise distance matrix between these points to solve then a dynamic program (DP) using Bellman’s recursion with a quadratic $(nm)$ cost.

The DTW geometry. Because it encodes efficiently a useful class of invariances, DTW has often been used in a discriminative framework (with a $k$ -NN or SVM classifier) to predict a real or a class label output, and engineered to run faster in that context (Yi et al. 1998). Recent works by Petitjean et al. 2011; Petitjean & Gançarski 2012 have, however, shown that DTW can be used for more innovative tasks, such as time series averaging using the DTW discrepancy (see Schultz & Jain 2017 for a gentle introduction to these ideas). More generally, the idea of synthetising time series centroids can be regarded as a first attempt to output entire time series using DTW as a fitting loss. From a computational perspective, these approaches are, however, hampered by the fact that DTW is not differentiable and unstable when used in an optimization pipeline.

Soft-DTW. In parallel to these developments, several authors have considered smoothed modifications of Bellman’s recursion to define smoothed DP distances (Bahl & Jelinek 1975; Ristad & Yianilos 1998) or kernels (Saigo et al. 2004; Cuturi et al. 2007). When applied to the DTW discrepancy, that regularization results in a soft-DTW score, which considers the soft-minimum of the distribution of all costs spanned by all possible alignments between two time series. Despite considering all alignments and not just the optimal one, soft-DTW can be computed with a minor modification of Bellman’s recursion, in which all $(\min,+)$ operations are replaced with $(+,\times)$ . As a result, both DTW and soft-DTW have quadratic in time & linear in space complexity with respect to the sequences’ lengths. Because soft-DTW can be used with kernel machines, one typically observes an increase in performance when using soft-DTW over DTW (Cuturi 2011) for classification.

Our contributions. We explore in this paper another important benefit of smoothing DTW: unlike the original DTW discrepancy, soft-DTW is differentiable in all of its arguments. We show that the gradients of soft-DTW w.r.t to all of its variables can be computed as a by-product of the computation of the discrepancy itself, with an added quadratic storage cost. We use this fact to propose an alternative approach to the DBA (DTW Barycenter Averaging) clustering algorithm of (Petitjean et al. 2011), and observe that our smoothed approach significantly outperforms known baselines for that task. More generally, we propose to use soft-DTW as a fitting term to compare the output of a machine synthesizing a time series segment with a ground truth observation, in the same way that, for instance, a regularized Wasserstein distance was used to compute barycenters (Cuturi & Doucet 2014), and later to fit discriminators that output histograms (Zhang et al. 2015; Rolet et al. 2016). When paired with a flexible learning architecture such as a neural network, soft-DTW allows for a differentiable end-to-end approach to design predictive and generative models for time series, as illustrated in Figure 1. Source code is available at https://github.com/mblondel/soft-dtw.

Structure. After providing background material, we show in §2 how soft-DTW can be differentiated w.r.t the locations of two time series. We follow in §3 by illustrating how these results can be directly used for tasks that require to output time series: averaging, clustering and prediction of time series. We close this paper with experimental results in §4 that showcase each of these potential applications.

The DTW and soft-DTW loss functions

DP Recursion. Sakoe & Chiba 1978 showed that the Bellman 1952 equation (Bellman 1952) can be used to compute DTW. That recursion, which appears in line 5 of Algorithm 1 (disregarding for now the exponent $\gamma$ ), only involves $(\min,+)$ operations. When considering kernel $k_{\text{GA}}^{\gamma}$ and, instead, its integration over all alignments (see e.g. Lasserre 2009), Cuturi et al. 2007 and the highly related formulation of Saigo et al. 2004 use an old algorithmic appraoch (Bahl & Jelinek 1975) which consists in (i) replacing all costs by their neg-exponential; (ii) replace $(\min,+)$ operations with $(+,\times)$ operations. These two recursions can be in fact unified with the use of a soft-minimum operator, which we present below.

Unified algorithm Both formulas in Eq. (1) can be computed with a single algorithm. That formulation is new to our knowledge. Consider the following generalized $\min$ operator, with a smoothing parameter $\gamma\geq 0$ :

With that operator, we can define $\gamma$ -soft-DTW:

The original DTW score is recovered by setting $\gamma$ to $0$ . When $\gamma>0$ , we recover $\sdtw=-\gamma\log k_{\text{GA}}^{\gamma}$ . Most importantly, and in either case, $\sdtw$ can be computed using Algorithm 1, which requires $(nm)$ operations and $(nm)$ storage cost as well . That cost can be reduced to $2n$ with a more careful implementation if one only seeks to compute $\sdtw(\mathbf{x},\mathbf{y})$ , but the backward pass we consider next requires the entire matrix $R$ of intermediary alignment costs. Note that, to ensure numerical stability, the operator $\mming$ must be computed using the usual log-sum-exp stabilization trick, namely that $\log\sum_{i}e^{z_{i}}=(\max_{j}z_{j})+\log\sum_{i}e^{z_{i}-\max_{j}z_{j}}.$

2 Differentiation of soft-DTW

A small variation in the input $\mathbf{x}$ causes a small change in $\dtw(\mathbf{x},\mathbf{y})$ or $\sdtw(\mathbf{x},\mathbf{y})$ . When considering $\dtw$ , that change can be efficiently monitored only when the optimal alignment matrix $A^{\star}$ that arises when computing $\dtw(\mathbf{x},\mathbf{y})$ in Eq. (1) is unique. As the minimum over a finite set of linear functions of $\Delta$ , $\dtw$ is therefore locally differentiable w.r.t. the cost matrix $\Delta$ , with gradient $A^{\star}$ , a fact that has been exploited in all algorithms designed to average time series under the DTW metric (Petitjean et al. 2011; Schultz & Jain 2017). To recover the gradient of $\dtw(\mathbf{x},\mathbf{y})$ w.r.t. $\mathbf{x}$ , we only need to apply the chain rule, thanks to the differentiability of the cost function:

With continuous data, $A^{\star}$ is almost always likely to be unique, and therefore the gradient in Eq. (3) will be defined almost everywhere. However, that gradient, when it exists, will be discontinuous around those values $\mathbf{x}$ where a small change in $\mathbf{x}$ causes a change in $A^{\star}$ , which is likely to hamper the performance of gradient descent methods.

An immediate advantage of soft-DTW is that it can be explicitly differentiated, a fact that was also noticed by Saigo et al. 2006 in the related case of edit distances. When $\gamma>0$ , the gradient of Eq. (1) is obtained via the chain rule,

3 Algorithmic differentiation

Differentiating algorithmically $\sdtw(\mathbf{x},\mathbf{y})$ requires doing first a forward pass of Bellman’s equation to store all intermediary computations and recover $R=[r_{i,j}]$ when running Algorithm 1. The value of $\sdtw(\mathbf{x},\mathbf{y})$ —stored in $r_{n,m}$ at the end of the forward recursion—is then impacted by a change in $r_{i,j}$ exclusively through the terms in which $r_{i,j}$ plays a role, namely the triplet of terms $r_{i+1,j},r_{i,j+1},r_{i+1,j+1}$ . A straightforward application of the chain rule then gives

in which we have defined the notation of the main object of interest of the backward recursion: $e_{i,j}\coloneqq\tfrac{\partial r_{n,m}}{\partial r_{i,j}}$ . The Bellman recursion evaluated at $(i+1,j)$ as shown in line 5 of Algorithm 1 (here $\delta_{i+1,j}$ is $\delta(x_{i+1},y_{j})$ ) yields :

which, when differentiated w.r.t $r_{i,j}$ yields the ratio:

The logarithm of that derivative can be conveniently cast using evaluations of $\mming$ computed in the forward loop:

Similarly, the following relationships can also be obtained:

We have therefore obtained a backward recursion to compute the entire matrix $E=[e_{i,j}]$ , starting from $e_{n,m}=\tfrac{\partial r_{n,m}}{\partial r_{n,m}}=1$ down to $e_{1,1}$ . To obtain $\nabla_{\mathbf{x}}\sdtw(\mathbf{x},\mathbf{y})$ , notice that the derivatives w.r.t. the entries of the cost matrix $\Delta$ can be computed by $\tfrac{\partial r_{n,m}}{\partial\delta_{i,j}}=\tfrac{\partial r_{n,m}}{\partial r_{i,j}}\tfrac{\partial r_{i,j}}{\partial\delta_{i,j}}=e_{i,j}\cdot 1=e_{i,j},$ and therefore we have that

Learning with the soft-DTW loss

Note that each $\sdtw(\mathbf{x},\mathbf{y}_{i})$ term is divided by $m_{i}$ , the length of $\mathbf{y}_{i}$ . Indeed, since $\dtw$ is an increasing (roughly linearly) function of each of the input lengths $n$ and $m_{i}$ , we follow the convention of normalizing in practice each discrepancy by $n\times m_{i}$ . Since the length $n$ of $\mathbf{x}$ is here fixed across all evaluations, we do not need to divide the objective of Eq. (9) by $n$ . Averaging under the soft-DTW geometry results in substantially different results than those that can be obtained with the Euclidean geometry (which can only be used in the case where all lengths $n=m_{1}=\dots=m_{N}$ are equal), as can be seen in the intuitive interpolations we obtain between two time series shown in Figure 4.

Non-convexity of $\sdtw$ . A natural question that arises from Eq. (9) is whether that objective is convex or not. The answer is negative, in a way that echoes the non-convexity of the $k$ -means objective as a function of cluster centroids locations. Indeed, for any alignment matrix $A$ of suitable size, each map $\mathbf{x}\mapsto\langle A,\Delta(\mathbf{x},\mathbf{y})\,\rangle$ shares the same convexity/concavity property that $\delta$ may have. However, both $\min$ and $\mming$ can only preserve the concavity of elementary functions (Boyd & Vandenberghe 2004, pp.72-74). Therefore $\sdtw$ will only be concave if $\delta$ is concave, or become instead a (non-convex) (soft) minimum of convex functions if $\delta$ is convex. When $\delta$ is a squared-Euclidean distance, $\dtw$ is a piecewise quadratic function of $\mathbf{x}$ , as is also the case with the $k$ -means energy (see for instance Figure 2 in Schultz & Jain 2017). Since this is the setting we consider here, all of the computations involving barycenters should be taken with a grain of salt, since we have no way of ensuring optimality when approximating Eq. (9).

Smoothing helps optimizing $\sdtw$ . Smoothing can be regarded, however, as a way to “convexify” $\sdtw$ . Indeed, notice that $\sdtw$ converges to the sum of all costs as $\gamma\rightarrow\infty$ . Therefore, if $\delta$ is convex, $\sdtw$ will gradually become convex as $\gamma$ grows. For smaller values of $\gamma$ , one can intuitively foresee that using $\mming$ instead of a minimum will smooth out local minima and therefore provide a better (although slightly different from $\dtw$ ) optimization landscape. We believe this is why our approach recovers better results, even when measured in the original $\dtw$ discrepancy, than subgradient or alternating minimization approaches such as DBA (Petitjean et al. 2011), which can, on the contrary, get more easily stuck in local minima. Evidence for this statement is presented in the experimental section.

2 Clustering with the soft-DTW geometry

The (approximate) computation of $\sdtw$ barycenters can be seen as a first step towards the task of clustering time series under the $\sdtw$ discrepancy. Indeed, one can naturally formulate that problem as that of finding centroids $\mathbf{x}_{1},\dots,\mathbf{x}_{k}$ that minimize the following energy:

To solve that problem one can resort to a direct generalization of Lloyd 1982’s algorithm (Lloyd 1982) in which each centering step and each clustering allocation step is done according to the $\sdtw$ discrepancy.

3 Learning prototypes for time series classification

One of the de-facto baselines for learning to classify time series is the $k$ nearest neighbors ( $k$ -NN) algorithm, combined with DTW as discrepancy measure between time series. However, $k$ -NN has two main drawbacks. First, the time series used for training must be stored, leading to potentially high storage cost. Second, in order to compute predictions on new time series, the DTW discrepancy must be computed with all training time series, leading to high computational cost. Both of these drawbacks can be addressed by the nearest centroid classifier (Hastie et al. 2001, p.670), (Tibshirani et al. 2002). This method chooses the class whose barycenter (centroid) is closest to the time series to classify. Although very simple, this method was shown to be competitive with $k$ -NN, while requiring much lower computational cost at prediction time (Petitjean et al. 2014). Soft-DTW can naturally be used in a nearest centroid classifier, in order to compute the barycenter of each class at train time, and to compute the discrepancy between barycenters and time series, at prediction time.

4 Multistep-ahead prediction

where $\{f_{\theta}\}$ is a set of parameterized function that take as input a time series and outputs a time series. Natural choices would be multi-layer perceptrons or recurrent neural networks (RNN), which have been historically trained with a Euclidean loss (Parlos et al. 2000, Eq.5).

Experimental results

Throughout this section, we use the UCR (University of California, Riverside) time series classification archive (Chen et al. 2015). We use a subset containing 79 datasets encompassing a wide variety of fields (astronomy, geology, medical imaging) and lengths. Datasets include class information (up to 60 classes) for each time series and are split into train and test sets. Due to the large number of datasets in the UCR archive, we choose to report only a summary of our results in the main manuscript. Detailed results are included in the appendices for interested readers.

In this section, we compare the soft-DTW barycenter approach presented in §3.1 to DBA (Petitjean et al. 2011) and a simple batch subgradient method.

Experimental setup. For each dataset, we choose a class at random, pick 10 time series in that class and compute their barycenter. For quantitative results below, we repeat this procedure 10 times and report the averaged results. For each method, we set the maximum number of iterations to 100. To minimize the proposed soft-DTW barycenter objective, Eq. (9), we use L-BFGS.

Qualitative results. We first visualize the barycenters obtained by soft-DTW when $\gamma=1$ and $\gamma=0.01$ , by DBA and by the subgradient method. Figure 5 shows barycenters obtained using random initialization on the ECG200 dataset. More results with both random and Euclidean mean initialization are given in Appendix B and C.

We observe that both DBA or soft-DTW with low smoothing parameter $\gamma$ yield barycenters that are spurious. On the other hand, a descent on the soft-DTW loss with sufficiently high $\gamma$ converges to a reasonable solution. For example, as indicated in Figure 5 with DTW or soft-DTW ( $\gamma=0.01$ ), the small kink around $x=15$ is not representative of any of the time series in the dataset. However, with soft-DTW ( $\gamma=1$ ), the barycenter closely matches the time series. This suggests that DTW or soft-DTW with too low $\gamma$ can get stuck in bad local minima.

When using Euclidean mean initialization (only possible if time series have the same length), DTW or soft-DTW with low $\gamma$ often yield barycenters that better match the shape of the time series. However, they tend to overfit: they absorb the idiosyncrasies of the data. In contrast, soft-DTW is able to learn barycenters that are much smoother.

Quantitative results. Table 1 summarizes the percentage of datasets on which the proposed soft-DTW barycenter achieves lower DTW loss when varying the smoothing parameter $\gamma$ . The actual loss values achieved by different methods are indicated in Appendix G and Appendix H.

As $\gamma$ decreases, soft-DTW achieves a lower DTW loss than other methods on almost all datasets. This confirms our claim that the smoothness of soft-DTW leads to an objective that is better behaved and more amenable to optimization by gradient-descent methods.

2 kk-means clustering experiments

We consider in this section the same computational tools used in §4.1 above, but use them to cluster time series.

Experimental setup. For all datasets, the number of clusters $k$ is equal to the number of classes available in the dataset. Lloyd 1982’s algorithm alternates between a centering step (barycenter computation) and an assignment step. We set the maximum number of outer iterations to $30$ and the maximum number of inner (barycenter) iterations to 100, as before. Again, for soft-DTW, we use L-BFGS.

Qualitative results. Figure 6 shows the clusters obtained when runing Lloyd’s algorithm on the CBF dataset with soft-DTW ( $\gamma=1$ ) and DBA, in the case of random initialization. More results are included in Appendix E. Clearly, DTW absorbs the tiny details in the data, while soft-DTW is able to learn much smoother barycenters.

Quantitative results. Table 2 summarizes the percentage of datasets on which soft-DTW barycenter achieves lower $k$ -means loss under DTW, i.e. Eq. (10) with $\gamma=0$ . The actual loss values achieved by all methods are indicated in Appendix I and Appendix J. The results confirm the same trend as for the barycenter experiments. Namely, as $\gamma$ decreases, soft-DTW is able to achieve lower loss than other methods on a large proportion of the datasets. Note that we have not run experiments with smaller values of $\gamma$ than 0.001, since $\mathbf{dtw}_{0.001}$ is very close to $\dtw$ in practice.

3 Time-series classification experiments

In this section, we investigate whether the smoothing in soft-DTW can act as a useful regularization and improve classification accuracy in the nearest centroid classifier.

Experimental setup. We use 50% of the data for training, 25% for validation and 25% for testing. We choose $\gamma$ from 15 log-spaced values between $10^{-3}$ and $10$ .

Quantitative results. Each point in Figure 7 above the diagonal line represents a dataset for which using soft-DTW for barycenter computation rather than DBA improves the accuracy of the nearest centroid classifier. To summarize, we found that soft-DTW is working better or at least as well as DBA in 75% of the datasets.

4 Multistep-ahead prediction experiments

In this section, we present preliminary experiments for the task of multistep-ahead prediction, described in §3.4.

Experimental setup. We use the training and test sets pre-defined in the UCR archive. In both the training and test sets, we use the first 60% of the time series as input and the remaining 40% as output, ignoring class information. We then use the training set to learn a model that predicts the outputs from inputs and the test set to evaluate results with both Euclidean and DTW losses. In this experiment, we focus on a simple multi-layer perceptron (MLP) with one hidden layer and sigmoid activation. We also experimented with linear models and recurrent neural networks (RNNs) but they did not improve over a simple MLP.

Implementation details. Deep learning frameworks such as Theano, TensorFlow and Chainer allow the user to specify a custom backward pass for their function. Implementing such a backward pass, rather than resorting to automatic differentiation (autodiff), is particularly important in the case of soft-DTW: First, the autodiff in these frameworks is designed for vectorized operations, whereas the dynamic program used by the forward pass of Algorithm 1 is inherently element-wise; Second, as we explained in §2.2, our backward pass is able to re-use log-sum-exp computations from the forward pass, leading to both lower computational cost and better numerical stability. We implemented a custom backward pass in Chainer, which can then be used to plug soft-DTW as a loss function in any network architecture. To estimate the MLP’s parameters, we used Chainer’s implementation of Adam (Kingma & Ba 2014).

Qualitative results. Visualizations of the predictions obtained under Euclidean and soft-DTW losses are given in Figure 1, as well as in Appendix F. We find that for simple one-dimensional time series, an MLP works very well, showing its ability to capture patterns in the training set. Although the predictions under Euclidean and soft-DTW losses often agree with each other, they can sometimes be visibly different. Predictions under soft-DTW loss can confidently predict abrupt and sharp changes since those have a low DTW cost as long as such a sharp change is present, under a small time shift, in the ground truth.

Quantitative results. A comparison summary of our MLP under Euclidean and soft-DTW losses over the UCR archive is given in Table 3. Detailed results are given in the appendix. Unsurprisingly, we achieve lower DTW loss when training with the soft-DTW loss, and lower Euclidean loss when training with the Euclidean loss. Because DTW is robust to several useful invariances, a small error in the soft-DTW sense could be a more judicious choice than an error in an Euclidean sense for many applications.

Conclusion

We propose in this paper to turn the popular DTW discrepancy between time series into a full-fledged loss function between ground truth time series and outputs from a learning machine. We have shown experimentally that, on the existing problem of computing barycenters and clusters for time series data, our computational approach is superior to existing baselines. We have shown promising results on the problem of multistep-ahead time series prediction, which could prove extremely useful in settings where a user’s actual loss function for time series is closer to the robust perspective given by DTW, than to the local parsing of the Euclidean distance.

MC gratefully acknowledges the support of a chaire de l’IDEX Paris Saclay.

References

Appendix A Recursive forward computation of the average path matrix

The average alignment under Gibbs distribution $p_{\gamma}$ can be computed with the following forward recurrence, which mimics closely Bellman’s original recursion. For each $i\in\llbracket n\rrbracket,j\in\llbracket m\rrbracket$ , define

Here terms $r_{ij}$ are computed following the recursion in Algorithm 2. Border matrices are initialized to $0$ , except for $E_{1,1}$ which is initialized to $ $. Upon completion, the average alignment matrix is stored in$ E_{n,m}$.

The operation above consists in summing three matrices of size $(i+1,j+1)$ . There are exactly $(nm)$ such updates. A careful implementation of this algorithm, that would only store two arrays of matrices, as Algorithm 1 only store two arrays of values, can be carried out in $nm\min(n,m)$ space but it would still require $(nm)^{2}$ operations.

Appendix B Barycenters obtained with random initialization

Appendix C Barycenters obtained with Euclidean mean initialization

Appendix D More interpolation results

Left: results obtained under Euclidean loss. Right: results obtained under soft-DTW ( $\gamma=1$ ) loss.