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

“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 = 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
時間序列的搜索
DTW(Dynamic Time Warping)動態時間規整
Time Series Classification and Clustering