摘要：這兩種順序的計(jì)算成本對(duì)比如下圖所示和模型詳解左圖表示從左到右的順序計(jì)算差值，它需要計(jì)算個(gè)時(shí)間步才能判斷是否提前結(jié)束，而右圖找到一個(gè)新的計(jì)算順序，這時(shí)候只需要計(jì)算個(gè)時(shí)間步就能判斷是否提前結(jié)束。

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

DTW(Dynamic Time Warping)及已知的優(yōu)化策略

計(jì)算兩個(gè)時(shí)間序列 Q和C之間的相似度，常用的度量方法是歐式距離(ED)，計(jì)算公式如下圖(1-1)所示：

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

從上圖可以看到歐式距離的局限性：歐式距離通過(guò)建立兩個(gè)序列之間的一一對(duì)應(yīng)關(guān)系，使得Q和C之間的波峰沒(méi)有對(duì)齊，因此計(jì)算得到的序列相似度存在較大的偏差。DTW 算法則可以很好的解決這個(gè)問(wèn)題。

大部分情況下，兩個(gè)序列整體上具有非常相似的形狀，但是這些形狀在時(shí)間軸上并不是對(duì)齊的。所以在計(jì)算兩個(gè)序列的相似度之前，需要將其中一個(gè)(或者兩個(gè))序列在時(shí)間軸上進(jìn)行warping，使得兩個(gè)序列波峰更好的對(duì)齊。

DTW就是實(shí)現(xiàn)這種warping的一種有效方法。換句話說(shuō)，DTW算法通過(guò)找到兩個(gè)序列之間的另外一種非一一對(duì)應(yīng)的映射關(guān)系，這個(gè)映射關(guān)系也稱(chēng)為warping path。以上面的Q,C為例，得到的對(duì)應(yīng)關(guān)系如下圖(1-2)中的灰色線段所示:

DTW算法的求解過(guò)程的直觀理解為構(gòu)建一個(gè)nxn的矩陣(此處假設(shè)Q和C的時(shí)間序列長(zhǎng)度都為n，矩陣中的元素(i,j)表示序列Q的時(shí)間點(diǎn)qi和序列C的時(shí)間點(diǎn)cj之間的歐式距離）目標(biāo)是在矩陣中找到一條從(0,0)到(n,n)的路徑，使得路徑上的所有元素之和最小。

如下圖(1-3)所示，紅色標(biāo)識(shí)的路徑即為warping path：

目前已知的DTW順序搜索的四種優(yōu)化方法如下:

去掉平方根計(jì)算
通過(guò)使用lower bounds來(lái)剪枝，因?yàn)閘ower bounds的計(jì)算時(shí)間復(fù)雜度都小于DTW的時(shí)間復(fù)雜度。比如：
LBKimFLLBKimFL的時(shí)間復(fù)雜度為O(1)
本文在實(shí)現(xiàn)時(shí)，由于對(duì)時(shí)間序列進(jìn)行標(biāo)準(zhǔn)化后，時(shí)間序列數(shù)據(jù)中的最大和最小值對(duì)于整個(gè)lower bound距離貢獻(xiàn)較小，因此，去掉原來(lái)LB_kim算法(時(shí)間復(fù)雜度為O(n))中提取的四個(gè)特征點(diǎn)中的最大值和最小值，使得時(shí)間復(fù)雜度降為O(1)。但是，實(shí)現(xiàn)中為了使此策略發(fā)揮最大作用，作者還提取了第2,3和倒數(shù)第2,3個(gè)時(shí)間點(diǎn)，來(lái)進(jìn)行級(jí)聯(lián)剪枝。(詳情參考算法實(shí)現(xiàn)lb_kim_hierarchy方法)

基于這條優(yōu)化策略，作者提出的Reorder early abandoning可以進(jìn)一步降低計(jì)算成本。

即在計(jì)算ED或者LBKeoghLBKeogh的時(shí)候，如果當(dāng)前的兩個(gè)序列的時(shí)間點(diǎn)(1,k)(注：k<=|Q|)之間的差值平方之和，大于當(dāng)前兩個(gè)序列最小距離值best-so-far，那么可以提前結(jié)束Q和C是否相似的判斷。計(jì)算過(guò)程如下圖(1-4)所示：

4.Early Abandoning of DTW

計(jì)算完整的LBKeoghLBKeogh的值，仍然需要計(jì)算完整的DTW的值，我們可以通過(guò)利用部分的LBKeoghLBKeogh 的值來(lái)減少DTW的計(jì)算量。

比如，先從左到右計(jì)算時(shí)間點(diǎn)[1,k]的DTW值，然后在時(shí)間點(diǎn)[k+1,n]，復(fù)用前面計(jì)算好的LBKeoghLBKeogh 值。最終得到的距離值依然為完整的DTW值的lower bound。這樣的話，我們就可以使用stop early策略，每計(jì)算當(dāng)前時(shí)間點(diǎn)的DTW值，就可以復(fù)用前面計(jì)算好的LBKeoghLBKeogh來(lái)獲得整個(gè)子序列的lower bound值。

