基于張量多線性PCA的多變量時間序列模式匹配

董紅玉，陳曉云，潘江山

(福州大學數學與計算機科學學院，福建福州 350116)

多變量時間序列(multivariate time series，MTS)是指多個變量依照時間先后順序得到的一系列觀察值的集合，常用一個n×m的矩陣表示，其中n代表序列長度，m代表變量個數［1，2］．MTS數據廣泛存在于視頻監控、金融和醫療等各種應用領域，故對其挖掘具有廣闊的應用前景［3］．依照挖掘目標的不同，MTS挖掘包括模式匹配［3－4］、主題發現［4］、異常檢測［5］、監督和無監督的分類［1］等．

MTS模式匹配是指將數據集的一個MTS樣本的模式表示作為輸入，計算與數據集中的其余模式相似度的過程［4］．MTS的模式表示和相似性度量是MTS模式匹配的兩大關鍵步驟．

1 相關研究

現有MTS模式匹配方法可分為參數法和非參數法．參數法是指對MTS建立模型，用模型參數度量其相似度，如文［6］對MTS建立向量自回歸模型，計算各模型參數間的相似度．而難以選擇恰當的模型是此類方法的不足．

非參數法則先提取MTS的模式表示，之后根據選取的相似性度量計算模式之間的匹配程度．文［7］提出基于點分布特征的模式匹配方法，提取MTS的局部極值點作為其模式，并由局部極值點的統計分布特征組成特征向量，之后計算特征向量的歐氏距離，此方法適用于小規模的MTS數據集，對大規模數據集，其匹配效果較差．文［8］采用多維分段線性擬合方法提取擬合線段的傾斜角、時間跨度與序列長度的比值作為多變量時間序列的模式表示，并用DTW距離做為模式相似性度量．其不足之處在于模式參數由經驗給定，缺乏自適應性且計算DTW距離的時間開銷較大．文［9］提出基于單線性PCA的模式匹配算法，利用PCA對MTS的序列長度進行降維，并提取降維后的前z個主成分作為原MTS的模式表示，而模式相似性度量則采用如下的夾角余弦平方公式:

其中:X和Y是兩個MTS樣本，θij表示X的第i個主成分與Y的第j個主成分的夾角，主成分的個數z由其貢獻率決定．該方法的不足之處在于利用PCA降維時僅考慮序列長度降維而忽略了多變量間的相關性，且夾角余弦平方不滿足三角不等式;文［10］將MTS的協方差矩陣的右奇異分解矩陣作為其模式表示，并利用奇異值對右奇異矩陣加權，以加權的右奇異矩陣列向量間的內積絕對值作為模式相似性度量，其公式如下:

2 基于張量多線性主成分的MTS模式匹配

文［9］采用的單線性PCA方法和文［10］采用的奇異值分解方法僅對MTS序列長度進行降維，MTS的多個變量間也可能存在冗余．文本將多變量時間序列看成一個二階張量(矩陣)，通過對MTS進行多線性PCA變換，變換后得到的序列模式不僅序列長度降低，變量個數也可減少，即對MTS進行雙向降維，降維后的MTS仍以矩陣的形式表示．為了度量矩陣形式表示的模式間相似度，本文采用Frobenius范數作為MTS模式的相似性度量．下面給出方法的相關定義與原理．

多變量時間序列的數學表達．多變量時間序列(multivariate time series，MTS)是指由多個變量按照時間順序所記錄的一系列觀察值的集合，記為Xi(t)，(i=1，2，…，m;t=1，2，…，n)．其中i表示變量個數(m≥2)，t表示序列長度．它可以用一個n×m的矩陣表示如下:其中:xit表示第i個變量在t時刻的觀察值．

圖1為一個來自JV數據集的MTS樣本，其變量個數為12，序列長度為21．

張量［11］是一個多維數組．用數學觀點來說，一個N階張量是N向量空間中一個實多重線性函數．張量Α∈RI1×I2×…×IN的階為 N．

圖1 MTS的3D圖像Fig．1 Tree－ dimensional graphics of MTS

2．1 基于張量多線性PCA的MTS模式表示

MTS樣本X∈Rn×m是張量空間Rn?Rm的一個二階張量．設(u1，u2，…，un)，(v1，v2，…，vm)分別是空間Rn，Rm的標準正交集，則MTS樣本X可唯一寫成

其中:{uivTj}構成張量空間Rn?Rm的基．定義二階張量U=(u1，u2，…，ul)∈Rn×l和V=(v1，v2，…，vs)∈Rm×s(l＜n，s＜m)．假設μ，ν分別是空間中由基向量{ui}li=1，{vj}sj=1張成的子空間，則二階張量X∈Rn×m在低維子空間μ?ν上的重構可寫為

