UCRDTW和UCRED模型详解

“UCR-DTW”和“UCR-ED”模型详解

DTW(Dynamic Time Warping)及已知的优化策略

计算两个时间序列 Q 和 C 之间的相似度，常用的度量方法是欧式距离 (ED)，计算公式如下图(1-1) 所示：

“UCR-DTW”和“UCR-ED”模型详解

从上图可以看到欧式距离的局限性：欧式距离通过建立两个序列之间的一一对应关系，使得 Q 和 C 之间的波峰没有对齐，因此计算得到的序列相似度存在较大的偏差。DTW 算法则可以很好的解决这个问题。

大部分情况下，两个序列整体上具有非常相似的形状，但是这些形状在时间轴上并不是对齐的。所以在计算两个序列的相似度之前，需要将其中一个 (或者两个) 序列在时间轴上进行 warping，使得两个序列波峰更好的对齐。

DTW 就是实现这种 warping 的一种有效方法。换句话说，DTW 算法通过找到两个序列之间的另外一种非一一对应的映射关系，这个映射关系也称为 warping path。以上面的 Q,C 为例，得到的对应关系如下图 (1-2) 中的灰色线段所示:

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

如下图 (1-3) 所示，红色标识的路径即为 warping path：

目前已知的 DTW 顺序搜索的四种优化方法如下:

去掉平方根计算
通过使用 lower bounds 来剪枝，因为 lower bounds 的计算时间复杂度都小于 DTW 的时间复杂度。比如：
LBKimFLLBKimFL 的时间复杂度为 O(1)
本文在实现时，由于对时间序列进行标准化后，时间序列数据中的最大和最小值对于整个 lower bound 距离贡献较小，因此，去掉原来 LB_kim 算法 (时间复杂度为 O(n)) 中提取的四个特征点中的最大值和最小值，使得时间复杂度降为 O(1)。但是，实现中为了使此策略发挥最大作用，作者还提取了第 2,3 和倒数第 2,3 个时间点，来进行级联剪枝。(详情参考算法实现 lb_kim_hierarchy 方法)

基于这条优化策略，作者提出的 Reorder early abandoning 可以进一步降低计算成本。

即在计算 ED 或者 LBKeoghLBKeogh 的时候，如果当前的两个序列的时间点 (1,k)(注：k<=|Q|) 之间的差值平方之和，大于当前两个序列最小距离值 best-so-far，那么可以提前结束 Q 和 C 是否相似的判断。计算过程如下图 (1-4) 所示：

4.Early Abandoning of DTW

计算完整的 LBKeoghLBKeogh 的值，仍然需要计算完整的 DTW 的值，我们可以通过利用部分的 LBKeoghLBKeogh 的值来减少 DTW 的计算量。

比如，先从左到右计算时间点 [1,k] 的 DTW 值，然后在时间点[k+1,n]，复用前面计算好的 LBKeoghLBKeogh 值。最终得到的距离值依然为完整的 DTW 值的 lower bound。这样的话，我们就可以使用 stop early 策略，每计算当前时间点的 DTW 值，就可以复用前面计算好的 LBKeoghLBKeogh 来获得整个子序列的 lower bound 值。

通过比较这个 lower bound 和当前最小距离值 best-so-far 进行比较，如果当前的 lower bound 值大于 best-so-far，那么可以提前结束 DTW 的计算。

这种方式的直观表示如下图(1-5)：

“UCR-DTW”和“UCR-ED”模型详解

以上，介绍完了前人对 DTW 的四种优化方案。

还有一种已知的提升 DTW 的计算速度的策略即使用多核计算资源。

UCR Suite 的优化策略

相关概念及定义

定义 1:

时间序列 T 是一个有序列表：T=t1,t2,⋯,tm。然而源数据是一个很长的时间序列，我们最终需要把它与一个更短的子序列进行相似度比较。

定义 2:

子序列 Ti,k 是时间序列 T 中的一个子序列，它起始于 ti，长度为 k，即:

Ti,k=ti,ti+1,…,ti+k−1,1≤i≤m−k+1。

这里，我们把 Ti,kTi,k 记为 C，作为与 query Q 比较相似的候选子序列。令 Q 的长度为 |Q|=n。

定义 3:

Q 和 C 之间的欧式距离 (|Q|=|C|) 定义为 (公式 1)：

路径 P 的第 t 个元素定义为 pt=(i,j)tpt=(i,j)t，则我们可以把 warping path 表示为 (公式 2)：

P=p1,p2,⋯,pt,⋯,pT,n≤T≤2n−1

优化策略:

1.Early Abandoning Z-normalization