通過(guò)比較這個(gè)lower bound和當(dāng)前最小距離值best-so-far進(jìn)行比較，如果當(dāng)前的lower bound值大于best-so-far，那么可以提前結(jié)束DTW的計(jì)算。

這種方式的直觀表示如下圖(1-5)：

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

以上，介紹完了前人對(duì)DTW的四種優(yōu)化方案。

還有一種已知的提升DTW的計(jì)算速度的策略即使用多核計(jì)算資源。

UCR Suite 的優(yōu)化策略

相關(guān)概念及定義

定義1:

時(shí)間序列T是一個(gè)有序列表：T=t1,t2,?,tm。然而源數(shù)據(jù)是一個(gè)很長(zhǎng)的時(shí)間序列，我們最終需要把它與一個(gè)更短的子序列進(jìn)行相似度比較。

定義2:

子序列Ti,k是時(shí)間序列T中的一個(gè)子序列，它起始于ti，長(zhǎng)度為k，即:

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

這里，我們把Ti,kTi,k記為C，作為與query Q比較相似的候選子序列。令Q的長(zhǎng)度為|Q|=n。

定義3:

Q和C之間的歐式距離(|Q|=|C|)定義為 (公式1)：

路徑P的第t個(gè)元素定義為pt=(i,j)tpt=(i,j)t，則我們可以把warping path表示為 (公式2)：

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

優(yōu)化策略:

1.Early Abandoning Z-normalization

在計(jì)算DTW距離之前都需要對(duì)Q和C進(jìn)行標(biāo)準(zhǔn)化處理，但是對(duì)整個(gè)數(shù)據(jù)集進(jìn)行標(biāo)準(zhǔn)化的復(fù)雜度太高了，因此這里使用online Z-normalization，這樣的話就可以采用early stop的策略來(lái)提前結(jié)束normalization的計(jì)算。

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

當(dāng)使用online Z-normalization的時(shí)候，當(dāng)前遍歷到源序列T中的第k個(gè)時(shí)間點(diǎn)，所計(jì)算得到的時(shí)間點(diǎn)元素累加和以及時(shí)間點(diǎn)元素的平方累加和表示為(公式4)：

那么對(duì)于k-m+1 到k之間的這m個(gè)時(shí)間點(diǎn)對(duì)應(yīng)的均值和方差的計(jì)算公式如下所示(公式5)

所以，基于online Z-normalization的abandon normalization early策略偽代碼如下圖(1-7)所示：

論文中，作者提到此處存在浮點(diǎn)計(jì)算誤差累加的問(wèn)題，這里通過(guò)每比對(duì)完100萬(wàn)子序列，就進(jìn)行一次完整的Z-normalization，進(jìn)而消除誤差累加的問(wèn)題。

Reorder early abandoning
前面early abandoning策略的計(jì)算方式都是從子序列的第一個(gè)時(shí)間點(diǎn)開(kāi)始，自左向右進(jìn)行計(jì)算的。本文提出一種策略是先快速找到Q和C之間差值之和最大的子序列，然后根據(jù)它來(lái)判斷這個(gè)子序列是否大于best-so-far值，從而可達(dá)到降低計(jì)算成本的目的。這兩種順序的計(jì)算成本對(duì)比如下圖(1-8)所示：

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

左圖表示從左到右的順序計(jì)算差值，它需要計(jì)算9個(gè)時(shí)間步才能判斷是否提前結(jié)束，而右圖找到一個(gè)新的計(jì)算順序，這時(shí)候只需要計(jì)算5個(gè)時(shí)間步就能判斷是否提前結(jié)束。
那么，現(xiàn)在的問(wèn)題就變成了，如何找到這些差值之和最大的子序列？
這里有一個(gè)疑問(wèn)，找到的這些子序列是否是連續(xù)的？通過(guò)閱讀源碼，可知子序列不一定是連續(xù)的。
本文中的做法是，首先對(duì)被Z-normalization 處理過(guò)的Q序列所有時(shí)間點(diǎn)元素的絕對(duì)值進(jìn)行排序，這樣做的理論依據(jù)是，在進(jìn)行DTW獲取序列之間的距離時(shí)，qi可以對(duì)應(yīng)多個(gè)序列C中的時(shí)間點(diǎn)。進(jìn)行z-normalization后的C序列服從高斯分布，意味著均值為0，因此，距離均值0最遠(yuǎn)的qiqi對(duì)距離值貢獻(xiàn)最大，因此對(duì)z-normalizated Q序列的絕對(duì)值進(jìn)行排序，從而可以快速差值之和最大的子序列。
作者通過(guò)實(shí)驗(yàn)證明，使用這樣的方式找到計(jì)算順序與真實(shí)的最好計(jì)算順序的相關(guān)度為0.999。
Reversing the query/data role in LBKeoghLBKeogh
基于Q，使用LBKeoghEQLBKeoghEQ來(lái)進(jìn)行剪枝，這里只需要對(duì)Q計(jì)算一次的U和L，從而可以節(jié)省很多時(shí)間和空間開(kāi)銷(xiāo)；如果全部采用LBKeoghECLBKeoghEC來(lái)進(jìn)行剪枝，基于每一個(gè)C，計(jì)算U和L，那么會(huì)增加很多的計(jì)算成本。
因此，LBKeoghECLBKeoghEC策略是可選項(xiàng)，只有當(dāng)LBKeoghEQLBKeoghEQ剪枝效果不太理想的時(shí)候，可以“just-in-time” 的策略來(lái)使用LBKeoghECLBKeoghEC來(lái)輔助LBKeoghEQLBKeoghEQ提高剪枝效率，從而大大降低空間開(kāi)銷(xiāo)。對(duì)于LBKeoghECLBKeoghEC的時(shí)間開(kāi)銷(xiāo)，可以通過(guò)剪枝來(lái)降低完整DTW的時(shí)間開(kāi)銷(xiāo)來(lái)抵消掉。對(duì)兩種計(jì)算策略的直觀理解如下圖(1-9)所示:
“UCR-DTW”和“UCR-ED”模型詳解

