论文复现：soft-DTW：a Differentiable Loss Function for Time-Series

“soft-DTW：a Differentiable Loss Function for Time-Series”原文下载

翻译

摘要：

我们在本文提出一种基于著名的动态时间规整(DTW)差异的可微学习损失函数，用于比较时间序列数据之间的差异。与欧氏距离不同，DTW可以比较具有不同长度的时间序列，并且对于时间维度上的平移或拉伸具有鲁棒性。为了计算DTW,一种常规方法是使用动态规划解决两个时间序列之间的最小成本校准。我们的工作利用了DTW的一种平滑形式，称为soft-DTW，它计算所有校准成本的soft最小值。我们在这篇论文中展示了soft-DTW是一个可微的损失函数，而且它的值和梯度都可以以二次时间/空间复杂度计算（而DTW具有二次时间复杂度但线性空间复杂度）。我们展示了这种正则化特别适用于在DTW几何空间下对时间序列进行平均化和聚类的任务，在这个任务中，我们的方法明显优于现有的基准方法。接下来，我们建议通过以soft-DTW的方式最小化输出时间序列机器与地面真实标签之间的拟合来调整其参数。

1.介绍

监督学习的目标是通过使用输入与输出对象的示例，学习将输入映射到输出的关系。当输出对象具有结构时，例如它们不是向量，这个任务明显更加困难。我们在这里研究每个输出对象都是时间序列的情况，即一系列随时间索引的观测值。虽然将时间视为另一个特征，并将向量时间序列处理为所有这些向量的串联是很诱人的，但采用这种简单的方法会导致几个实际问题：时间索引的现象通常在时间轴上的某些区域被拉伸（例如，以稍慢于平常的速度发出的单词），但其特征不受影响；不同的采样条件可能导致它们具有不同的长度；时间序列可能不同步。

DTW模型 ｜　用于时间序列的生成模型通常在考虑上述不变性的基础上构建：这些属性通常通过潜在变量和/或马尔可夫假设来处理。一种受到几何学启发的更简单的方法，是通过直接定义时间序列之间的差异来实现这些不变性，例如DTW score。DTW计算两个时间序列之间的最佳对齐（最优对齐本身也可能是感兴趣的），这两个时间序列的长度分别为$n$和$m$。首先，通过计算这些点之间的$n\times m$成对距离矩阵，然后使用Bellman递推法进行动态规划（DP），其计算成本为二次复杂度$(nm)$。动态时间规整（DTW）的几何性质。由于DTW有效地编码了一类有用的不变性，它经常被用于判别性框架中（例如，使用k-NN或SVM分类器），以预测实数或类别标签的输出，并在该背景下进行了优化以提高运行速度。

DTW的几何性质 ｜　由于DTW有效地编码了一类有用的不变性，它经常被用于判别性框架中（例如，使用k-NN或SVM分类器），以预测实数或类别标签的输出，并在该背景下进行了优化以提高运行速度。最近Petitjean等人在2011年的研究和Petitjean同Ganc¸arski在2012年的研究表明，DTW可以用于更具创新性的任务，例如使用DTW差异进行时间序列平均化。更普遍地说，将时间序列的中心合成为一个整体可以被视为使用DTW作为拟合损失输出整个时间序列的首次尝试。然而，从计算角度来看，这些方法受到DTW不可微分且在优化过程中不稳定的影响。

Soft-DTW ｜　与这些发展平行的是，一些作者考虑使用Bellman的递归的平滑修改来定义平滑的DP距离或核函数。当应用于DTW差异时，这种正则化导致了soft-DTW得分，它考虑了两个时间序列之间所有可能对齐所涵盖的所有成本的软最小值。尽管考虑了所有对齐而不仅仅是最优对齐，但soft-DTW可以通过对Bellman递归进行轻微修改来计算，在这种修改中，所有的$(\mathbf{min}, +)$运算被替换为$(+,\times)$ 运算。因此，无论是DTW还是soft-DTW，都具有相对于序列长度的时间复杂度为二次的特点，而空间复杂度为线性。由于soft-DTW可以与核机器一起使用，在分类任务中使用soft-DTW相对于DTW通常会观察到性能的提升。