在计算 DTW 距离之前都需要对 Q 和 C 进行标准化处理，但是对整个数据集进行标准化的复杂度太高了，因此这里使用 online Z-normalization，这样的话就可以采用 early stop 的策略来提前结束 normalization 的计算。

首先计算序列 C 的均值和方差的公式如下所示(公式 3)：

当使用 online Z-normalization 的时候，当前遍历到源序列 T 中的第 k 个时间点，所计算得到的时间点元素累加和以及时间点元素的平方累加和表示为(公式 4)：

那么对于 k -m+1 到 k 之间的这 m 个时间点对应的均值和方差的计算公式如下所示(公式 5)

所以，基于 online Z-normalization 的 abandon normalization early 策略伪代码如下图 (1-7) 所示：

论文中，作者提到此处存在浮点计算误差累加的问题，这里通过每比对完 100 万子序列，就进行一次完整的 Z -normalization，进而消除误差累加的问题。

Reorder early abandoning
前面 early abandoning 策略的计算方式都是从子序列的第一个时间点开始，自左向右进行计算的。本文提出一种策略是先快速找到 Q 和 C 之间差值之和最大的子序列，然后根据它来判断这个子序列是否大于 best-so-far 值，从而可达到降低计算成本的目的。这两种顺序的计算成本对比如下图 (1-8) 所示：

“UCR-DTW”和“UCR-ED”模型详解

左图表示从左到右的顺序计算差值，它需要计算 9 个时间步才能判断是否提前结束，而右图找到一个新的计算顺序，这时候只需要计算 5 个时间步就能判断是否提前结束。
那么，现在的问题就变成了，如何找到这些差值之和最大的子序列？
这里有一个疑问，找到的这些子序列是否是连续的？通过阅读源码，可知子序列不一定是连续的。
本文中的做法是，首先对被 Z -normalization 处理过的 Q 序列所有时间点元素的绝对值进行排序，这样做的理论依据是，在进行 DTW 获取序列之间的距离时，qi 可以对应多个序列 C 中的时间点。进行 z -normalization 后的 C 序列服从高斯分布，意味着均值为 0，因此，距离均值 0 最远的 qiqi 对距离值贡献最大，因此对 z -normalizated Q 序列的绝对值进行排序，从而可以快速差值之和最大的子序列。
作者通过实验证明，使用这样的方式找到计算顺序与真实的最好计算顺序的相关度为 0.999。
Reversing the query/data role in LBKeoghLBKeogh
基于 Q，使用 LBKeoghEQLBKeoghEQ 来进行剪枝，这里只需要对 Q 计算一次的 U 和 L，从而可以节省很多时间和空间开销；如果全部采用 LBKeoghECLBKeoghEC 来进行剪枝，基于每一个 C，计算 U 和 L，那么会增加很多的计算成本。
因此，LBKeoghECLBKeoghEC 策略是可选项，只有当 LBKeoghEQLBKeoghEQ 剪枝效果不太理想的时候，可以“just-in-time”的策略来使用 LBKeoghECLBKeoghEC 来辅助 LBKeoghEQLBKeoghEQ 提高剪枝效率，从而大大降低空间开销。对于 LBKeoghECLBKeoghEC 的时间开销，可以通过剪枝来降低完整 DTW 的时间开销来抵消掉。对两种计算策略的直观理解如下图 (1-9) 所示:
“UCR-DTW”和“UCR-ED”模型详解

Use cascading lower bounds
目前存在很多种 lower bound 的计算方式。每一种 lower bound 都可以用于对 DTW 进行剪枝而且时间复杂度可估计。截止目前为止，至少有 18 种 lower bound 机制，作者把它们都实现了一遍，然后分别在 50 个不同的数据集上进行测试和对比，得到的结果如下图 (1-10) 所示：
根据以上实验结果，作者通过级联各种 lower bound 的方式来进行 ED 和 DTW 进行剪枝：
首先，使用时间复杂度为 O(1) 的 LBKimFLLBKimFL，这样可以过滤掉很多的 candidate subsequence，接着，基于 Q，使用 LBKeoghEQLBKeoghEQ 来进行剪枝，
如果，LBKeoghEQLBKeoghEQ 剪枝效果不太理想的时候，使用 LBKeoghECLBKeoghEC 来辅助 LBKeoghEQLBKeoghEQ 提高剪枝效率，
最后，如果以上的剪枝策略全部失效，则依然可以通过 early abandoning 来计算完整 DTW
实验证明，上面使用的每一个 lower bound 策略都能帮助提升 DTW 的速度，去掉任意一个 lower bound 策略都会使得搜索速度加倍。在大规模搜索中，以上的剪枝策略可以节省 99.9999% 的 DTW 算法的时间开销。
实验结果分析

