含噪聲的海量多元時間序列降維方法研究

劉 博，郭建勝

空軍工程大學 裝備管理與安全工程學院，西安 710051

1 引言

多元時間序列（Multivariate Time Series，MTS）是在同一個時間點上擁有多個觀測值的時間序列。為了得到更加豐富、準確的信息，MTS觀測量的維數不斷增加，在信息更加準確的同時，也帶來了“維災難”[1]的問題。因此，MTS的降維問題以成為目前時間序列分析領域的一個熱點問題。

主成分分析（Principal Component Analysis，PCA）[2]作為一種經典的線性降維算法，在數據挖掘領域得到了廣泛應用。PCA方法實際上就是尋找一組相互正交的方向向量（主成分），構成一個低維特征空間，使得原始數據在低維特征空間上的投影能夠最大程度地體現差異性[3]。但由于其產生的低維特征空間易受噪聲的影響，降低了其學習能力。為了降低噪聲對低維特征空間的影響，文獻[4]提出了Robust PCA方法，運用迭代法對經典的PCA加權，進而減小噪聲的影響。但Robust PCA法時間復雜度大，降低了算法對大樣本數據的學習能力。文獻[5]提出了一種基于角度優化的全局嵌入算法（Angle Optimized Global Embedding，AOGE），該方法雖然解決了噪聲造成的子空間偏移問題，但在非線性數據的降維方面還有待提高。

在MTS的降維研究中，文獻[6-9]均是先用PCA方法對數據集中的MTS分別進行降維，然后再進行后續分析。文獻[10]引入了共同特征空間的概念，提出了共同核主成份法（Common Kernel Principal Component Analysis，CKPCA）。但是這幾種方法均忽視了MTS的高噪聲性，影響了其降維效果。

針對上述問題，本文在結合CKPCA和AOGE方法優點的基礎上，提出了一種基于角度優化的共同核主成份法（Angle Optimized Common Kernel Principal Component Analysis，AOCKPCA），并將這種算法部署到Hadoop平臺上。該方法提高了含噪聲的非線性MTS降維的效果，同時也提高了算法對海量MTS降維的能力。最后，通過對比實驗，驗證了所提方法的有效性。

2 相關研究

2.1 MapReduce編程模式介紹

Hadoop平臺是目前主流的云計算平臺，為解決海量數據的降維問題提供了一種新的并行解決方案。Hadoop平臺的MapReduce編程模式，是Google Map-Reduce算法的開源實現，追求“Moving computation is cheaper than moving data”。在MapReduce編程模式中，開發者僅需設計輸入輸出的key/value格式，即可實現分布式運算，降低了開發人員的開發難度和強度。MapReduce具體執行過程見圖1。

圖1 MapReduce運行機制

由于MTS是由各時間點的觀察值組成的，且又相互獨立，所以MTS數據集可以方便的進行并行分析。依托Hadoop平臺，可以將海量的MTS存儲在不同的節點上，有利于提高分析效率。

2.2 CKPCA降維模型介紹

CKPCA降維方法是運用MTS數據集的共同特征空間，來消除各MTS特征空間差異給降維帶來的影響。共同特征空間是通過共同協方差矩陣得到的，共同協方差矩陣的計算公式：

其中，Ci是各MTS的協方差矩陣，l是MTS數據集的容量。該方法的具體步驟：首先通過核函數的非線性映射，將原始序列映射到高維特征空間；再運用點積與核函數的關系，求得共同特征空間；然后按照PCA方法進行降維。

2.3 AOGE降維模型介紹

AOGE方法是利用中心化樣本偏離其在低維空間的正交投影的角度來刻畫數據之間的關系，消除了基于距離度量對子空間帶來的誤差。AOGE降維算法構造了有約束的優化模型模型[5]：

3 AOCKPCA降維方法

AOCKPCA方法基本思路：先通過非線性映射，將原始序列映射到高維特征空間；再在高維特征空間上計算基于角度優化的共同子空間；然后再按照PCA方法進行降維。其中，在計算基于角度優化的共同子空間過程中得到的共同主成分被稱為基于角度優化的共同核主成分。

由于具體的映射形式和高維空間的維數未知，實現基于角度優化的共同核主成分降維的關鍵在于：如何將映射到高維特征空間中的數據向基于角度優化的共同核主成分方向的投影，轉化為高維空間的點積運算，從而通過核函數實現非線性降維。