我们的贡献 ｜　在本论文中，我们探索了平滑动态时间规整（DTW）的另一个重要优势：与原始的DTW差异不同，soft-DTW在所有参数上都是可微分的。我们展示了soft-DTW相对于其所有变量的梯度可以作为计算差异本身的副产品得到，但会增加二次存储成本。我们利用这一事实提出了一种替代方法来处理（Petitjean等人，2011年）提出的DTW质心聚类算法（DBA），并观察到我们的平滑方法在该任务中明显优于已知的基准方法。更一般地说，我们建议使用soft-DTW作为拟合项来比较合成时间序列片段的机器输出与地面真实观测的差异，就像例如，正则化的Wasserstein距离用于计算质心（Cuturi＆Doucet，2014年），后来用于拟合输出直方图的判别器（Zhang等人，2015年；Rolet等人，2016年）一样。当与灵活的学习架构（如神经网络）配对时，soft-DTW允许进行可微的端到端方法，用于设计时间序列的预测和生成模型，如图1所示。在我们提出的方法中，soft-DTW作为一个拟合项，可以用于将合成的时间序列片段与地面真实观测进行比较，从而设计预测和生成时间序列的机器学习模型。在使用灵活的学习架构（如神经网络）时，soft-DTW允许实现端到端的可微分学习方法。这种方法可以应用于多种时间序列相关的任务，如分类、聚类等。我们的实验结果表明，使用soft-DTW可以在一些任务中明显提高性能，超过了传统的DTW方法。

本文中我们探讨了DTW平滑的另一个重要优点：与原始的DTW差异相比，soft-DTW在其所有参数上都是可微的。我们展示了soft-DTW相对于其所有变量的梯度可以作为计算差异本身的副产品进行计算，但增加了二次存储成本。我们利用这个事实来提出一种替代方法来处理DBA（DTW质心平均）聚类算法（Petitjean等，2011），并观察到我们的平滑方法在该任务中显著优于已知的基准。更一般地，我们建议使用soft-DTW作为拟合项来比较机器合成时间序列段与地面真实观察结果之间的差异，就像我们在计算质心（Cuturi和Doucet，2014）或拟合输出直方图的鉴别器（Zhang等，2015；Rolet等，2016）时使用正则化的Wasserstein距离一样。当与灵活的学习架构（如神经网络）配对使用时，soft-DTW允许采用可微的端到端方法来设计时间序列的预测和生成模型，如图1所示。源代码可在 https://github.com/mblondel/soft-dtw 上找到。

结构｜　在提供背景材料后，我们在第2节中展示了如何针对两个时间序列的位置对soft-DTW进行微分。接着在第3节中，我们说明了这些结果如何直接用于需要输出时间序列的任务，如平均、聚类和时间序列的预测。在第4节中，我们通过实验结果展示了每个潜在应用的效果，从而结束了本文。

符号｜　接下来，我们考虑取值范围在$\Omega \subset \mathbb{R}^p$的多变量离散时间序列，这些时间序列的长度可以不同。因此，一个时间序列可以表示为一个包含$p$行和可变列数的矩阵。我们考虑一个可微的替代成本函数$\delta :\mathbb{R}^p \times \mathbb{R}^p \rightarrow \mathbb{R}+$，在大多数情况下，这个函数将是两个向量之间的平方欧氏距离。对于一个整数$n$，我们用$[n]$表示整数集${1,…,n}$。给定两个序列的长度$n$和$m$，我们用$\mathcal{A}{n,m} \subset {0,1}^{n\times m}$表示(二进制)对齐矩阵的集合，即在一个$n\times m$上，通过使用$\downarrow,\rightarrow,\And$移动将左上角$(1, 1)$矩阵元素连接到右下角$(n, m)$矩阵元素的路径。这些对齐矩阵是由$0$和$1$组成的二进制矩阵，表示序列之间的匹配关系。集合$\mathcal{A}_{n,m}$的基数被称为 delannoy$(n-1,m-1)$，这个数随着$m$和$n$的增长呈指数级增加。

2.DTW和soft-DTW损失函数

在本节中，我们提出了原始的DTW差异（Sakoe＆Chiba，1978）和全局对齐核（Global Alignment kernel，简称GAK）（Cuturi等，2007）的统一形式。这两者都可用于比较两个时间序列 $x = (x1, . . . , xn) ∈ \mathbb{R}^{p*n}$ 和 $y = (y1, . . . , ym) ∈ \mathbb{R}^{p*m}$。

2.1.对齐成本：最优性与求和

给定成本矩阵$\Delta(x,y) := [\delta(x_i,y_j)]{ij} \in \mathbb{R}^{n\times m}$，该矩阵与一个在$\mathcal{A}{n,m}$中的对齐矩阵$A$的内积$\langle A,\Delta(x,y) \rangle$给出了$A$的分数，如图2所示。无论是 DTW 还是 GAK 都考虑了所有可能的对齐矩阵的成本，但它们的计算方式不同：