论文中，针对以下这几种实现方式进行性能的比较分析：

Naive: 每个子序列都是从零开始归一化 z 的。每一步都使用完整的欧氏距离或 DTW。(大约有 2 / 3 的文章是基于这种思想来进行相似度计算的)
State-of-the-art: 当前最好的模型基于 Z -normalization,early abandoning 以及使用 lower bound 来辅助完整 DTW 计算这些策略来实现的。(大约有 1 / 3 的文章基于这种思想来进行相似度计算的)
UCR Suite
GOd’s ALgorithm (GOAL) 直接基于均值、方差来进行比较计算相似度的，时间复杂度为 O(1)
GOAL 模型相当于所有解决长度未知无限制的序列搜索问题的最快模型的一个 baseline model。
在用于对比实验的 4 个模型都是使用 UCR Suite 的代码，模型之间的区别只在于把相应的加速代码注释掉而已。

基于随机生成数据集的实验效果对比

由上图可以看到，对于 128 长度的 query，SOFA 和 UCR Suite 算法集之间的性能差异很大。

不同长度 query 的实验对比

接下来，看看对于不同长度 query，这几种模型的性能对比情况：

UCR-DTW python 实现

UCR-DTW 应用了以上所有的优化策略

GitHub:ucr-suite-python

UCR-ED python 实现

UCR-ED 应用的优化策略为：

Early Abandoning of ED
Reorder early abandoning

import time
import math
class UCR_ED(object):
def __init__(self,input_file, query_file,m=128):
self.fp = open(input_file,’r’)
self.qp = open(query_file,’r’)
self.m = m #length of query
self.Q = [None]*self.m#query array
self.T = [0.0](self.m2)#array of current data
self.order = [] #ordering of query by |z(q_i)|
self.bsf = float(‘inf’)
self.loc = 0 #answer:location of the best-so-far match

self.ex,self.ex2,self.mean,self.std=0.0,0.0,0.0,0.0
#用于统计运行时间
self.t1 = time.time()
self.t2 = 0.0

self.Q_normalize()

#Sort the query data
self.sort_query_order()

#read data file, one value at a time
ex = 0.0
ex2 = 0.0
i = 0
while True:
try:
line = self.line_to_float(next(self.fp))
ex += line
ex2 += line*line
self.T[i%m] = line
self.T[(i%m)+m] = line
except:
break
# if there is enough data in T, the ED distance can be calculated
if i>=m-1:
#the current starting location of T
j = (i+1)%m
#Z-norm(T[i]) will be calculated on the fly
mean = ex/self.m
std = ex2/self.m
std = math.sqrt(std-mean*mean)
#Calculate ED distance
dist = self.distance(self.Q, self.T, j, self.m, mean, std, self.order, self.bsf)
if dist<self.bsf:
self.bsf = dist
self.loc = i-m+1
ex -= self.T[j]
ex2 -= self.T[j]*self.T[j]
i+=1

self.fp.close()
self.t2 = time.time()

print(“Location: “, self.loc)
print(“Distance: “,math.sqrt(self.bsf))
print(“Data Scanned: “, i)
print(“Total Execution Time: “,(self.t2-self.t1),’ sec’)

def line_to_float(self, s):
return ConvertELogStrToValue(s.strip())[1]

def sort_query_order(self):
self.Q_tmp = {}
for i in range(self.m):
self.Q_tmp[i] = self.Q[i]
self.Q_tmp = dict(sorted(self.Q_tmp.items(),key=lambda x:x[1]))
#also create another arrays for keeping sorted envelop
self.order = list(self.Q_tmp.keys())

def Q_normalize(self):
i = 0
ex = 0.0
ex2 = 0.0
while i<self.m:
try:
line = self.line_to_float(next(self.qp))
ex += line
ex2 += line*line
self.Q[i] = line
i+=1
except:
break
self.qp.close()

mean = ex/self.m
std = ex2/self.m
std = math.sqrt(std-mean*mean)

#Do z-normalization on query data
for i in range(self.m):
self.Q[i] = (self.Q[i]-mean)/std

def distance(self,Q,T,j,m,mean,std,order,bsf):
”’
Main function for calculating ED distance between the query Q and current data T
Q is already sorted by absolute z-normalization value |Z-normalize(Q[i])|
”’
distance_sum = 0.0
for i in range(m):
if distance_sum>=bsf:
break
x = (T[order[i]+j]-mean)/std
distance_sum += (x-Q[i])*(x-Q[i])

return distance_sum
参考资料

Searching and mining trillions -blog
时间序列的搜索
DTW(Dynamic Time Warping) 动态时间规整
Time Series Classification and Clustering