3.1 AOCKPCA降維模型

在介紹AOCKPCA降維模型之前，先給出MTS的基本概念[6]：一系列觀測值xt(j)稱為多元時間序列，其中xt(j)表示第j(j=1，2，…，m)個觀測量在第t(t=1，2，…，n)個時間點上的觀測值。

下面給出AOCKPCA方法的數學推導。

再對已經中心化的觀測向量進行單位化，即

則第k個MTS在特征空間F中的M矩陣為：

在MTS數據集中，MTS通常具有相同的觀測量維數，且各觀測量一一對應，具有相同的含義，但時間長度一般不同[10]，所以Mk是一個m×m的矩陣。對l個MTS的M矩陣取平均為：

其中，αj為相關系數。不失一般性，令式（6）中k=1，將式（4）、（6）帶入式（5），同時左乘以轉置向量(x[1]r)可得：

為簡化表達式，下面引入兩個矩陣定義：

（1）N×N矩陣 Kh，稱為核矩陣。其第i行j列上的元素如式（8）所示，其中Kx，y表示φ(x)與φ(y)的內積。

（2）N×1向量a，第j個元素為參數aj。

通過引入核矩陣，式（7）可以簡化為：

一般情況下，式（1）不成立[10]，通常采用式（11）中的來代替 Kh。

其中IN是每一元素都為1N的N階方陣。

由上述模型的數學推導不難發現，若所有觀測向量的模都相同（不妨設所有觀測值的模均為常數），則核矩陣第i行第j列的元素將變為，此時模型（10）就是CKPCA模型?？梢姡珹OCKPCA算法是對CKPCA算法的推廣。

3.2 基于AOCKPCA的MTS降維

對高維共同特征空間F中的特征向量α進行規范化，計算公式為：

設x是一個待降維的觀測向量，在高維公共特征空間F中的映射后為φ(x)，則其向特征向量α投影為：

其中為ak第k各特征向量，是ak的第i個系數。

向高維公共特征空間F中的p個向量依次投影，從而形成降維后的P維向量：

下面給出AOCKPCA降維方法的具體步驟：

（1）選擇合適的核函數、確定核函數的形狀參數，確定方差貢獻率σ。

（2）根據式（8）計算數據集中每個MTS的核矩陣Kh={Kx(h)i，x(1)j}。

（4）求解式（10）表示的特征值問題，并找出全部的實數特征值λi和對應的特征向量αi，根據方差貢獻率σ得到p個基于角度優化的共同核主成分。

（5）由式（13），對步驟（4）中求得的p個基于角度優化的共同核主成分進行規范化處理。

（6）p個規范化后的基于角度優化的共同核主成分構成公共特征空間，根據式（14）和（15）將一個MTS依次往公共特征空間上投影，得到p(p＜m)維的MTS序列。

4 AOCKPCA降維方法的并行化改進

由第2.1節的介紹可知，基于Hadoop平臺的并行算法設計最主要的工作就是設計和實現Map和Reduce函數，包括輸入和輸出鍵值對的類型以及Map和Reduce函數的具體邏輯等。由3.2節的介紹可知，AOCKPCA降維方法的主要計算可以歸并成兩部分：一部分是計算共同核矩陣，另一部分是求得特征空間后的降維。因此本文分兩部分設計Map和Reduce函數，并將這兩對函數進行串連，達到并行降維的目的，具體流程如圖2所示，圖中每一個方框均包含獨立的MapReduce任務。

圖2 并行AOCKPCA的實現流程

下面給出并行AOCKPCA降維方法的主要設計。

4.1 并行計算共同核矩陣

（1）降維前期處理。在降維前需要首先選定核函數和確定方差貢獻率σ。同時為了適合MapReduce計算模型處理，須將待降維的MTS以時間點的向量形式存儲，即存儲為 ＜h，i，list:v＞ 的形式，其中h表示第h個MTS樣本，i表示第h個MTS樣本的第i個時間點，list:v為第h個MTS的第i個時間點的上的觀察值列表。

（2）Map函數設計。Map函數的任務是根據選定的核函數計算核矩陣各個元素的值。輸出為 ＜h， ＜i，j，v′＞＞，其中h是key， ＜i，j，v′＞ 是value，v′是第h個核矩陣第i行第j列元素的值。函數偽代碼描述如下：