$$ \mathbf{DTW}(x,y) := \mathbf{min}{A \in \mathcal{A}{n,m}} \langle A,\Delta(x,y) \rangle , $$

$$ k_{GA}^\gamma(\mathbb{x},\mathbb{y}) := \sum_{A \in \mathcal{A}_{n,m}} e^{-\langle A,\Delta(\mathbb{x},\mathbb{y}) \rangle / \gamma} \tag{1} $$

DP递归 ｜　Sakoe和Chiba在1978年的论文中证明了Bellman方程可以用于计算DTW。这个递归在算法1的第5行中出现（暂时忽略指数$\gamma$），仅涉及$(\mathbf{min}, +)$运算。当考虑核函数$k^{\gamma}_{GA}$ 及其在所有对齐中的积分时（详见 Lasserre 2009），Cuturi 等人（2007 年，定理 2）以及与之高度相关的 Saigo 等人的公式（2004 年，第 1685 页）采用了一种旧的算法方法（Bahl 和 Jelinek，1975 年），其中包括$(i)$将所有成本替换为其负指数；$(ii)$将$(\mathbf{min}, +)$运算替换为$(+,\times)$运算。事实上，这两个递归可以通过使用soft最小值运算符进行统一，下面我们将介绍这个方法。

统一算法 ｜　这两个公式在公式（1）中可以使用一个单一的算法来计算。据我们所知，这种表述是新的。考虑以下带有平滑参数$\gamma \geq 0$的广义$\mathbf{min}$运算符：

$$ \mathbb{min^\gamma}{a_1,…,a_n} := \begin{cases} \mathbb{min_{i \leq n}}a_i, & \gamma=0,\ -\gamma \mathbb{log} \sum_{i=1}^n e^{-a_i/{\gamma}}, & \gamma < 0. \end{cases} $$ 使用上面提到的泛化最小值运算符，我们可以定义$\gamma\mathbb{-soft-}\mathbf{DTW}$： $$ \boxed{\mathbf{dtw}\gamma(\mathbb{x},\mathbb{y}) := \mathbf{min}^{\gamma} {\langle A,\Delta(\mathbb{x},\mathbb{y}) \rangle ,A \in \mathcal{A}{n,m} }.} $$

将$\gamma$设置为0，可以恢复原始的DTW分数。当$\gamma > 0$时，我们可以恢复$\mathbb{dtw_{\gamma}} = -\gamma \mathbb{log}k_{GA}^\gamma$。最重要的是，在任何情况下，都可以使用算法1来计算$\mathbb{dtw_{\gamma}}$，该算法需要$(nm)$的操作和$(nm)$的存储成本。如果只想计算$\mathbb{dtw_{\gamma}}(\mathbb{x},\mathbb{y})$，则可以通过更仔细的实现将这个成本减少到$2n$，但我们接下来考虑的反向传播需要中间对齐成本矩阵$R$的整个矩阵。请注意，为了确保数值稳定性，$\mathbb{min^{\gamma}}$操作符必须使用常规的对数-求和-指数稳定化\技巧进行计算，即：$\mathbb{log}\sum_ie^{z_i} = (\mathbb{max_j}z_j) + \mathbb{log}\sum_i e^{z_j-\mathbb{max_j}z_j}.$

2.2.对 soft-DTW 求导

输入$\mathbb{x}$的微小变化会导致$\mathbb{dtw_{0}}(\mathbb{x},\mathbb{y})$或$\mathbb{dtw_{\gamma}}(\mathbb{x},\mathbb{y})$的微小变化。当考虑$\mathbb{dtw_{0}}$时，只有在计算方程1的$\mathbb{dtw_{0}}(\mathbb{x},\mathbb{y})$时产生的最优对齐矩阵$A^$是唯一的情况下，这种变化才能高效地被监测到。由于是$\Delta$的有限线性函数的最小值，因此$\mathbb{dtw_{0}}$对于成本矩阵$\Delta$在局部是可微的，其梯度为$A^$。这一事实已经被用于所有基于 DTW 距离的时间序列平均算法（Petitjean et al., 2011；Schultz & Jain, 2017）。为了计算$\mathbb{dtw_{0}}(\mathbb{x},\mathbb{y})$对于$\mathbb{x}$的梯度，我们只需要应用链式法则，得益于成本函数的可微性： $$ \nabla_{\mathbb{x}} \mathbb{dtw_{0}}(\mathbb{x},\mathbb{y}) = \Bigl(\frac{\partial \Delta(\mathbb{x},\mathbb{y})}{\partial\mathbb{x}}\Bigl)^T A^*,\tag2 $$