Use cascading lower bounds
目前存在很多種lower bound的計(jì)算方式。每一種lower bound都可以用于對(duì)DTW進(jìn)行剪枝而且時(shí)間復(fù)雜度可估計(jì)。截止目前為止，至少有18種lower bound機(jī)制，作者把它們都實(shí)現(xiàn)了一遍，然后分別在50個(gè)不同的數(shù)據(jù)集上進(jìn)行測(cè)試和對(duì)比，得到的結(jié)果如下圖(1-10)所示：
根據(jù)以上實(shí)驗(yàn)結(jié)果，作者通過(guò)級(jí)聯(lián)各種lower bound的方式來(lái)進(jìn)行ED和DTW進(jìn)行剪枝：
首先，使用時(shí)間復(fù)雜度為O(1) 的LBKimFLLBKimFL，這樣可以過(guò)濾掉很多的candidate subsequence，接著，基于Q，使用LBKeoghEQLBKeoghEQ來(lái)進(jìn)行剪枝，
如果，LBKeoghEQLBKeoghEQ剪枝效果不太理想的時(shí)候，使用LBKeoghECLBKeoghEC來(lái)輔助LBKeoghEQLBKeoghEQ提高剪枝效率，
最后，如果以上的剪枝策略全部失效，則依然可以通過(guò)early abandoning 來(lái)計(jì)算完整DTW
實(shí)驗(yàn)證明，上面使用的每一個(gè)lower bound策略都能幫助提升DTW的速度，去掉任意一個(gè)lower bound策略都會(huì)使得搜索速度加倍。在大規(guī)模搜索中，以上的剪枝策略可以節(jié)省99.9999%的DTW算法的時(shí)間開(kāi)銷(xiāo)。
實(shí)驗(yàn)結(jié)果分析

論文中，針對(duì)以下這幾種實(shí)現(xiàn)方式進(jìn)行性能的比較分析：

Naive:每個(gè)子序列都是從零開(kāi)始?xì)w一化z的。每一步都使用完整的歐氏距離或DTW。(大約有2/3的文章是基于這種思想來(lái)進(jìn)行相似度計(jì)算的)
State-of-the-art: 當(dāng)前最好的模型基于Z-normalization,early abandoning以及使用lower bound來(lái)輔助完整DTW計(jì)算這些策略來(lái)實(shí)現(xiàn)的。(大約有1/3的文章基于這種思想來(lái)進(jìn)行相似度計(jì)算的)
UCR Suite
GOd’s ALgorithm (GOAL) 直接基于均值、方差來(lái)進(jìn)行比較計(jì)算相似度的，時(shí)間復(fù)雜度為O(1)
GOAL模型相當(dāng)于所有解決長(zhǎng)度未知無(wú)限制的序列搜索問(wèn)題的最快模型的一個(gè)baseline model。
在用于對(duì)比實(shí)驗(yàn)的4個(gè)模型都是使用UCR Suite的代碼，模型之間的區(qū)別只在于把相應(yīng)的加速代碼注釋掉而已。

基于隨機(jī)生成數(shù)據(jù)集的實(shí)驗(yàn)效果對(duì)比

由上圖可以看到，對(duì)于128長(zhǎng)度的query，SOFA和UCR Suite算法集之間的性能差異很大。

不同長(zhǎng)度query的實(shí)驗(yàn)對(duì)比

接下來(lái)，看看對(duì)于不同長(zhǎng)度query，這幾種模型的性能對(duì)比情況：

UCR-DTW python實(shí)現(xiàn)

UCR-DTW應(yīng)用了以上所有的優(yōu)化策略

GitHub:ucr-suite-python

UCR-ED python 實(shí)現(xiàn)

UCR-ED 應(yīng)用的優(yōu)化策略為：

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
#用于統(tǒng)計(jì)運(yùn)行時(shí)間
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 = 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 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
時(shí)間序列的搜索
DTW(Dynamic Time Warping)動(dòng)態(tài)時(shí)間規(guī)整
Time Series Classification and Clustering