基于Matrix Profile的時間序列變長模體挖掘

朱 旭，朱曉曉，王繼民

(河海大學計算機與信息學院,江蘇 南京 211100)

0 引 言

作為模式發現和相似性搜索的交叉主題，模體挖掘最先由加利福尼亞大學河濱分校的Lin和Keogh等[1]在2002年首次提出，稱這樣的模式為“模體”。在相關文獻中，模體被定義為重復模式、頻繁出現的趨勢或近似和重復的序列、形狀、片段、子序列等。模體挖掘也可為分類的核心[2]、群集[3-4]、異常檢測[5-6]及規則發現[7-8]算法。時間序列模體挖掘技術在遺傳學、金融、工業等諸多領域也得到應用。例如，金融領域的證券交易數據[9]、工業領域的用電數據[10]等。

定長模體挖掘僅挖掘由用戶指定長度的模體。目前對于定長模體挖掘研究較多。文獻[1]提出的Brute-Force算法是第一個模體挖掘算法，算法采取枚舉的方式計算所有子序列之間的距離，然后利用最近鄰算法找到模體，它具有完全的覆蓋率和準確率，被廣泛用作基準算法以衡量其他算法的準確性。Chiu等[11]結合SAX算法與概率思想提出了Random Projection算法。Mueen等[12]提出了MK算法，使用早棄技術對Brute-Force算法進行改進，并提出通過選擇參考序列得到更緊密的下界距離，從而提高算法效率。Yeh等[13]引入STAMP算法，使用快速相似性搜索算法[14]挖掘時間序列模體。Zhu等[15]對STAMP算法進行改進提出了STOMP算法，其計算速度更快，通過快速傅里葉變換實現子序列間的點積來計算子序列之間的距離，通過重用前一個相鄰子序列的點積來加速當前子序列點積的計算。相較于STAMP的隨機順序計算，STOMP是一種有序搜索。

變長模體挖掘無需用戶預先指定模體長度，相比時間序列定長模體挖掘更難以解決，針對此問題的研究也相對較少。變長模體挖掘一般以定長模體挖掘為基礎，通過在可行長度范圍內迭代，得到不同長度的模體。Nunthanid等[16]提出了VLMD算法，利用MK算法計算所有可能長度的模體，然后通過模體分組和提取模體分組代表找到自適應長度的模體。VLMD算法的時間復雜度較高，可擴展性較差。Mueen等[7]提出了MOEN算法，采用下界距離加速Brute-Force算法的效率，以加速枚舉不同長度模體。Tang等[17]擴展了K-Motif算法,基于RP算法，定義搜索空間，并使用分段重疊和分段等價類條件從模體分組中獲取變長模體。

以上算法雖然在不同方面提高了變長模體的發現速度，但仍然存在參數過多、算法可擴展性差以及時間復雜度高等問題。本文提出一種基于Matrix Profile[15]的時間序列變長模體挖掘算法稱為Matrix Profile Variable-Length Motif Discovery (MPVLMD)，基于VLMD算法的計算框架，采用STOMP算法挖掘所有可能長度下的定長模體，使用結合增量距離的下界距離技術提升計算速度，引入模體重疊條件、長度相似性條件和模體分組等價類劃分，該算法能有效地發現變長模體?；赨CR數據集，針對MPVLMD的運行效率和準確性進行實驗，實驗結果表明，本文提出的算法具有較高的計算效率和準確性。

1 相關定義

定義1時間序列[18]。時間序列T是在時間上有序的值序列，可記為T={t1,t2,…,tn}。其中，ti表示一個觀察值，|T|=n表示時間序列T的長度。

定義2時間序列子序列[19]。時間序列子序列是T中一個更小的時間序列，記為Ti,m={ti,ti+1,…,ti+m-1}。其中，i是子序列在T中的開始位置，m是子序列的長度，并且m>1。

定義3時間序列模體[20]。時間序列模體Mw是指時間序列T中，長度為w且彼此相似度最高的一對子序列?？蓪⑵涠x為一個四元組：Mw=(d,L1,L2,w),其中，L1和L2為2個子序列的起始位置，滿足L1

定義4平凡匹配。給定時間序列T，假設存在子序列C和F，其起始位置分別為p和q。如果p=q，那么C和F屬于平凡匹配。

定義5模體重疊[16]。2個長度分別為i和j的模體Mi和Mj，當且僅當它們滿足條件：Mi.L1≤Mj.L1

定義6模體分組[16]。模體分組MG是模體的集合，且包含的2個連續模體互相重疊。

定義7模體分組代表[16]。模體分組代表MR，指模體分組MG所包含的所有不同長度的模體中d值最小的模體。