其中，${\partial \nabla(\mathbb{x},\mathbb{y})}/{\partial{x}}$是关于$\mathbb{x}$的$\Delta$的雅可比矩阵，是从$\mathbb{R}^{p \times n}$到$\mathbb{R}^{n \times m}$ 的线性映射。当$\delta$是平方欧几里得距离时，该雅可比矩阵对矩阵$B \in \mathbb{R}^{n \times m}$进行转置（$\circ$表示逐元素乘积）： $$ ({\partial{\Delta(\mathbb{x},\mathbb{y})}}/{\partial{\mathbb{x}}})^T B = 2((\boldsymbol{1}_p \boldsymbol{1}_m^T B^T)\circ{\mathbb{x}} - \mathbb{y}B^T) $$

在连续数据的情况下，$A^$几乎总是唯一的，因此方程（2）中的梯度几乎在任何地方都是有定义的。然而，当存在梯度时，在那些小的$\mathbb{x}$变化引起$A^$变化的值$\mathbb{x}$周围，梯度将是不连续的，这可能会影响梯度下降方法的性能。

$\gamma > 0$的情况 ｜　soft-DTW 的一个直接优势是它可以被明确地求导，这也是Saigo 等人在编辑距离相关案例中注意到的事实（2006年）。当$\gamma > 0$时，通过链式法则可以得到方程（1）的梯度：

$$ \nabla_{\mathbb{x}} \mathbb{dtw_{\gamma}}(\mathbb{x},\mathbb{y}) = \Bigl(\frac{\partial \Delta(\mathbb{x},\mathbb{y})}{\partial\mathbb{x}}\Bigl)^T \mathbb{E}{\gamma}[A],\tag3 $$ $$ \mathbb{where}\quad \mathbb{E}{\gamma}[A] := \frac{1}{k_{GA}^\gamma(\mathbb{x},\mathbb{y})} \quad \sum_{A \in \mathcal{A}{n,m}}e^{-\langle A,{\Delta(\mathbb{x},\mathbb{y})/{\gamma}} \rangle}A, $$ 这里的$A^*$是在 Gibbs 分布$p{\gamma} \propto e^{-\langle A,{\Delta(\mathbb{x},\mathbb{y})}\rangle /{\gamma}}$下的平均对齐矩阵$A$，该分布在所有的对齐矩阵$\mathcal{A}{n,m}$上进行定义。因此，核函数$k{GA}^\gamma(\mathbb{x},\mathbb{y})$可以被解释为$p_{\gamma}$的归一化常数。当然，由于$\mathcal{A}{n,m}$在$n$和$m$中呈指数级增长，简单的求和是不可行的。虽然存在一个用于计算该平均对齐矩阵$\mathbb{E}{\gamma}[A]$的 Bellman 递归（见附录 A），但这个计算具有四次方的复杂性$n^2m^2$。

需要注意的是，这与 Saigo 等人（2006）对编辑距离得到的二次复杂性形成了鲜明的对比，这是因为他们考虑的序列只能取有限字母表中的值。为了计算 soft-DTW 的梯度，我们提出了一种算法，其复杂性在计算上仍然是二次的$(nm)$。实现这种降低的关键是按照算法1中给出的 Bellman 递归的逆序应用链式规则，即反向传播。最近在（Blondel 等人，2016）中也使用了类似的思想来计算 ANOVA 核的梯度。

2.3.算法微分

通过算法微分来区分地计算$\mathbb{dtw_{\gamma}}(\mathbb{x},\mathbb{y})$需要首先进行Bellman方程的前向计算，以存储所有的中间计算结果并在运行算法1时恢复出$R = [r_{i,j}]$。在前向递归的末尾，存储在$r_{n,m}$中的$\mathbb{dtw_{\gamma}}(\mathbb{x},\mathbb{y})$的值会受到$r_{i,j}$的变化的影响，这主要是通过$r_{i,j}$在其中起到作用的三个项来实现的，即$r_{i+1,j},r_{i,j+1},r_{i+1,j+1}$。然后，直接应用链式法则可以得到如下结果： $$ ….. $$ 其中我们定义了后向递归的主要对象的符号表示：$e_{i,j} := \frac{\partial r_{n,m}}{\partial r_{i,j}}.$在算法1的第5行所示的$(i+1,j)$处计算的Bellman递归（这里$\delta_{i+1,j}$表示$\delta(x_{i+1},y_j)$得到如下结果： $$ r_{i+1,j} = \delta_{i+1,j} + \mathbb{min}^{\gamma}{r_{i,j-1},r_{i,j},r_{i+1,j-1} }, $$