如因單個MTS的觀測量維數過高，可以將觀測量分段，依靠多個并行的Map函數并行處理。例如一個MTS的觀測向量有900維，現想用3個Map函數并行處理，則可以用第1個Map函數處理第1～300號觀測量的相關運算；用第2個Map函數處理第301～600號觀測量的相關運算；用第3個Map函數處理第601～900號觀測量的相關運算。

同樣，如因MTS的時間維數太高，也可按同樣的方法進行再并行處理。當然，也可以將MTS同時從兩個方向劃分，每個Map僅計算其中的一塊。

由于MapReduce編程模式的特點，Map函數僅需考慮劃分的問題，合并的問題應交于Reduce函數處理。

（3）Reduce函數設計。Hadoop平臺可以自動將不同key的Map函數輸出傳遞給不同的Reduce函數[13]。這里因Map輸出的key是MTS的編號，所以Hadoop平臺自動將的不同的MTS的核矩陣的元素值傳遞給不同的Reduce函數。然后Reduce函數的主要工作是合并Map函數的運算結果，根據式（11）計算----Kh各元素的值，函數偽代碼描述如下：

以上這個Reduce函數是在未建立并行的Map函數的前提下運行的。如前面建立了并行的Map函數，可在開始計算前先將并行的Map運算結果進行合并，然后再運算。例如前面建立了并行的3個Map函數，則在Reduce開始運算前，可先將3個Map函數的對應的輸出結果合并，然后再開始計算。

4.2 降維

（1）計算公共特征空間。這一步是在求得共同核矩陣之后進行的，僅需按常規的方法調用API函數編寫一個函數即可，這里僅給出具體執行步驟：

①求解式（10）表示的特征值問題，并找出全部的實數特征值λi和對應的特征向量αi。

②根據式（12），得到p個基于角度優化的共同核主成分。

③由式（13）對基于角度優化的共同核主成分進行規范化，構成共同特征空間F。

（2）Map函數設計。在得到公共特征空間后在串聯一個Map函數，其任務是根據式（14）和式（15）得到降維序列。輸入為 ＜h，i，list:v＞ ，輸出為 ＜null， ＜h，i，list:val′＞＞ ，其中list:val′是降維后第h個核矩陣第i行的各值列表，因是最終輸出結果，未設key值。函數偽代碼描述如下：

因Map函數降維后的輸出結果為獨立的向量，有利于進一步的并行分析，所以本文未將輸出結果用Reduce函數合并成MTS。

5 實驗結果與分析

5.1 實驗環境

本實驗選取5臺計算機搭建Hadoop平臺，每臺機器都是雙核IntelCorei3處理器，1 GB內存，操作系統采用國產麒麟操作系統，Hadoop版本選用Hadoop 0.20.2；jdk的版本是1.6.30；一臺機器作為Master和JobTracker服務節點，其他機器作為Slave和TaskTracker服務節點。

5.2 實驗數據集選取

Australian Sign Language（記為 ASL）數據集包含95種語意（95個類），每種語意都有27組序列，每個序列包含22個變量。實驗中σ=80%，選用RBF核函數。

5.3 實驗

（1）降維性能實驗：降維的目的是為了更有效地利用數據，所以本文采用模式識別的方法測試算法降維的性能。實驗采用動態時間彎曲（Dynamic Time Warping，DTW）模式匹配模型，進行k-近鄰模式識別。實驗在10種語意中，隨機選取20組序列作為訓練樣本，剩下7組序列作為測試樣本。測試指標為分類準確率e=n0/n，其中n0為識別準確數，n為待識別樣本數。為了測試算法對含噪聲數據的分類效果，首先對數據加以均值為0，方差為0.1的高斯白噪聲，實驗共做了10次，實驗結果見圖3。

從實驗結果可以看出，在對于含有噪聲的MTS進行模式識別中，AOCKPCA方法較CKPCA方法有效，有較強的抗噪聲能力。

圖3 含噪聲的MTS模式識別率對比

（2）集群加速比實驗：加速比是同一個任務在單處理器系統和并行處理器系統中運行消耗時間的比率，是衡量并行系統或并行化程序性能和效果的主要指標。本文對比了AOCKPCA算法和基于Hadoop平臺的AOCKPCA算法處理不同大小數據集時的加速比，實驗選取了10種語意，30種語意，60種語意，90種語意的4組數據進行測試，結果如表1所示。

表1 AOCKPCA方法的加速比