張量多線性PCA的目標［12］是尋找少數幾個最大特征值對應特征向量所張成的張量子空間，使重構誤差最?。畬τ贛TS數據集D={X1，X2，…，Xd}，其張量多線性PCA重構的目標函數為［13］

其中:Yi是MTS樣本Xi在張量子空間μ?ν上的重構，為所有MTS樣本在子空間μ?ν的重構均值

由矩陣Frobenius范數的定義 Q2=tr(QQT)=tr(QTQ)，可推出［12］

其中:

則最佳二階投影張量U，V分別是DV，DU的最大特征值對應特征向量所組成的二階張量．不失一般性，取 U，V 為正交矩陣，即 UUT=I，VVT=I，故［12］

從而最佳二階投影張量U，V分別是D'V，D'U的最大特征值對應特征向量所組成的二階張量．

由式(5)可知，MTS樣本在低維子空間的重構由二階投影張量U，V確定．因此，二階投影張量U，V的維數選擇至關重要．分別求得D'V，D'U的特征值和特征向量，本文選擇特征值貢獻率占所有特征值98%的特征值所對應的特征向量構成U，V的列向量．

2．2 MTS的相似性度量

提取MTS的低維重構作為其模式表示之后，需要一個恰當的相似性度量去描述兩兩模式間的相似度以完成模式匹配．對于兩個MTS樣本X1=(x1ij)n×m和 X2=(x2ij)n×m，經多線性PCA重構后的模式表示為Y1=(y1ij)l×s和 Y2=(y2ij)l×s(l＜n，s＜m)，樣本X1和X2的相似性度量定義如下:

其中: ·F為Frobenius范數．由于Frobenius范數滿足三角不等式，故SDLPC也滿足三角不等式．

3 實驗結果與分析

3．1 實驗數據

本文選取四組 MTS 數據集進行實驗，分別為 Wafer［9－10］、ECG［9－10］、JV［14］和 LP1［15］，表 1 概要描述了這四組實驗數據．

表1 實驗數據概述Tab．1 The Summary of datasets used in experiments

Wafer是硅晶體生產過程中用6個真空傳感器所記錄的半導體微電子序列，每個硅晶體用一個MTS刻畫，分為normal和abnormal兩類，327個樣本中200個屬于normal類，127個屬于abnormal類;ECG是用2個電極測得的一組心跳數據，包括normal和abnormal兩類，200個樣本中133個屬于normal類，67個屬于abnormal類;JV是記錄9個測試者的12個同態譜數據描述的日文元音發音過程，每次發音作為一個MTS;LP1是Robot Execution Failure數據集的一個子數據集，它通過6個傳感器對Robot進行異常監控．

3．2 實驗方法

采用k－近鄰法進行實驗．假定數據集包含d個樣本，任意選取一個樣本X作為輸入，對數據集中所有樣本進行重構，得到其模式表示，采用k－近鄰法尋找與輸入樣本X最相似的k個樣本，本文分別取k=1，5，10，最后統計這k個樣本中與輸入樣本X具有相同類標簽的樣本個數n，由式(15)計算準確率c

依次將其它樣本作為輸入，重復以上過程，得到d個準確率．根據本文k的取值，準確率c的取值應為{0，0．1，0．2，…，1}，記為 ci(i=1，2，…，11)，計算準確率的比率 p如下:

進一步由式(17)計算整個數據集的準確率期望值

準確率期望值越高，表明該模式匹配方法越好．此外，由于在進行樣本重構時，張量多線性PCA要求樣本序列長度一致．因此，對于不等長數據集，通過截取使所有樣本序列長度與數據集中最短的序列長度一致．

3．3 三種模式匹配方法準確率比較

實驗主要比較三種多變量時間序列模式匹配方法的準確率:本文方法(張量多線性PCA+Frobenius范數)、文［9］方法(單線性PCA+夾角余弦平方和)及文［10］方法(奇異值分解+加權內積絕對值)的模式匹配準確率．分別采用這三種模式匹配方法在Wafer、ECG、LP1和JV四組數據集上測試，統計近鄰數k=1，5，10時三種方法各準確率上的匹配樣本數．下文給出ECG數據集的實驗結果，如表2所示，其余數據集的實驗結果詳見表3．

表2 三種模式匹配方法在ECG數據集上各準確率的匹配樣本數Tab．2 Matching number of accuracy of 3 pattern matching methods on ECG dataset