在对$r_{i,j}$微分时得到的比率为（w.r.t是 “with respect to” 的缩写，表示针对某个变量或者参数进行微分或者考虑。在数学和科学领域中常常使用这个缩写来指示微分或者导数是关于某个特定的变量）： $$ \frac{\partial r_{i+1,j}}{\partial r_{i,j}} = e^{-r_{i,j}/\gamma}/(e^{-r_{i,j-1}/\gamma} + e^{-r_{i,j}/\gamma} + e^{-r_{i+1,j-1}/\gamma}). $$

该导数的对数可以通过在前向循环中计算的$\mathbb{min}{\gamma}$的评估来方便地表示为： $$ \gamma \mathbb{log} \frac{\partial r{i+1,j}}{\partial r_{i,j}} = \mathbb{min}^{\gamma}{r_{i,j-1},r_{i,j},r_{i+1,j-1} } - r_{i,j} \=r_{i+1,j} - \delta_{i+1,j} - r_{i,j}. $$

类似地，还可以得到以下关系： $$ \gamma \mathbb{log} \frac{\partial r_{i,j+1}}{\partial r_{i,j}} =r_{i,j+1} - \delta_{i,j+1} - r_{i,j}, \\gamma \mathbb{log} \frac{\partial r_{i+1,j+1}}{\partial r_{i,j}} =r_{i+1,j+1} - \delta_{i+1,j+1} - r_{i,j}. $$

因此，我们已经得到了一个向后递归来计算整个矩阵$E = [e_{i,j}]$，从$e_{n,m} = \frac{\partial r_{n,m}}{\partial r_{n,m}} = 1$开始，一直到$e_{1,1}$。为了获得$\nabla_{\mathbb{x}}\mathbb{dtw_{\gamma}}(\mathbb{x},\mathbb{y})$，注意到关于代价矩阵$\Delta$的各个条目的导数可以通过$\frac{\partial r_{n,m}}{\partial \delta_{i,j}} = \frac{\partial r_{n,m}}{\partial r_{i,j}} \frac{\partial r_{i,j}}{\partial \delta_{i,j}} = e_{i,j} \cdot 1 = e_{i,j}$来计算，因此我们有：

$$ \nabla_{\mathbb{x}}\mathbb{dtw_{\gamma}}(\mathbb{x},\mathbb{y}) = \Big(\frac{\partial \Delta(\mathbb{x},\mathbb{y})}{\partial \mathbb{x}}\Big)^T E, $$ 其中，$E$正好是方程（3）中的平均对齐$\mathbb{E}{\gamma}[A]$。这些计算总结在算法2中，一旦$\Delta$被计算出来，其复杂度在时间和空间上为$nm$。由于$\mathbb{min}^{\gamma}$具有$1/\gamma-$Lipschitz 连续梯度，当 $\delta$是平方欧氏距离时，$\mathbb{dtw}{\gamma}$的梯度是$2/\gamma-$Lipschitz 连续的。

3.学习使用 soft-DTW 损失

3.1.使用soft-DTW几何方法进行平均化

在本节中，我们将算法2直接应用于计算相对于$\mathbf{dtw}\gamma$偏差的时间序列的Frechet均值(1948)。给定一组$N$个时间序列$y_1,…,y_N$，即$N$个矩阵，每个矩阵包含$p$行和可变列数$m_1,…,m_N$，我们希望在一组标准化权重 $\lambda_1,...,\lambda_N \in R^+$下定义一个单一的重心时间序列$\mathbf{x}$，满足$\sum{i=1}^N{\lambda_i = 1}$。我们的目标是近似解决以下问题，其中我们假设$\mathbf{x}$的长度为固定值$n$： $$ \mathbf{min}{\mathbf{x} \in \mathbb{R}^{p\times n}} \sum^N{i=1}{\frac{\lambda_i}{m_i}\mathbf{dtw}_\gamma(\mathbf{x},\mathbf{y}i)}.\tag{4} $$ 请注意，每个$\mathbf{dtw}\gamma(\mathbf{x},\mathbf{y}_i)$项都除以$m_i$，即$\mathbf{y}_i$的长度。实际上，由于$\mathbf{dtw}_0$函数与输入长度$n$和$m_i$呈递增（大致线性）关系，我们遵循将每个偏差通过$n\times m_i$进行标准化的惯例。由于$\mathbf{x}$的长度$n$在所有评估中是固定的，因此我们不需要将公式（4）的目标除以$n$。在 soft-DTW 几何中进行平均会得到与欧氏几何截然不同的结果（欧氏几何只适用于$n = m_1 = … = m_N$所有长度相等的情况），如图4所示，我们可以看到在两个时间序列之间得到的直观插值。