定義8模體分組重疊。2個模體分組MGi和MGj，假設任意Mw∈MGi，Mw∈MGj當且僅當它們滿足條件Mw.L1=Mx.L1‖Mw.L1=Mx.L2‖Mw.L2=Mx.L1‖Mw.L1=Mx.L2時，稱模體分組MGi和MGj重疊。

定義9變長模體集合。給定時間序列T和模體長度范圍[lmin,…,lmax]，變長模體集合是所有滿足模體長度約束，且彼此不重疊的模體的集合。

定義10距離矩陣。距離矩陣是時間序列中所有的子序列之間歐氏距離的矩陣。距離矩陣的特征表示為：

其中，di,j為時間序列中第i個子序列與第j個子序列之間的歐氏距離。

定義11Matrix Profile[15]。Matrix Profile是時間序列中所有子序列與其最相似的子序列間距離的向量。即距離矩陣中每一列的最小值。向量表示為：(P1,P2,…,Pn-m+1),其中，n表示時間序列長度，m表示子序列長度，Pi表示第i個子序列與所有子序列間距離的最小值Min(Di)。

定義12Matrix Profile Index[15]。Matrix Profile Index是時間序列中所有子序列與其最相似的子序列的索引向量。向量索引表示為：(I1,I2,…,In-m+1),其中Ii表示第i個子序列的最近的鄰居。

2 STOMP算法

2016年，Yeh等提出了基于Matrix Profile的時間序列模體挖掘算法STAMP[13]，在其基礎上，Zhu等進一步提出了STOMP算法[15]，在距離矩陣基礎上提出了新的數據結構Matrix Profile，該向量中的最小距離值對應的是彼此相似度最高的2個子序列，即為該時間序列的模體。STOMP算法通過快速傅里葉變換計算距離矩陣，從而可以有效地計算給定時間序列的Matrix Profile和Matrix Profile Index[15]。算法包含提取子序列、子序列全連接和模體發現這3個步驟。

2.1 提取子序列

根據定義2，使用長度為m的滑動窗口對時間序列T={t1,t2,…,tn}進行分段，截取所有長度為m的子序列，T1,m={t1,t2,…,tm},T2,m={t2,t3,…,tm+1},…,Tn-m+1,m={tn-m+1,tn-m+2,…,tn}，獲取所有子序列Ti,m={ti,ti+1,…,ti+m-1}(1≤i≤n-m+1)。

2.2 子序列全連接

子序列全連接用來計算時間序列的距離矩陣。計算所有子序列的均值μ和標準差σ，通過快速傅里葉變換計算子序列Ti,m與所有子序列之間的點積，使用z歸一化歐氏距離公式(1)計算2個子序列間的距離。

(1)

其中,m是子序列的長度,μi和μj分別是Ti,m和Tj,m的均值,σi和σj分別是Ti,m和Tj,m的標準差，QTi,j是Ti,m和Tj,m的點積。計算di,j所需的時間只取決于QTi,j計算所需的時間，QTi-1,j-1可以分解為：

(2)

并且QTi,j可以分解為：

(3)

由式(2)和式(3)可獲得式(4)：

QTi,j=QTi-1,j-1-ti-1tj-1+ti+m-1tj+m-1

(4)

當計算第i個子序列與T中所有子序列的點積時，可以利用第i-1個子序列與T中所有子序列的點積結果加速計算。

求得第i個子序列與T中所有子序列點積之后，使用z歸一化歐氏距離計算子序列之間的距離。

2.3 模體發現

根據定義10，使用Matrix Profile存儲距離矩陣中每列元素的最小值。同時根據定義11，將最小值的行號作為索引存入Matrix Profile Index。提取Matrix Profile中的最小值對應的索引所指向的子序列是時間序列T的模體。其中，Matrix Profile的結構如圖1所示。

圖1 Matrix Profile結構圖

3 基于Matrix Profile的時間序列變長模體挖掘算法MPVLMD

本文算法基于VLMD算法的計算框架，首先，使用STOMP算法作為子程序，并引入結合增量計算的下界距離加速策略，提升模體提取速度；其次，加入模體重疊和長度相似性條件進行模體分組，踢除過短和存在平凡匹配的模體；再次，加入模體分組重疊條件對模體分組進行等價類劃分，剔除提取模體代表時存在的過長模體；最后，提取每個分組等價類中的模體代表，模體代表集合即為變長模體?；贛atrix Profile的時間序列變長模體挖掘算法的流程如圖2所示。

圖2 MPVLMD算法流程圖

3.1 模體提取

STOMP算法只能提取固定長度的模體，本文算法以STOMP作為子程序，在所有可能長度范圍內迭代該程序，結合增量計算的下界距離加速策略，加速模體提取。

