综合转载于以下文章：
动态时间规整/规划（Dynamic Time Warping，DTW）是一个比较老的算法，大概在1970年左右被提出来，最早用于处理语音方面识别分类的问题。
简介
简单来说，给定两个离散的序列（实际上不一定要与时间有关），DTW能够衡量这两个序列的相似程度，或者说两个序列的距离。同时DTW能够对两个序列的延展或者压缩能够有一定的适应性，举个例子，不同人对同一个词语的发音会有细微的差别，特别在时长上，有些人的发音会比标准的发音或长或短，DTW对这种序列的延展和压缩不敏感，所以给定标准语音库，DTW能够很好得识别单个字词，这也是为什么DTW一直被认为是语音处理方面的专门算法。实际上，DTW虽然老，但简单且灵活地实现模板匹配，能解决很多离散时间序列匹配的问题，视频动作识别，生物信息比对等等诸多领域都有应用。
例如下图，有两个呈现正弦规律序列，其中蓝色序列是稍微被拉长了。即使这两个序列，不重合，但是我们也可以有把握说这两个序列的相似程度很高（或者说这两个序列的距离很小）。
DTW能够计算这两个序列的相似程度，并且给出一个能最大程度降低两个序列距离的点到点的匹配。见下图，其中黑色与红色曲线中的虚线就是表示点点之间的一个对应关系。
也就是说，两个比对序列之间的特征是相似的，只是在时间上有不对齐的可能，这个算法名中的Time Warping，指的就是对时间序列进行的压缩或者延展以达到一个更好的匹对。
简单的例子
比如说，给定一个样本序列X和比对序列Y，Z：
1 | X：3，5，6，7，7，1 |
请问是X和Y更相似还是X和Z更相似？
DTW首先会根据序列点之间的距离（欧氏距离），获得一个序列距离矩阵 $M$，其中行对应X序列，列对应Y序列，矩阵元素为对应行列中X序列和Y序列点到点的欧氏距离：
X和Y的距离矩阵：
| X/Y | 3 | 6 | 6 | 7 | 8 | 1 | 1 |
|---|---|---|---|---|---|---|---|
| 3 | 0 | 3 | 3 | 4 | 5 | 2 | 2 |
| 5 | 2 | 1 | 1 | 2 | 3 | 4 | 4 |
| 6 | 3 | 0 | 0 | 1 | 2 | 5 | 5 |
| 7 | 4 | 1 | 1 | 0 | 1 | 6 | 6 |
| 7 | 4 | 1 | 1 | 0 | 1 | 6 | 6 |
| 1 | 2 | 5 | 5 | 6 | 7 | 0 | 0 |
然后根据距离矩阵生成损失矩阵（Cost Matrix）或者叫累积距离矩阵 $M_c$，其计算方法如下：
- 第一行第一列元素为 $M$ 的第一行第一列元素，在这里就是0；
- 其他位置的元素 $M_c(i,j)$ 的值则需要逐步计算，具体值的计算方法为 $M_c(i,j)=\text{Min}(M_c(i-1,j-1),M_c(i-1,j),M_c(i,j-1))+M(i,j)$ ，得到的 $M_c$ 如下：
| X/Y | 3 | 6 | 6 | 7 | 8 | 1 | 1 |
|---|---|---|---|---|---|---|---|
| 3 | 0 | 3 | 6 | 10 | 15 | 17 | 19 |
| 5 | 2 | 1 | 2 | 4 | 7 | 11 | 15 |
| 6 | 5 | 1 | 1 | 2 | 4 | 9 | 14 |
| 7 | 9 | 2 | 2 | 1 | 2 | 8 | 14 |
| 7 | 13 | 3 | 3 | 1 | 2 | 8 | 14 |
| 1 | 15 | 8 | 8 | 7 | 8 | 2 | 2 |
最后，两个序列的距离，由损失矩阵最后一行最后一列给出，在这里也就是2。同样的，计算X和Z的距离矩阵：
| X/Z | 2 | 5 | 7 | 7 | 7 | 7 | 2 |
|---|---|---|---|---|---|---|---|
| 3 | 1 | 2 | 4 | 4 | 4 | 4 | 1 |
| 5 | 3 | 0 | 2 | 2 | 2 | 2 | 3 |
| 6 | 4 | 1 | 1 | 1 | 1 | 1 | 4 |
| 7 | 5 | 2 | 0 | 0 | 0 | 0 | 5 |
| 7 | 5 | 2 | 0 | 0 | 0 | 0 | 5 |
| 1 | 1 | 4 | 6 | 6 | 6 | 6 | 1 |
和损失矩阵：
| X/Z | 2 | 5 | 7 | 7 | 7 | 7 | 2 |
|---|---|---|---|---|---|---|---|
| 3 | 1 | 3 | 7 | 11 | 15 | 19 | 20 |
| 5 | 4 | 1 | 3 | 5 | 7 | 9 | 12 |
| 6 | 8 | 2 | 2 | 3 | 4 | 5 | 9 |
| 7 | 13 | 4 | 2 | 2 | 2 | 2 | 7 |
| 7 | 18 | 6 | 2 | 2 | 2 | 2 | 7 |
| 1 | 19 | 10 | 8 | 8 | 8 | 8 | 3 |
所以，X和Y的距离为2，X和Z的距离为3，X和Y更相似。
定义
有一个具体例子作为帮助，我们再来定义DTW算法。
假设给定两个序列，样本序列 $X=(x_1,…,x_N)$ 和测试序列 $Y=(y_1,…,y_N)$，同时给定一个序列中点到点的距离函数 $d(i,j)=f(x_i,y_j)\geq 0$（一般为欧氏距离，实际上也可以是别的函数）。
那么DTW的核心在于求解扭曲曲线（Warping Curve）或者说扭曲路径，也就是点点之间的对应关系。我们表示为 $\phi(k)=(\phi_x(k),\phi_y(k))$，其中 $\phi_x(k)$ 的可能值为 $1,2,…,N$，$\phi_y(k)$ 的可能值为 $1,2,…,M$，$k=1,…,T$。也就是说，求出 $T$ 个从X序列中点到Y序列中点的对应关系，例如若 $\phi(k)=(1,1)$，那么就是说X曲线的第一个点与Y曲线的第一个点是一个对应。
给定了 $\phi(k)$，我们可以求解两个序列的累积距离（Accumulated Distortion）：
DTW的最后输出，就是要找到一个最合适的 $\phi(k)$ 扭曲曲线，使得累积距离最小，也就是损失矩阵的最后一行最后一列的值：
换句话说，就是给定了距离矩阵，如何找到一条从左上角到右下角的路径，使得路径经过的元素值之和最小。这个问题可以由动态规划（Dynamic Programming）解决（时间复杂度 O(N+M)），也就是上面例子中，计算损失矩阵的过程，实际上不需要把整个矩阵都求解出来，大致将对角线上的元素求解出来即可。
讨论
实际上，虽然这个算法简单，但是有很多值得讨论的细节。
约束条件
首先，路径的寻找不是任意的，一般来说有三个约束条件：
- 单调性：$\phi_x(k+1)\geq\phi_x(k)$ 且 $\phi_y(k+1)\geq\phi_y(k)$ 也就是说扭曲曲线不能往左或者往上后退，否则会出现无意义的循环；
- 连续性：$\phi_x(k+1)-\phi_x(k)\leq1$，即扭曲曲线不能跳跃，必须是连续的，保证两个序列里的所有点都被匹配到，但这个条件可以一定程度上被放松；
- 边界条件确定性：$\phi_x(1)=\phi_y(1)=1$，$\phi_x(T)=N$，$\phi_y(T)=M$，即路径一定从左上开始，结束于右下，这个条件也可以被放松，以实现局部匹配。
除此之外，我们还可以增加别的约束：
- 全局路径窗口（Warping Window）：$|\phi_x(s)-\phi_y(s)|\leq r$，比较好的匹配路径往往在对角线附近，所以我们可以只考虑在对角线附近的一个区域寻找合适路径（$r$ 就是这个区域的宽度）；
- 斜率约束（Slope Constrain）：$\dfrac{\phi_x(m)-\phi_x(n)}{\phi_y(m)-\phi_y(n)}\leq p$ 和 $\dfrac{\phi_y(m)-\phi_y(n)}{\phi_x(m)-\phi_x(n)}\leq q$， 这个可以看做是局部的Warping Window，用于避免路径太过平缓或陡峭，导致短的序列匹配到太长的序列或者太长的序列匹配到太短的序列。
步模式
实际上，这些步模式（Step Pattern）一定程度上涵盖了不同的约束，步模式指的是生成损失矩阵时的具体算法，例如在例子中使用的是：
准对称步模式。实际上还有很多其他步模式，不同的步模式会影响最终匹配的结果。关于不同的步模式，可以参见[2]第四章。常用的有对称，准对称和非对称三种。
标准化
序列的累积距离，可以被标准化，因为长的测试序列累积距离很容易比短的测试序列累积距离更大，但这不一定说明后者比前者与样本序列更相似，可以通过标准化累积距离再进行比较。不同的步模式会需要的不同的标准化参数。
点与点的距离函数
除了测试序列以外，DTW唯一需要的输入，就是距离函数 $d$（除了欧氏距离，也可以选择Mahalanobis距离等），所以不需要考虑输入的具体形式（一维或多维，离散或连续），只要能够给定合适的距离函数，就可以DTW比对。前面说到，DTW是对时间上的压缩和延展不敏感，但是对值的大小是敏感的，可以通过合理选取距离函数来让DTW适应值大小的差异。
具体应用场景
这里讨论两个具体应用DTW的可能场景：
分类
气象指数在旱季和雨季的样本序列分别为 $X_1$ 和 $X_2$，现有一段新的气象指数 $Y$，要判断该气象指数测得时，是雨季还旱季？
算出DTW$(X_1,Y)$ 和 DTW$(X_2,Y)$，小者即为与新测得气象指数更贴近，根据此作判断。
DTW就是一个很好的差异比较的工具，给出的距离（或标准化距离）能够进一步输入到KNN等分类器里（KNN就是要找最近的邻居，DTW能够用于衡量“近”与否），进行进一步分类，比对。
点到点匹配
给定标准语句的录音X，现有一段新的不标准的语句录音Y，其中可能缺少或者掺入了别的字词。如何确定哪些是缺少的或者哪些是掺入别的？
通过DTW的扭曲路径，我们可以大致得到结论：
DTW的输出是很丰富的，除了距离外，还提供了扭曲路径，可用于点到点的匹配，这个信息是非常丰富的，能够看到序列的比对，发现异常的序列。
优化算法（2012以前）
论文中还提到了利用多核并行化优化。因为在这里不涉及到具体数学优化算法，所以不做介绍。
去平方根操作
在 DWT 和 ED 计算中都有带有开平方根的操作，为了方便理解后边优化在这里就去掉了平方根，因为这两种累加计算是单调递增的，所以去掉平方根对排序相似度结果上没有任何影响。
使用下界
通过使用下界可以提前筛选掉部分不符合条件的候选结果。我们可以这样理解下界这个概念：如果想找到班里最高的同学，常规做法是准确测量全班同学的身高并按身高排序，再择测得最高身高数据的那个。这样做当然可行，但想得到最终结果我们必须测量每一个同学的身高，这个工作量是相对比较大的。但也可以换个方法考虑：比如通过目测估计，选出一个中等偏上高的同学作为标杆，所有目测比他矮的同学肯定都不是要找的最高的那个，所以就不用花时间准确测量这些目测都不达标同学的身高。这样就通过这个中等偏上身高的同学（下界）排除了一大批是最高那个同学的可能。最后只需要测量小部分比这个下界高的同学，就能很快找到最终结果。
计算下界有很多种不同的方法，有些虽然复杂度小，但筛选效果未必很理想，有些效果好，但是计算量却很高。本文中提到了下面两种被工业广泛应用的计算下界方法。
LB_Kim_FL
该算法其实是 LB_Kim 算法的一个在归一化处理后的简化版本。原始的 LB_Kim 会一共取了4个观测点，即：初始点，终止点，时序数据最大值，时序数据最小值。我们在这里提到的 LB_Kim_FL 只用了起始和终止值，因为在归一化后的极值在这里几乎不起任何作用，但是在省略它们后整个算法的时间复杂度却从 O(n) 降为 O(1)。复杂度的下降是省略了遍历整个数组寻找全局最大和最小值的过程而做到的。经实践证明，这种方法在大多数情况下十分地粗暴有效，但在起始值和终止值差别较大时会对筛选结果造成一定的误差。
计算公式为：
LB_Keogh
与简单粗暴的 LB_Kim 不同，LB_Keogh 计算时不只关注四个关键值，对查询子序列 Q 引入了 Upper Bound 和 Lower Bound 的概念。这种方法很好的解决了 LB_Kim 只关注起始值与终止值的问题，从而达到了更加精确地定义下界的目的，但由于它需要遍历整个样本才能完成计算，故复杂度为O(n)。
具体计算步骤如下：
定义一个时间序列子序列 $Q:\{q_1,\dots,q_n\}$ 以及窗口值 $r$ 且在其范围内寻找 max 和 min 分别作为 Upper Bound，Lower Bound。
确定 Upper Bound ：$U_i=\text{max}(q_{i-r},q_{i+r})$
确定 Lower Bound ：$L_i=\text{min}(q_{i-r},q_{i+r})$
利用公式计算 LB_Keogh：
引用一位大牛对这个算法的理解： LB_Keogh 方法首先求出 Q 序列的上下包络线，然后对data序列 ( C ) 与上下包络线进行比较，如果不在上下包络线的范围内，就对该点与对应的包络线上的点求欧几里得距离（此处的欧几里得距离主要是指Y轴上的距离，并非二维），最终求和得到误差，与之前得到的误差进行比较，整条data中误差最小的就是目标序列。
Early Abandoning of ED and LB_Keogh
我们在之前讨论的 Lower Bound 剪枝其实就是一种 Early Abandoning。在计算中涉及到一个 best-so-far 的概念，在初始化时这个值被设置为无穷大。本着一旦找到更好的结果就替换的原则，然后用下界（Lower Bound）与这个值对比，如果下界大于它，那么就可以在此终止后面计算且筛选掉这条不符合条件的数据，反之则取代（更新）当前 best-so-far。
Early Abandoning of DTW
对于计算完 LB_Keogh 再计算 DTW 的复杂度也有相应的优化方法，在原文中作者也介绍了他们自己的方法。文章指出在 LB_Keogh 已经计算完一次后，这个结果可以被重复使用在计算 DTW 中，以省略重复计算带来的开销。
如上图所示，在计算 DTW 时从 K = 0 位置开始，计算至 K = 10 结束，这时我们得到了部分 DTW 距离累加的结果。在此可以用该部分结果与之前 LB_Keogh 的 K = 11 后的结果拼接，并把这个拼接的结果当做 Lower Bound 对比 best-so-far 决定是否要提前终止计算。
这种混合 Lower Bound 计算可以表示为：
本文的优化算法
Searching and Mining Trillions of Time Series Subsequences under Dynamic Time Warping
Early Abandoning Z-Normalization
作者在这里吐槽了一下之前人研究成果的疏忽，他指出尽管计算归一化的开销远大于计算欧氏距离的开销，但迄今为止（2012年）还没人提出过对这个归一化过程做优化，并在此引出他们的创意：在计算欧式距离或者 LB_Keogh 的过程中同时进行 online Z-Normalization，而不是分别做。而且一旦发现 C 子序列已经可以排除在最优结果之外了也没有必要继续对该序列做剩余归一化的操作。这个方案可以节省的对时间序列做完整归一化操作带来的巨大开销，因为这里的归一化操作只是针对部分子序列的计算而不是全部。
常规来讲想要对某个子序列进行 Z-Normalization 需要先计算其整体的均值和标准差。即：
这里的标准差计算和公式略有不同，因为作者展开了平方差项。然后用遍历样本后求得的值带入归一化公式：
用这种常规方法计算均值和标准差是需要遍历整个子序列的，这个造成了一定的开销。在此原文作者提出了他改良过的方法，专门计算子序列归一化涉及到的均值和标准差的计算。即：
且在计算过程中储存中间输出的部分子序列加和结果和加和平方结果，这样做可以避免重复计算。这里的 $m$ 指的是时间子序列的长度，$k$ 为整个时间序列中从起点到截取的第 $k$ 个时间段序列。
此处需要用图解释 $m$ 和 $k$ 在整个时序数据中的位置 TBD
论文中，作者提到此处存在浮点计算在上亿条数据下回出现累加误差的问题，并提出这里每比对100万子序列后，就要进行一次完整的Z-normalization。
Reordering Early Abandoning
前面 Early Abandoning 的计算方式都是从子序列的起点自左向右进行计算的。作者在这里提出一种先快速找到子序列 Q 和 C 之间差值之和最大的部分，然后以它为中心计算 LB 来判断这个子序列是否大于best-so-far值，从而可达到降低计算成本的目的。引用下原文的图：
左图中计算的顺序是从左向右的，在计算到第9个时间点时发现 Lower Bound 已经超过了 best-so-far 值，然后终止了计算并筛选掉 C。 右图中将顶点放在离 Q 峰值很近的位置，且交替的从初始点两侧分别取值累加，这样只运算了 5 次就超过了 best-so-far 值，这里就可以提前终止运算并排除当前 C是最优解的可能性。紧接着作者引出他的优化问题，如何选择合适的初始点来做优化，并给出了他的经验：在已经做完归一化的时序子序列中，由于本文讨论的是 Z-Normalization 也就是说归一化操作后的序列均值为0，方差为1 的高斯分布。使用离均值 0 最远距离的点一般应该是对计算距离贡献最大的点，也是起始的最佳位置，并在文章后边的实验中给予了证实。
Reversing the Query/Data Role in LB_Keogh
在讨论 LB_Keogh 算法时我们是围绕 Q 进行展开的，而不是全部的时间序列 C。这样只需要对 Q 进行一次计算，而不是针对每一个可能的结果进行逐次计算。但有时候对 Q 的计算有可能效果并不理想，这里作者提出了可以 just-in-time 的互换 Q 和 C 的角色来加速计算。
上文讲过 LB_Keogh 的计算复杂度是 O(n) 而 DTW 是动规问题，高于 O(n) 。比起算所有的 C 的 DTW值，这里计算所有 C 的 LB_Keogh 还是比计算 DTW 节省开销的。具体看后边逐层优化，这个优化算法是在比较后面才被用到的。越往后剩下的 C 子序列越少，因为大部分不符合最优解的 C 已经被简单的 LB 方法排除出去了。
Cascading Lower Bounds
原文作者在50个数据集上使用18种不同的计算 Lower Bound 的方式测试了筛选效果，并绘制了下图。
从图中可以看出每个 LB 算法都可以从计算复杂度和筛选宽松度来评价。在此基础上作者引出了“层级筛选”的概念，即：先使用复杂度较低的算法快速过滤掉一批不符合最优结果的候选 C，再逐渐使用复杂度高且比复杂度低算法更严谨的筛选算法过滤掉剩余数据。准确的说：先使用 O(1) 复杂度的 LB_Kim_FL，然后对 Q 进行 LB_Keogh_EQ 剪枝，此处的算法复杂度可能介于O(1) 与 O(n) 之间，如果该剪枝算法执行完后 LB 仍没有达到 best-so-far 再用 3.3 章节提到的互换 Q 和 C 的 LB_Keogh_EC 算法辅助剪枝效果，最后使用 提前终止 DTW 的算法处理剩余数据。利用这个顺序剪枝可以节省99.9999%的DTW算法的时间开销。
总的来说，这个层级筛选的方法是整个 Paper 的精华所在，之前作者在考虑到计算复杂度和筛选精准的前提下分别介绍了几种高效计算 LB 的方法，最后从18种中选出四个组合效果最好的，并按复杂度从小到大依次就情况使用。
总结
通过使用URC优化的 DTW 可以大大提升搜索速度，这方便了大文本时间序列相似度对比在工业中的落地。在写完这个总结时原作者脑子里从日常的需求中浮现出如下应用场景：
- 实时系统监控： 当今的分布式系统例如 Hadoop， Kubernetes 等每天都在工业界产生大量时序日志数据。如果可以将这些数据导入分布式时序数据库中的话，可以通过对特定指标的搜索来更好地做日志分析，甚至由于改善后的速度提升做系统的实时监控。
- 股票 K 线分析：假设可以观察股票 K 线的走势预测未来股票的涨跌。那么可以通过查找历史数据中某个盈利股票在暴涨前的一段时序数据，将这个数据作为 Q 去搜索当前所有股票 C 与这个股票最相近的一个，如果他们的相似度高，也许可以预测这个股票在未来的一段时间会涨。
- 客服中心 Trend 预测： 以本人所在保险行业为例，客服中心（Call Center）的工作量与时间周期，突发性事件（如：自然灾害，政治变动，传染病）有极大的相关性。一旦某事件发生，会引起 Hotline 在未来的一段时间段打入电话数量暴增，导致需要同时在线的客服人员数量翻倍增长。为了节约成本，优化安排客服人员工作量，可以使用 DTW 这类算法对历史数据进行搜索，如发现某个 Trend 会在未来的一段时间造成工作量非正常的变化，则可以提前组织相应工作以便做到节约成本（如：通过优化客服人员数量、工时）和提高客服质量的目的。
- 销售 Trend 预测： 方法同3，可以预测下个阶段销售数量，以便优化供应链管理。
- 工业中可以结合分布式流式处理中间件 Kafka 以及 Kafka Streams，Spark Streaming，Flink 等分布式框架，结合时序 NoSQL 数据库实现海量数据高并发实时监控/预测。
参考文献
- Giorgino, Toni. “Computing and visualizing dynamic time warping alignments in R: the dtw package.” Journal of statistical Software 31.7 (2009): 1-24.
- Rabiner, Lawrence R., and Biing-Hwang Juang. Fundamentals of speech recognition. Vol. 14. Englewood Cliffs: PTR Prentice Hall, 1993.
