基于DTW探索挖掘万亿级时间序列
【前言】：时间序列数据大量充斥在人们的生活当中, 如医疗，金融，科学等诸多领域。而时间序列的数据挖掘需要相似度比较作为子程序。常用方法为欧氏距离ED(Euclidean distance)和动态时间规整DTW(Dynamic Time Warping)。因为考虑到时间序列存在长短不一，时间步没有对齐等问题， 所以DTW在时间序列探索过程中效果更好。该文章是2012年发表的ACM SIGKDD的文章，主要基于DTW进行了优化(ED也同样优化)，该文章称之为UCR套件。优化后的算法对比该paper之前算法速度有非常巨大的提升，并获得了SIGKDD 2012年最佳论文奖。该文章作者为Thanawin Rakthanmanon，Bilson Campana等参考文献：
一. 要解决的问题
本文解决的问题--探索挖掘万亿级(Trillions)数据。即在一个长序列中查找与目标序列最相近的一段子序列。其中长序列为，代表trillion级别数据；目标序列即需要查找的quary； 代表长序列中截取的子序列，也可以称为候选序列C。在长序列找到一段使其与最相似。具体示例:
: 4, 1, 1, 3, 5, 8, 9, 9, 7, ...........], 长度为n.
: , 长度 m=6
从i=0->n-m， https://www.zhihu.com/equation?tex=C_%7Bi%7D%3DT_%7Bi%2C+i%2Bm%7D , 依次比较https://www.zhihu.com/equation?tex=dist%28Q%2C+C_%7Bi%7D%29 距离, 搜索出与距离最近的一段子序列。
这个相似度匹配的过程即可以使用前言中提到的ED和DTW(其他距离比较也可以尝试使用)。DTW效果更好一些，但其运算速度比使用ED要慢。因此文章中提出了相应的解决办法，并命名为UCR suite。另外trillion是非常大体量的数据, 之前论文中从未涉及到。
另外,解决问题过程中文章中提出了三个假定的声明:
1> 时间序列子序列必须标准化(normalize)， 即上文提到的或需要normalize。当然Q也是需要normalize的，Q的长度相应较短，不影响计算效率。因为对两个时间序列进行有意义的比较，必须对他们进行规范化， 以保证比较的尺度是统一的。作者在论文中针对“枪支分类为题”做了验证，改变偏移量和比例5%， 错误率至少提升50%；
2> DTW是最好的时间序列相似度度量方法. 后面会详细介绍DTW；
3> quary是任意长度的, 无法被索引；
二. DTW算法原理.
欧氏距离, 只能用于两个相等长度的时间序列之间的距离， 且时间步之间需要一对一匹配:
https://www.zhihu.com/equation?tex=ED%28Q%2CC%29+%3D+%5Csqrt%7B%5Csum_%7Bi+%3D+1%7D%5E%7Bm%7D%7B%28q_%7Bi%7D-c_%7Bi%7D%29%5E2%7D%7D ,
ED示意图如左上图，其算法复杂度O(n)。但当两个序列不等长， 或序列之间的时间步不是一一对应的情况下 ，就难以比较两个序列的距离，例如相同时间段内按照不停的采样频率采样得到的数据，两个人说的相同内容的语音等。
图1 ED和DTW示意图
因此，DTW被提出来，它出现的目的也比较单纯，是一种衡量两个长度不同的时间序列的相似度的方法——在对齐两个序列的过程中通过定义的距离计算公式计算序列的相似度。如左下图。DTW算法定义如下，假设我们有两条序列 Q（query）和 C (candidates), 分别长度为 https://www.zhihu.com/equation?tex=n 和 https://www.zhihu.com/equation?tex=m ，他们不一定相等。
https://www.zhihu.com/equation?tex=Q%3Dq_1%2C+q_2%2C......%2Cq_i%2C......%2Cq_m
https://www.zhihu.com/equation?tex=C%3Dc_1%2C+c_2%2C......%2Cc_j%2C......%2Cc_n+
为了比对Q和C，组成一个n by m的矩阵，每两个对齐后的元素之间的距离下标定义为 https://www.zhihu.com/equation?tex=%EF%BC%88i%5E%7Bth%7D%2C+j%5E%7Bth%7D%EF%BC%89 距离定义为比如 https://www.zhihu.com/equation?tex=d%28q_i%2C+c_j%29%3D%28q_i+-+c_j%29%5E2 。这个矩阵称为距离矩阵,有一条 warping path 从(0,0)->(m,n)穿越这个矩阵， 定义为。的第 https://www.zhihu.com/equation?tex=K 个元素，定义为 https://www.zhihu.com/equation?tex=w_k%3D%28i%2C+j%29_k。
DTW需要满足相应的约束：
1> 路径必须开始和结束在对角的单元
2> warping路径的步骤被限制在相邻的单元
3> warping路径上的点必须在时间上连续.如图1右侧,每个方格代表距离, 红色线代表最优的W路径。
则DTW即求从矩阵左下角到右上角的最优路径：
https://www.zhihu.com/equation?tex=DTW%28Q%2C+C%29+%3D+min+%5Cleft%5C%7B+++++%5Csqrt%7B%5Csum_%7Bk-1%7D%5E%7BK%7D%7Bw_k%7D%7D+%5Cright.+++
这条最优路径可以通过动态规划的算法求解， 每个时间步的累计距离 ,即称为损失矩阵MCMC, 其中r(i,j)为损失矩阵中的对应元素。
https://www.zhihu.com/equation?tex=%5Cgamma%28i%2C+j%29+%3D+d%28q_i%2C+c_j%29+%2B+min+%5C%7B+%5Cgamma%28i-1%2Cj-1%29%2C+%5Cgamma%28i-1%2C+j%29%2C+%5Cgamma%28i%2C+j-1%29%5C%7D
整体算法复杂度 O(mn)。
为减少计算， 通常会限定偏离对角的长度R， 以达到优化算法的目的。也有可以采用下界方式进行优化。具体代码可参照my github。
以下为DTW计算的示例:
我们假设 Q: ,m=7
C :,n=8
则计算构建距离矩阵,其中红色为DTW最优路径, Q和C的DTW距离为所有红色方框的和, DTW(Q,C) = 1:
下面我们采用动态规划的方法求解r(i, j)，并构建MCMC矩阵如下，则DTW(Q,C)值为图中绿色位置的值, DTW(Q,C) = r(m-1,n-1):
三. 已知的算法优化(该论文发表之前)
1.使用平方距离代替均方根
是否进行开根号，不影响搜索过程中距离的排序。该过程只在算法内部使用，输出距离的时候会进行开根号呈现给用户。(该方法适用于DTW和ED算法)。
2. 采用下界(Lower Bounding)
由于搜索的数量 https://www.zhihu.com/equation?tex=T 很大，因此，针对每一个候选子序列与Q进行DTW距离计算时非常耗时的。因此提出了下界的思想，使用计算成本较低的下界来删除没有前途的候选子序列。该trick仅适用于DTW距离计算。伪代码如下：
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 Bounding具体的计算方法这些年学界有很多积累。 这里只提文中提到的两种。 LB_Kim 和 LB_keogh。
1> LB_Kim
LB_Kim即取两个序列Q和Ci的起始点，终止点，最大点，最小点之间距离的最大值。公式：
https://www.zhihu.com/equation?tex=LB%5C__%7BKim%7D%28Q%2C+C_%7Bi%7D%29+%3D%5Cmax_%7Bk%3DA%2CB%2CC%2CD%7D%5C+d%28q_%7Bk%7D%2C+c_%7Bk%7D%29
该算法时间复杂度为O(n), 因为Q和 在进行比较之前都进行过标准化, 因此最大值B点和最小值C点的差值相对较小。本文中直接忽略，提出了LB_Kim_FL方法。这样时间复杂度降低为O(1)。公式如下：
https://www.zhihu.com/equation?tex=LB%5C__%7BKim%5C_%7B%7DFL%7D%28Q%2C+C_%7Bi%7D%29+%3D%5Cmax_%7Bk%3DA%2CD%7D%5C+d%28q_%7Bk%7D%2C+c_%7Bk%7D%29
2> LB_keogh
第一部分为Q的{U， L} 包络曲线（具体如图）， 给Q序列的每个时间步定义上下界。 定义如下:
https://www.zhihu.com/equation?tex=U_i+%3D+max%28q_%7Bi-r%7D%3Aq_%7Bi%2Br%7D%29 ， https://www.zhihu.com/equation?tex=L_%7Bi%7D+%3D+min%28q_%7Bi-r%7D%3Aq_%7Bi%2Br%7D%29，
https://www.zhihu.com/equation?tex=s.t%3Am-r+%5Cle+i+%5Cle+m%2Br+
U 为上包络线，就是把每个时间步为Q当前时间步前后r的范围内最大的数。L 下包络线同理。
那么LB_Keogh定义如下：
https://www.zhihu.com/equation?tex=LB%5C_Keogh%28Q%2C+C%29+%3D+%5Csqrt+%7B+%5Csum_%7Bn%7D%5E%7Bi%3D1%7D+++++%5Cbegin%7Bequation%7D+++++++%5Cleft%5C%7B++++++++++++++++%5Cbegin%7Barray%7D%7Blr%7D++++++++++++++++%28c_i+-+u_i%29%5E2++%26+%5Ctext%7Bif+%7D+c_i+%5Cgt+u_i%2C+%5C%5C++++++++++++++%28c_i+-+l_i%29%5E2++%26+%5Ctext%7Bif+%7D+c_i+%5Clt+l_i+%2C%5C+%5C%5C+++++++++++++++0+%26+%5Ctext%7Botherwise%7D++++++++++++++++++%5Cend%7Barray%7D+++++++%5Cright.+++++++%5Cend%7Bequation%7D+++%7D
如图阴影部分对DTW距离的贡献。
3. Early Abandoning of ED and LB_Keogh
ED和LB_Keogh的早弃思想。在ED和LB_Keogh的计算过程中, 如果发现当前累加的距离值已经超过当前最优值(best so far)时,可以停止计算。如图，当计算到箭头位置时，累积的距离已经超过当前最优值，则停止计算。
4. Early Abandoning of DTW
DTW的早弃思想。当计算全部的LBkeogh后发现还是需要计算该子序列C的DTW值, 这时DTW计算过程中我们依然可以采用Early Abandoning的思想。由LB_keogh的算法可知, 对于相同长度的序列, DTW_dist≥LB_keogh。因此我们可以结合DTW和LB_keogh思想达到Early Abandoning的目的。 可以计算DTW(Q1:K,C1:K) + LB_Keogh(QK+1:n,CK+1:n)，为DTW(Q, C)在k时刻的下界，即在k时刻DTW(Q1:K,C1:K) + LB_Keogh(QK+1:n,CK+1:n) >DTW(Q, C)时可以直接停止当前子序列的相似度计算。
5. 采用多进程思想加速
四. 本文提出的优化思想-UCR suite
论文中共提出四种优化思想, 基于ED/DTW, 统称他们为UCR suite.
1. Early Abandoning Z-Normalization
前面提到了对于每一个候补子序列都需要进行Z-Normalization。Z-Normalization可以理解为每个时间部的z值(即或)需要进行normalization。
Z-Normalization的过程比计算ED距离的时间复杂度稍长一点，因此我们可以递归的去计算Z-Normalization，当ED/LB_keogh计算早弃(Early Abandoning)时, 也是对Z-Normalization的Early Abandoning.
Z-Normalization计算公式:
https://www.zhihu.com/equation?tex=%5Cmu+%3D+%5Cdfrac%7B1%7D%7Bm%7D%5Csum%7Bx_i%7D , https://www.zhihu.com/equation?tex=%5Csigma%5E2+%3D+%5Cdfrac%7B1%7D%7Bm%7D%5Csum%7Bx_i%5E%7B2%7D+-+u%5E2%7D .
其中方差计算参照公式 https://www.zhihu.com/equation?tex=D%28X%29+%3D+E%28X%5E2%29-%5BE%28X%29%5D%5E2 。
在相似度查找过程中，每个子序列在与查询进行比较之前需要进行normalization。子序列的均值可以通过保存当前时间k之前m的时刻的和求得，即 https://www.zhihu.com/equation?tex=%5Csum_%7Bi%3D1%7D%5E%7Bk-m%7D%7Bx_i%7D 。子序列的平方和可以用类似的方法计算。因此，Z-Normalization递归更新成为了可能， 递归计算时采用公式如下：
https://www.zhihu.com/equation?tex=%5Cmu+%3D+%5Cdfrac%7B1%7D%7Bm%7D%28%5Csum_%7Bi%3D1%7D%5E%7Bk%7D%7Bx_i+-+%5Csum_%7Bi%3D1%7D%5E%7Bk-m%7D%7Bx_i%7D%7D%29，https://www.zhihu.com/equation?tex=%5Csigma%5E2+%3D+%5Cdfrac%7B1%7D%7Bm%7D%28%5Csum_%7Bi%3D1%7D%5E%7Bk%7D%7Bx%5E2%7D+-+%5Csum_%7Bi%3D1%7D%5E%7Bk-m%7D%7Bx%5E2%7D%29+-+%5Cmu%5E2
2. Reordering Early Abandoning
对序列Q进行重新排序以达到Early Abandoning效果。如左图，顺序进行ED/LB_keogh时需要花费9次累加计算才能超过当前最优值,而右图仅需要5次计算即可超过当前最优值，从而放弃当前的候补子序列C。
因此可以考虑在 和距离较大的位置开始计算，那么这个最大的值怎么找呢？考虑到已经对Q和进行标准化，其分布应该服从高斯分布, 即数据均值为0。因此取Q中距离均值0最远的位置，可能会使 和 更大。这个想法是直觉的，无法证明。在实验中有一定效果。
3. Reversing the Query/Data Role in LB_Keogh
之前提到的LB_Keogh都是基于Q求出上下界 https://www.zhihu.com/equation?tex=U_Q , https://www.zhihu.com/equation?tex=L_Q+ , 然后计算LB_Keogh, 也可以叫做LB_Keogh_EQ, 同理我们也可以计算候补子序列的上下界, 并计算LB_Keogh_EC。两种方式得到的结果是不同的。 我们可以优先检查LB_Keogh_EQ，如果LB_Keogh_EQ<当前最优距离值时, 再检查LB_Keogh_EC。
4. Cascading Lower Bounds
级联下界思想, 即可以采用多个下界, 按照其时间复杂度和检出能力逐次使用。如果都不能检出当前子序列, 则再进行计算dist(Q, C).
如上图, 不同Lower Bounds算法对应的时间复杂度和紧密程度(检出Ci的能力)。作者依次采用了LB_KimFL, LB_keoghEQ, LB_keoghEC 三种下界算法.
五. 实验内容
[*]基于Random Walk date的验证
上表为基于Random Walk date的验证结果, query长度为128。可以看到, 不论是哪个量级的数据，UCR-ED和UCR-DTW在查询速度上均有较大的提升(10倍左右)。
随着Query序列长度的增加相应查询时间的变化。作者稳重还提到UCR-DTW比UCR-ED的查训时间比例是1.18, 但是并没有找到相应的实验数据。
2. 基于心电图和DNA等进行了实验验证, 且均取得了很好的效果.
六. 文献阅读总结
个人认为该文章较大的创新点:1>提出了normalization过程递归计算, 可适用于各种下界裁剪;2>提出了采用级联下界的思想, 其中LB_kim_FL具有O(1)的时间复杂度, 能够很快的过滤一大部分子序列.我觉得这两个创新也是该文章查找速度有如此大提升并能适用亿万级数据的根本原因.
附录:
[*]作者备注了code地址, 为C语言编写:
Project Website: http://www.cs.ucr.edu/~eamonn/UCRsuite.html
2. 个人github整理了DTW, LB_kim, LB_keogh的python实现, 欢迎参考.
https://github.com/Taiyx/UCR_DTW-ED
页:
[1]