3.1.1 求解z值數組

依據下界距離的計算流程，首先需要計算對應模體長度的z值。在提取所有可能長度模體的時候，模體長度范圍默認為從2～n/2(其中n為時間序列長度)，則需要計算長度為n/2-1的z值數組，其具體計算流程如算法1所示。

算法1求解z值數組

forj←m0+1 tomxdo

fori←1 ton-jdo

Maxj=Y

return Max

Output:z值數組Max//用于計算下界距離

3.1.2 模體提取

1)提取模長為m0的模體。

使用STOMP算法提取模長為m0的模體，生成按距離升序排列的彼此最相似的子序列對列表List。

2)結合增量計算的下界距離加速策略提取模長為m0+1的模體。

使用緩存技術[21]保存時間序列值的累積和與累積平方和，便于下界距離和增量距離計算。

3)提取所有可能長度的模體。

重復上述步驟，直到獲取所有可能長度的模體。

下面闡述下界距離和增量距離的具體求解過程。

通過式(5)求模體長度為m0+1的下界距離：

(5)

其中，z=maxi(ti+m0-μi,m0)/σi,m0，(1≤i≤n-m0，當模長為m0時，時間序列中第i+m0個元素歸一化后的最大值)，d表示模長為m0的子序列間的z歸一化歐氏距離集合List中的最大值。

(6)

基于長度為m的tx、ty、(tx)2、(ty)2和txty的運行和，就可以增量計算長度為m+1的距離函數的值。

算法2模體提取

Input:時間序列T,z值數組Max

List←STOMP(T,m0)//提取長度為m0的模體，以及n-m0+1對相似子序列集合List

z←Maxm0+1

forj←m0+1 tomxdo

forp←1 to List.Count do

i←List(P).i;k←List(P).k;

distance(Ti,j,Tk,j)←List(P).d+

Sum(Ti,j-1,Tk,j-1,(Ti,j-1)2,(Tk,j-1)2,(Ti,j-1,Tk,j-1))；

d←distance(Ti,j,Tk,j);

NewList.Add (d,i,k)//將距離d,位置i、k存入NewList中

best←NewList.Min

z←Maxj

L1j←i,L2j←k

output List.Min for lengthj

else

List←STOMP(T,j)

Output:模體候選集合List

3.2 模體分組

模體提取階段產生的模體候選集合List中可能存在較長模體覆蓋較短模體的問題，因此對List中的模體進行分組。

遍歷List集合中的模體，按照定義5，將滿足模體重疊條件的2個模體置入相同模體分組中，反之創建新的模體分組，并將其中未分組的一個模體作為首個元素存儲到其中。

遍歷所有的模體分組，并對同一個分組中的模體，使用長度相似性條件HMw=?n/2w」，修剪過短模體。如果模體Mw的HMw值與其他模體的HMother值不同，剔除模體Mw。其中，n表示時間序列長度，w表示模體長度。

算法3模體分組

Input:模體候選集合List

for eachMi∈List do

for eachMj∈List，j>i

ifMioverlaps withMj

Mj.groupid=Mi.groupid

for each groupID groupid do

Put all motifs that have the groupidiinto groupi

for eachMw∈groupido

if HMwnot equals to HMother

deleteMw

return group

Output:模體分組集合group

3.3 模體分組等價類劃分

完成模體分組后，模體分組中的模體可能存在以下情況：(si,ti)∈groupi，(sj,tj)∈groupj，并且ti=tj。由于si與ti相似，sj與tj相似，則si與sj有很大概率是相似的。對模體分組集合group進行模體分組重疊判斷，遍歷group中的每一個模體分組，將滿足定義8的模體分組置入同一個模體分組等價類中。

算法4模體分組等價類劃分

Input:模體分組集合group

for each groupi∈group do

for each groupj∈group，j>i

if ?Mx∈groupioverlaps with?My∈groupj

//如果Mx和My滿足模體分組重疊條件

groupj.equalclassid=groupi.equalclassid

return equalclass

Output:模體分組等價類集合equalclass

3.4 變長模體提取

完成模體分組等價類劃分后，每個模體分組等價類將包含平凡匹配的模體。

1)提取等價類模體代表。

每個模體分組等價類中有多個模體分組，每個模體分組中有多個模體。提取模體分組等價類中每個模體分組中z歸一化歐氏距離最小的模體，作為每個模體分組的模體代表，每個模體等價類將包含多條模體代表。將這些模體分組的模體代表按照z歸一化歐氏距離正序排列，選擇中間位置模體代表的z歸一化歐氏距離(如果模體分組個數為奇數即為中間位置模體代表的z歸一化歐氏距離，如果是偶數取中間2個模體代表的z歸一化歐氏距離的均值)作為距離最大值，刪除z歸一化歐氏距離大于該最大距離的模體代表(修剪具有較大z歸一化歐氏距離的過長模體)。