$\mathbf{dtw}_\gamma$的非凸性 ｜　从公式(4)中可以自然地引出一个问题，即该目标函数是否是凸函数。答案是否定的，这种情况类似于$k$-means目标函数作为群集质心位置的非凸性。确实，对于合适大小的任何对齐矩阵$A$，每个映射$\mathbb{x} \mapsto {\langle A,\Delta(\mathbb{x},\mathbb{y}) \rangle}$具有与$\delta$相同的凸凹特性。然而，$\mathbb{min}$和$\mathbb{min^\gamma}$只能保持基本函数的凹性。因此，如果$\delta$是凹函数，$\mathbb{dtw_\gamma}$将只是凹的，但是如果$\delta$是凸函数，则会成为凸函数的（非凸）（软）最小值。当$\delta$是平方欧氏距离时，$\mathbb{dtw_0}$是$\mathbb{x}$的分段二次函数，这也与 $k$-means 能量的情况相同（例如，参见 Schultz & Jain 2017 中的图2）。由于这是我们在这里考虑的设置，所有涉及重心计算的过程应该谨慎对待，因为在近似等式 (4) 时我们无法确保最优解。

平滑有助于优化$\mathbf{dtw}\gamma$ ｜　然而，平滑可以被视为"凸化"$\mathbf{dtw}\gamma$的一种方法。事实上，注意到当$\gamma \rightarrow \infty$，$\mathbf{dtw}\gamma$收敛到所有成本的总和。因此，如果$\delta$是凸的，随着$\gamma$的增长，$\mathbf{dtw}\gamma$将逐渐变得凸。对于较小的$\gamma$，可以直观地预见使用$\mathbb{min}_{\gamma}$而不是最小值将平滑局部极小值，因此提供更好（与$\mathbf{dtw}_0$稍有不同）的优化环境。我们相信这就是为什么我们的方法比次梯度或交替最小化方法（例如DBA，Petitjean等人，2011年）在原始的$\mathbf{dtw}_0$差异度量下表现更好的原因。后者相反，更容易陷入局部最小值。对于这一说法，我们在实验部分中提供了证据。

3.2.使用soft-DTW几何方法进行聚类

$\mathbf{dtw}_\gamma$重心的(近似)计算可以看作是在$\mathbf{dtw}_\gamma$差异度量下对时间序列进行聚类的第一步。事实上，可以自然地将该问题表述为找到最小化以下能量的质心$\mathbb{x}1,…,\mathbb{x}k$的问题： $$ \mathbf{min}{\mathbb{x}1,…,\mathbb{x}k \in \mathbb{R}^{p\times n}} \sum^N{i=1}{\frac{1}{m_i} \mathbf{min}{j \in [[k]]} \mathbf{dtw}\gamma(\mathbf{x}_j,\mathbf{y}i)}.\tag{5} $$ 为了解决这个问题，我们可以采用Lloyd算法（1982）的直接推广，其中每个中心步骤和每个聚类分配步骤都是根据$\mathbf{dtw}\gamma$差异度量进行。

3.3.学习时间序列分类中的原型

学习时间序列分类的实际基准之一是$k$-NN算法，结合DTW作为时间序列之间的差异度量。然而，$k$-NN有两个主要缺点。首先，用于训练的时间序列必须被存储，导致潜在的高存储成本。其次，在对新的时间序列进行预测时，必须使用所有训练时间序列计算DTW差异度量，导致计算成本高昂。这两个缺点都可以通过最近质心分类器来解决（Hastie et al.，2001，p.670）（Tibshirani et al.，2002）。这种方法选择距离待分类时间序列最近的类别的质心（中心点）。尽管非常简单，但已被证明在预测时需要更低的计算成本，与$k$-NN方法相比具有竞争力（Petitjean et al.，2014）。在最近质心分类器中，可以自然地使用soft-DTW来在训练时计算每个类别的质心，并在预测时计算质心与时间序列之间的差异。