由表2可知，對于ECG數據集，當準確率較小時，本文方法的匹配個數遠小于文［9－10］的方法;當準確率較高時，本文方法的匹配個數遠大于文［9－10］方法．特別指出的是，當準確率為1時，本文方法在ECG數據集上三種近鄰對應的匹配個數比文［9－10］方法多60個左右．四組數據集的實驗結果表明本文方法對數據集的序列長度不限制，且匹配效果遠好于文［9－10］方法．進一步由式(16)～(17)計算得到四個數據集各自的準確率期望值，結果見表3．

表3 三種模式匹配算法在四組數據集上的準確率期望值Tab．3 Accuracy expectations of 3 pattern matching methods on 4 datasets

觀察表3可知，對四組數據集而言，在三種不同的近鄰數模式匹配情況下，本文方法的匹配期望準確率都比文［9－10］方法高．特別指出:對于ECG數據集，三種模式匹配狀態下，本文方法的準確率期望值比文［9－10］方法提高了約20%．對于LP1數據集，本文方法的準確率期望值約為77%，而文［9－10］的準確率期望值僅約為44%，這說明文［9－10］方法不適于序列長度較短的數據集，而本文方法則沒有這方面的局限．

3．4 三種模式匹配方法時間開銷比較

在四組數據集分別測試采用本文方法、文［9－10］方法的時間開銷，實驗結果如圖2所示．可以直觀地看出，在四組數據集上，本文方法的時間開銷都遠少于文［9－10］方法．

圖2 三種模式匹配算法在四組數據集上的時間開銷比較Fig．2 Computational time comparision of 3 pattern matching methods on 4 datasets

4 結論

提出基于張量多線性PCA的模式匹配方法．該方法使用張量多線性PCA對多變量時間序列進行低維重構，得到序列長度及變量個數均被減少的模式表示，并采用Frobenius范數度量矩陣形式表示的模式的相似度．在四組真實MTS數據集上實驗結果表明該方法比現有多變量時間序列模式匹配方法的匹配準確率更高，匹配時間更短且適用不同序列長度的數據集．

［1］Jeong Y S，Jeong M K，Omitaomu O A．Weighted dynamic time warping for time series classification［J］．Pattern Recognition，2011，44(9):2 231－2 240．

［2］Pham T D．Possibilistic nonlinear dynamical analysis for pattern recognition［J］．Pattern Recognition，2013，46(3):808 －816．

［3］丁永偉，楊小虎，陳根才，等．基于弧度距離的時間序列相似度量［J］．電子與信息學報，2011，33(1):122－128．

［4］Frank J，Mannor S，Pineau J，et al．Time series analysis using geometric template matching［J］．IEEE Transactions on Pattern Ananlysis and Machine Intelligence，2013，35(3):740－754．

［5］劉博寧，張建業，張鵬，等．基于曲率距離的時間序列相似性搜索方法［J］．電子與信息學報，2012，34(9):2200－2 207．

［6］Raquel P，Francisco M，Gabriel H．Multivariate time series modeling and classification via hierarchical VAR mixtures［J］．Computational Statistics＆ Data Analysis，2006，51(3):1 445－1 462．

［7］管河山，姜青山，王聲瑞．基于點分布特征的多元時間序列模式匹配方法［J］．軟件學報，2009，20(1):67－79．

［8］李正欣，張鳳鳴，李克武．基于DTW的多元時間序列模式匹配方法［J］．模式識別與人工智能，2011，24(3):425－430

［9］Singhal A，Seborg D E．Pattern matching in historical batch data using PCA［J］．IEEE Control Systems Magazine，2002，22(5):53－63．

［10］Yang Kiyong，Shahabi C．A PCA －based similarity measure for multivariate time series［C］//Proceeding of the Second ACM International Workshop on Multimedia Databases．Washington:ACM Press，2004:65 －74．

［11］王甦菁．流形上的張量子空間人臉識別算法的研究［D］．長春:吉林大學，2012．

［12］Lu Haiping，Plataniotis K N，Venetsanopoulos A N．MPCA:Multilinear principal component analysis of tensor objects［J］．IEEE Transactions on Neural Networks，2008，19(1):18 －39．

［13］Lee K，Park H．Probabilistic learning of similarity measures for tensor PCA［J］．Pattern Recognition Letters，2012，33(10):1 364－1 372．

［14］Mineichi K，Jun T，Masaru S．Japanese vowels［EB/OL］．2009．http://kdd．ics．uci．edu/databases/Japanese Vowels/Japanese Vowels．html．

［15］Lopes L S，Camarinha － Matos L M．Robot execution failures［EB/OL］．2009．http://kdd．ics．uci．edu/databases/robotfailure/robotfailure．html．