2)輸出變長模體。

經過模體分組，模體等價類劃分等操作剔除過長、過短和平凡匹配模體后，輸出每個模體分組等價類中長度最長的模體代表的集合，所輸出的所有不同長度模體代表的集合即為時間序列的變長模體。

算法5變長模體提取

Input:模體分組等價類集合equalclass

for each equalclassid∈equalclass do

for each groupi∈equalclassid do

for eachMw∈groupido//提取每個分組中每個等價類中的模體代表

//按照z歸一化歐氏距離正序排列

//提取中間位置的模體作為模體分組等價類的模體代表

return VMotifList

Output:變長模體集合VMotifList

4 實驗與結果分析

4.1 實驗數據

以UCR[22]的部分數據集作為實驗數據，數據集信息如表1所示。

表1 數據集中所有植入模式的詳細信息

UCR數據集是由事先確定好模式長度的已知模式組成，將UCR數據集中已知模式隨機植入到隨機游走數據中，創建實驗所用數據集Dataset。

4.2 實驗方法

本文從模體發現的時間效率和準確率這2個方面來分析算法。

實驗1驗證本文算法的準確性。將本文提出的方法與MN方法[23]以及原始VLMD算法進行比較。采用準確性檢測方法AoD[24]作為準確性評價指標，該方法最初用于計算子序列匹配問題中任意子序列對的重疊百分比。本實驗用其計算算法輸出模體與植入模體間的重疊比，以衡量各算法的準確性。

(7)

其中，

max(Lx,Ly)+1

(8)

min(LX,Ly)+1

(9)

其中，Lx、Ly為模體開始位置，Sx、Sy為模體長度。

實驗2驗證不同算法的時間效率。增加算法中待挖掘的時間序列模體長度，對比VLMD的子程序MK算法和MPVLMD的子程序STOMP算法,獲取不同長度模體所需的時間。

4.3 實驗結果與分析

1)準確性分析。

為了驗證MPVLMD算法的準確性，基于2個不同的數據集，分別使用本文算法、MN算法和VLMD算法進行多次模體挖掘實驗，統計各算法挖掘模體的結果。并基于本文選用的準確性衡量方法AoD，計算各算法輸出模體與預先植入模體的重疊比。所得計算結果如表2所示，具體展示見圖3。

表2 各數據集中不同算法發現植入模體的準確率

圖3 不同算法發現植入模體的準確率

由圖3可知，MPVLMD算法能夠發現所有的植入模體，其發現模體的準確率要優于VLMD算法。同時，針對Dataset中的Beef數據，VLMD算法會出現不能發現植入模體的情況，這是因為模體候選集中存在過長、過短和平凡匹配模體影響算法發現模體的整體準確性，這也表明本文算法引入的長度相似性和模體分組等價類條件在篩除無效模體上的積極作用。雖然MN算法也能有效地發現所有的植入模體，但是其發現模體的準確率整體來看要略低于MPVLMD算法，表明本文算法不僅能夠有效地發現不同長度的模體，而且具有較高的準確率。

2)時間效率分析。

為了驗證MPVLMD算法所選用的模體發現算法STOMP較MK算法的優勢，基于Dataset數據集，取算法中待挖掘的時間序列模體長度分別為64、128、256、512、1024進行實驗，記錄STOMP算法和MK算法獲取不同長度模體所需的平均時間，如圖4所示。

圖4 模體發現耗時

由圖4可知，在模體長度相同的情況下，STOMP算法的效率遠優于MK算法。并且隨著模體長度和數據集長度的增加，MK算法運行所需時間呈現指數增長趨勢，其效率快速降低，相反STOMP算法始終能夠保持較高的運行速度，其運行時間對模體和數據集長度變化不敏感。

綜合上述2個實驗的結果可知，MPVLMD算法能夠有效地發現時間序列中不同長度的模體，在準確率和時間效率方面均優于原始的VLMD算法。

5 結束語

本文提出了一種基于Matrix Profile的變長時間序列模體挖掘算法MPVLMD，該算法從所有可能長度的模體中自適應返回合適長度的模體。通過模體提取、模體分組、模體分組等價類劃分和變長模體提取這4個步驟，將大量可能的滑動窗口長度縮減到真正有意義的幾個模體長度，能夠有效地發現時間序列中不同長度的模體。

基于UCR數據集的實驗結果表明，本文算法具有較高的準確率和時間性能。