3.4.多步领先预测

Soft-DTW非常适合用作需要时间序列输出的任何任务的损失函数。作为这类任务的一个示例，我们考虑以下问题：给定时间序列的前$1,…,t$个观测值，预测剩余的$(t+1),…,n$个观测值。令$\mathbb{x}^{t,t’} \in \mathbb{R}^{p \times (t’-t+1)}$之间的所有列的子矩阵，其中。学习预测时间序列的片段可以被视为以下问题： $$ \mathbb{min}{\theta \in \Theta} \sum{i=1}^N \mathbb{dtw}\gamma(f\theta(\mathbb{x}_i^{1,t}),\mathbb{x}i^{t+1,n}), $$ 其中${ f\theta }$是一组参数化函数，它们以时间序列作为输入并输出时间序列。自然的选择可能是多层感知器或循环神经网络（RNN），它们在历史上使用欧氏损失进行训练（Parlos et al., 2000, Eq.5）。

4.实验结果

在本节中，我们使用了加州大学河滨分校（UCR）时间序列分类数据库（Chen等，2015）。我们使用了一个包含79个数据集的子集，涵盖了广泛的领域（天文学、地质学、医学影像学）和不同长度的时间序列。这些数据集包含了每个时间序列的类别信息（最多60个类别），并被分为训练集和测试集。由于UCR数据库中数据集的数量较大，我们选择在主要文稿中仅报告结果摘要，详细结果可以在附录中找到，供有兴趣的读者查阅。

4.1.平均实验

在这一节中，我们将3.1中介绍的soft-DTW质心方法与DBA（Petitjean等，2011）和简单的批次次梯度方法进行比较。

实验设置 ｜　对于每个数据集，我们随机选择一个类别，在该类别中选择10个时间序列，并计算它们的重心。对于下面的定量结果，我们重复这个过程10次，并报告平均结果。对于每种方法，我们将最大迭代次数设置为100次。为了最小化提出的soft-DTW重心目标，即公式(4)，我们使用L-BFGS算法。

定性结果 ｜　我们首先通过soft-DTW在$\gamma = 1$和$\gamma = 0.01$时，通过DBA和次梯度方法获得的重心进行可视化。图5展示了在ECG200数据集上使用随机初始化得到的重心。在附录B和C中提供了更多随机和欧几里德均值初始化的结果。

我们观察到，无论是DBA还是具有较低平滑参数$\gamma$的soft-DTW，都会产生虚假的重心。另一方面，当$\gamma$足够高时，对spft-DTW损失进行下降会收敛到一个合理的解。例如，如图5所示，使用DTW或soft-DTW$(\gamma = 0.01)$，在$x = 15$附近的小弯曲不代表数据集中任何时间序列的特征。然而，使用soft-DTW$(\gamma = 1)$，重心与时间序列紧密匹配。这表明具有过低$\gamma$的DTW或soft-DTW可能会陷入糟糕的局部最小值中。

当使用欧几里得平均初始化（仅在时间序列具有相同长度时可能）时，具有较低$\gamma$的DTW或soft-DTW往往会产生更好地匹配时间序列形状的重心。然而，它们倾向于过拟合：它们会吸收数据的特异性。相比之下，soft-DTW能够学习到更平滑的重心。

定量结果 ｜　表1总结了在不同的平滑参数$\gamma$下，提出的soft-DTW重心方法在百分之多少的数据集上实现了较低的DTW损失。不同方法实际达到的损失值在附录G和附录H中有所指示。

随着$\gamma$的减小，soft-DTW在几乎所有数据集上都实现了比其他方法更低的DTW损失。这证实了我们的观点，即soft-DTW的平滑性使得目标函数更加良好，更容易通过梯度下降方法进行优化。

4.2.$k$均值聚类实验

在本节中，我们使用与上面第4.1节中相同的计算工具，但将它们用于对时间序列进行聚类。

实验设置 ｜　对于所有数据集，簇的数量$k$等于数据集中可用的类别数量。Lloyd算法在中心步骤（质心计算）和分配步骤之间交替进行。我们将外部迭代的最大次数设置为30次，将内部（质心）迭代的最大次数设置为100次，与之前一样。同样地，对于soft-DTW，我们使用L-BFGS方法。

定性结果 ｜　图6展示了在CBF数据集上运行Lloyd算法时，使用soft-DTW$(\gamma = 1)$和DBA进行随机初始化时获得的聚类结果。更多结果在附录E中。显然，DTW吸收了数据中的微小细节，而soft-DTW能够学习到更平滑的质心。