從實驗結果可以看出，當集群節點數增加時，系統在處理相同數量的作業時性能逐漸提升，當集群節點不變時，系統的性能也穩步提升。但實驗中發現，當實驗數據較少（10中語意）時，系統性能不如單機性能，這是因為MapReduce編程模型本身存在在大量的I/O操作，影響了運算速度。

（3）集群伸縮性實驗：為了驗證算法的可伸縮性，分別在由5臺和3臺機器組成的集群上進行測試并與傳統AOCKPCA算法進行對比，實驗選取了10種語意，30種語意，60種語意，90種語意的4組數據進行測試，結果如表2所示。

表2 運算時間對比 ms

從表中的數據可以看出，隨著數據量的增加，AOCKPCA方法的效率優勢會逐步體現，數據集規模越大，優勢越明顯。當集群中節點數的增多時，算法效率隨之增加，節點越多，對數據集規模的增大越不敏感。

5.4 結果分析

通過以上實驗可見，AOCKPCA方法較CKPCA方法在針對含噪聲的MTS降維運算時，具有更好的表現；在部署到Hadoop平臺后，本文所提的算法擁有較好的運行效率和伸縮性，可以有效地運用于含噪聲的海量MTS的降維運算。

6 結束語

針對含噪聲的海量MTS的特點，提出了AOCKPCA降維算法，并將該算法部署到Hadoop平臺，通過實驗進行了檢驗。實驗結果表明，在對于含噪聲的海量MTS進行降維時，AOCKPCA算法較CKPCA算法更有效，部署到Hadoop平臺上后，并行的AOCKPCA方法的執行效率得到了大幅提高，具有良好的加速比和可移植性。但由于MapReduce運算時存在大量的I/O操作，對算法的運行效率影響較大，需進一步優化Map和Reduce兩個函數的輸入輸出，進一步提高運行效率，并需將AOCKPCA算法運行在更大規模的Hadoop集群中以檢驗其處理海量數據的性能。

[1]Keogh E，Pazzani M.A simple dimensionality reduction technique for fast similarity search in large time series databases[C]//Proceedings of the 4th Pacific-Asia Conference on Knowledge Discovery and Data Mining，Kyoto，2000.

[2]Jolliffe I T.Principal component analysis[M].New York：Springer，2002.

[3]Boente G，Pires A M，Rodrigues I M.Detecting influential observations in principal components and common principalcomponents[J].ComputationalStatistics and Data Analysis，2010，54（12）：2967-2975.

[4]Dela Torre F，Black M J.Robust principal component analysis for computer vision[C]//Proceedings of the 8th IEEE InternationalConference on ComputerVision，Vancouver，Canada.[S.l.]：IEEE，2001：362-369.

[5]劉勝藍，閆德勤.一種新的全局嵌入降維算法[J].自動化學報，2011，37（7）：828-835.

[6]張軍，吳紹春，王煒.多變量時間序列模式挖掘的研究[J].計算機工程與設計，2006，27（18）：3364-3366.

[7]楊興江，周勇.多元時間序列相似性研究[J].西南民族大學學報，2007，33（4）：864-869.

[8]周大鐲，吳曉麗，閆紅燦.一種高效的多變量時間序列相似查詢算法[J].計算機應用，2008，28（10）：2541-2543.

[9]周大鐲，姜文波，李敏強.一個高效的多變量時間序列聚類算法[J].計算機工程與應用，2010，46（1）：137-139.

[10]李正欣，張鳳鳴，張曉豐，等.多元時間序列特征降維方法研究[J].小型微型計算機系統，2013，34（2）：338-344.

[11]BhandarkarM.MapReduceprogrammingwithapache Hadoop[C]//Proceedings of IEEE International Symposium on Parallel＆Distributed Processing，2010.

[12]Cheema S，Henne T，Koeckemann U，et al.Applicability of feature selection on multivariate time series data for robotic discovery[C]//Proc of 3rd International Conference on Advanced Computer Theory and Engineering，2010：592-597.

[13]Ziyatdinov A，Marco S，Chaudry A，et al.Drift compensation of gassensor array data by common principal component analysis[J].Sensors and Actuators B Chemical，2010，146（2）：460-465.

[14]Lam C.Hadoop實戰[M].北京：人民郵電出版社，2011：85-112.

[15]White T.Hadoop：the definitive guide[M].[S.l.]：O’Reilly Media，2009.