DTW算法挖掘亿万级时序数据，其优化能耐几何？

概述

常见的距离度量尺度有：欧式距离，马氏距离、曼哈顿距离。

欧式距离在n维空间中表示为：

欧几里得距离

常用于计算两个等长序列的相似性，序列可以包含股票、DNA等。对于不等长的序列，计算相似性关键是在于计算之前或者过程中，找到两个时序数据上点的对应关系，因为它不像等长数据已经存在了数据点的一一对应关系。

等长数据，比如：比较近10天两只开盘股票走势k线图
不等长数据，比如：相同时间内不同抽样频率（Hz）的心电图、两段“麻烦请开门”的语音音频

因此，日本学者Itakura最早提出 Dynamic Time Warping（下文简称DTW，中文常翻译为“动态时间规整”）算法，它出现的目的也比较单纯，是一种衡量两个长度不同的时间序列的相似度的方法——在对齐两个序列的过程中通过定义的距离计算公式计算序列的相似度。其应用广泛，主要在于模版匹配，如孤立词语音识别、手势识别、DNA序列配对等。

这里，距离计算公式包括（但不限于）上文提到的欧式距离。

DTW通过把时间序列进行延伸和缩短，来计算两个时间序列性之间的相似性，规整后两个序列之间的距离应该在所有对齐方案中是最小的，即达到feature-to-feature的对齐。

原理

'dtw算法原理'

以上面的Q,C序列为例，DTW算法的求解过程的直观理解为构建一个nxn的矩阵(此处假设Q和C的时间序列长度都为n——当然也可以不等长——矩阵中的元素(i,j)表示序列Q的时间点qi和序列C的时间点cj之间的欧式距离），目标是在矩阵中找到一条从(0,0)到(n,n)的路径，使得路径上的所有元素之和最小。

“矩阵中的元素(i,j)表示序列Q的时间点qi和序列C的时间点cj之间的欧式距离”，下图进行说明。

1.png

“目标是在矩阵中找到一条从(0,0)到(n,n)的路径，使得路径上的所有元素之和最小”，下图进行说明。

2.png

这里提及dtw的几个约束：

local constraint（在矩阵中表示为定义的“步”的方向都得朝向右上角，每一步都得离终点更近，否则会导致crossing line）
global constraint（不可跨越一定限度的数据点进行对齐，否则会导致对齐密度不均衡）
start&ending contraint（头尾数据点各自对齐）
weight（路径权重的设置，平衡业务偏好和local内的距离偏好）
distance（距离度量方式）

允许在对齐过程中，有些点被跳过，没有被对齐——取决于定义的“步”。
由于矩阵中路径上的每个点都能分解为一个“多点到达”的子问题结构，因此dtw就是通过动态规划法（dp，dynamic programming）进行求解的一个例子。

在后面给出的DTW算法的python简单实现中，通过在循环中约束i和j的关系实现global constraint——通过简单画图可以了解。

问题与优化

那简单介绍完DTW的原理之后便引入了这篇论文《KDD2012 Best Paper-"Searching and Mining Trillions of Time Series Subsequences under Dynamic Time Warping”》[1]。

上文提到，DTW在多个领域都有所应用。随着互联网的到来和数据量的爆发（很多生产环境数据量早已突破万亿级别，而学术界仍停留在百万、十亿级数据集的研究上），原始DTW的实现弊病暴露，场景的应用对算法的性能提出了更高的要求，而该论文的核心正是通过对其他论文的review，现有优化方式的review以及对提出的计算优化方案的review来告诉我们，优化后的DTW（论文称为UCR suite）仍然是最强最快的时序相似度计算方式。

论文首先指出以下假设或事实：

标准化（Z-normalization）非常重要，不仅需要在整个数据集上做，在计算两序列相似度之前，还需要在两个序列上单独做。
在巨量数据集数据库中检索变长序列，理论上可行，但实际上不可行。

已知的优化方法：

计算过程中欧式距离使用平方代替平方根，直到获得最小距离（的平方）时再开方获得最终结果。
使用lower bound技术，伪代码思想如下：

Algorithm Lower_Bounding_Sequential_Scan(Q)
best_so_far = infinity; 
for all sequences in database 
    LB_dist = lower_bound_distance(Ci,Q); 
        if LB_dist < best_so_far
            true_dist = DTW(Ci,Q);
            if true_dist < best_so_far
                best_so_far = true_dist;
                index_of_best_match= i;
            endif
        endif
endfor

lower bound具体计算方法有多种，如稳重提到的LB_kim，LB_keogh。

在计算欧式距离或lower bound时采用早停技术。
因为计算DTW真实距离（最优路径）时采用 DTW(Q1:K,C1:K) + LB_Keogh(QK+1:n,CK+1:n)作为实时的距离（为真实距离的lower bound，K为任一中间数据点序号），在此采用早停技术。
使用多核机器并行计算（众所周知）。

论文提出的新的优化方法：

标准化比计算欧几里得距离的耗时还要长一些，因此考虑在标准化过程中结合计算欧式距离或lower bound，引入早停技术。
对时序进行重排序

We conjecture that the universal optimal ordering is to sort the indices based on the absolute values of the Z-normalized Q.
For this we simply take each Ci and sort them, largest first, by their sum of their contributions to the Euclidean distance.
We compared this empirically optimal ordering with our predicted ordering (sorting the indices on the absolute values of Q) and found the rank correlation is 0.999.

对备查序列建立包络，而不是查询序列。
将多种lower bound计算方式串联。
- 因为lower bound技术实证有助于加速搜索过程，因此在研究界对此的研究很繁荣。每一种的紧度（tightness）和计算速度都不尽相同。下图为证：[图片上传失败...(image-6d35f9-1571736906097)]
  可见，紧度和计算复杂度是成正比的，因此没有最佳的技术，只有最合适的技术。作者的想法是将这些技术从复杂度低到高串联起来，如果计算完和lower bound计算完距离后没有超过目前最小距离（the best-so-far），那就接着计算上述方案（早停版的DTW）——不管是O(1)还是O(nR)，都远小于O(n2)。
另外，对于存在重抽样的序列（如[2.34, 2.34, 2.34, 2.01, 2.01, 2.01, 1.99, 1.99, 1.99,... ]），除了使用以上UCR-dtw方法，作者有以下建议：
1. 先对原数据的子集按照比例n下采样，使用UCR-dtw算法计算出距离rD；
2. 再在完整数据集上运行，提前设置参照距离为rD*10。
3. 如此之后，在特定数据集上能提升3倍的速度。

未完待续，欢迎留言交流。
此处可获取代码，结合理解