定量结果 ｜　表2总结了在哪些数据集上，soft-DTW质心在DTW下（即$\gamma = 0$时的方程（5））能够实现更低的$k$-means损失。所有方法实际实现的损失值在附录I和附录J中有所说明。结果证实了与质心实验相同的趋势。也就是说，随着$\gamma$的减小，soft-DTW能够在大部分数据集上实现比其他方法更低的损失。请注意，我们没有对$\gamma$小于0.001的更小值进行实验，因为实际上$\mathbb{dtw}{0.001}$非常接近$\mathbb{dtw}{0}$。

4.3.时间序列分类实验

在这一部分，我们研究了soft-DTW中的平滑是否能够作为一种有效的正则化方法，提高最近质心分类器的分类准确率。

实验设置 ｜　我们将50%的数据用于训练，25%用于验证，另外25%用于测试。我们从$10^{-3}$到$10$之间取15个对数间隔的$\gamma$值。

定量结果 ｜　在上方对角线上的每个点代表了一个数据集，对于这些数据集，使用soft-DTW进行质心计算而不是DBA可以提高最近质心分类器的准确性。总结起来，我们发现在75%的数据集中，soft-DTW的效果要么更好，要么至少与DBA相当。

4.4.多步领先预测实验

在这一部分，我们探究了soft-DTW中的平滑是否能够作为一种有效的方法，来提升多步预测任务的性能。

实验设置 ｜　我们使用预先定义的UCR存档中的训练集和测试集。在训练集和测试集中，我们将时间序列的前60%作为输入，剩余的40%作为输出，忽略类别信息。然后我们使用训练集来学习一个模型，用于从输入预测输出，并使用测试集来评估使用欧氏距离和DTW距离的结果。在这个实验中，我们专注于一个简单的多层感知器（MLP），包含一个隐藏层和Sigmoid激活函数。我们还尝试了线性模型和循环神经网络（RNN），但它们并没有在简单的MLP上取得更好的效果。

实验细节 ｜　深度学习框架，如Theano、TensorFlow和Chainer，允许用户为自己的函数指定自定义的反向传播。在soft-DTW的情况下，实现这样的反向传播，而不是采用自动微分（autodiff），尤其重要：首先，这些框架中的自动微分是为向量化操作而设计的，而正向传播中使用的动态规划本质上是逐元素的；其次，正如我们在第2.2节中解释的那样，我们的反向传播能够重新使用正向传播中的log-sum-exp计算，从而降低计算成本并提高数值稳定性。我们在Chainer中实现了自定义的反向传播，然后可以将soft-DTW作为损失函数插入任何网络架构中。为了估计MLP的参数，我们使用了Chainer对Adam（Kingma＆Ba，2014）的实现。

定性结果 ｜　我们在图1中展示了在欧氏距离和soft-DTW距离下获得的预测结果的可视化，附录F中也有类似的结果。我们发现对于简单的一维时间序列，MLP效果非常好，显示出其捕捉训练集中模式的能力。虽然在欧氏距离和soft-DTW距离下的预测通常是一致的，但有时它们可能在视觉上有所不同。在soft-DTW距离下的预测可以自信地预测突然和尖锐的变化，因为只要这种尖锐的变化在地面实况中以一个小的时间偏移存在，它的DTW成本就会很低。

定量结果 ｜　我们在表3中提供了我们在欧氏距离和soft-DTW距离下对UCR存档中的MLP的比较摘要。详细的结果在附录中给出。毫不奇怪，当使用soft-DTW损失进行训练时，我们实现了较低的DTW损失，当使用欧氏距离损失进行训练时，我们实现了较低的欧氏距离损失。由于DTW对于一些有用的不变性是鲁棒的，用soft-DTW意义上的一个小错误可能在许多应用中比欧氏距离意义上的错误更为明智的选择。

5.结论

本文提出将流行的时间序列之间的DTW差异转化为一个完整的损失函数，用于衡量基准实际时间序列与学习机输出之间的差异。我们通过实验表明，在现有的计算时间序列数据的贝叶斯中心和簇问题时，我们的计算方法优于现有的基准方法。在多步预测时间序列问题上，我们展示了有希望的结果，这在用户的实际时间序列损失函数更接近于DTW所提供的鲁棒性而不是欧氏距离的局部解析时，可能会非常有用。

致谢｜　非常感谢IDEX Paris Saclay的支持。