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

劉 博，郭建勝

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

1 引言

多元時(shí)間序列（Multivariate Time Series，MTS）是在同一個(gè)時(shí)間點(diǎn)上擁有多個(gè)觀測值的時(shí)間序列。為了得到更加豐富、準(zhǔn)確的信息，MTS觀測量的維數(shù)不斷增加，在信息更加準(zhǔn)確的同時(shí)，也帶來了“維災(zāi)難”[1]的問題。因此，MTS的降維問題以成為目前時(shí)間序列分析領(lǐng)域的一個(gè)熱點(diǎn)問題。

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

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

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

2 相關(guān)研究

2.1 MapReduce編程模式介紹

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

圖1 MapReduce運(yùn)行機(jī)制

由于MTS是由各時(shí)間點(diǎn)的觀察值組成的，且又相互獨(dú)立，所以MTS數(shù)據(jù)集可以方便的進(jìn)行并行分析。依托Hadoop平臺(tái)，可以將海量的MTS存儲(chǔ)在不同的節(jié)點(diǎn)上，有利于提高分析效率。

2.2 CKPCA降維模型介紹

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

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

2.3 AOGE降維模型介紹

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

3 AOCKPCA降維方法

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

由于具體的映射形式和高維空間的維數(shù)未知，實(shí)現(xiàn)基于角度優(yōu)化的共同核主成分降維的關(guān)鍵在于：如何將映射到高維特征空間中的數(shù)據(jù)向基于角度優(yōu)化的共同核主成分方向的投影，轉(zhuǎn)化為高維空間的點(diǎn)積運(yùn)算，從而通過核函數(shù)實(shí)現(xiàn)非線性降維。

3.1 AOCKPCA降維模型

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

下面給出AOCKPCA方法的數(shù)學(xué)推導(dǎo)。

再對已經(jīng)中心化的觀測向量進(jìn)行單位化，即

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

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

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

為簡化表達(dá)式，下面引入兩個(gè)矩陣定義：

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

（2）N×1向量a，第j個(gè)元素為參數(shù)aj。

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

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

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

由上述模型的數(shù)學(xué)推導(dǎo)不難發(fā)現(xiàn)，若所有觀測向量的模都相同（不妨設(shè)所有觀測值的模均為常數(shù)），則核矩陣第i行第j列的元素將變?yōu)?，此時(shí)模型（10）就是CKPCA模型?？梢?，AOCKPCA算法是對CKPCA算法的推廣。

3.2 基于AOCKPCA的MTS降維

對高維共同特征空間F中的特征向量α進(jìn)行規(guī)范化，計(jì)算公式為：

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

其中為ak第k各特征向量，是ak的第i個(gè)系數(shù)。

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

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

（1）選擇合適的核函數(shù)、確定核函數(shù)的形狀參數(shù)，確定方差貢獻(xiàn)率σ。

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

（4）求解式（10）表示的特征值問題，并找出全部的實(shí)數(shù)特征值λi和對應(yīng)的特征向量αi，根據(jù)方差貢獻(xiàn)率σ得到p個(gè)基于角度優(yōu)化的共同核主成分。

（5）由式（13），對步驟（4）中求得的p個(gè)基于角度優(yōu)化的共同核主成分進(jìn)行規(guī)范化處理。

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

4 AOCKPCA降維方法的并行化改進(jìn)

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

圖2 并行AOCKPCA的實(shí)現(xiàn)流程

下面給出并行AOCKPCA降維方法的主要設(shè)計(jì)。

4.1 并行計(jì)算共同核矩陣

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

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

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

同樣，如因MTS的時(shí)間維數(shù)太高，也可按同樣的方法進(jìn)行再并行處理。當(dāng)然，也可以將MTS同時(shí)從兩個(gè)方向劃分，每個(gè)Map僅計(jì)算其中的一塊。

由于MapReduce編程模式的特點(diǎn)，Map函數(shù)僅需考慮劃分的問題，合并的問題應(yīng)交于Reduce函數(shù)處理。

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

以上這個(gè)Reduce函數(shù)是在未建立并行的Map函數(shù)的前提下運(yùn)行的。如前面建立了并行的Map函數(shù)，可在開始計(jì)算前先將并行的Map運(yùn)算結(jié)果進(jìn)行合并，然后再運(yùn)算。例如前面建立了并行的3個(gè)Map函數(shù)，則在Reduce開始運(yùn)算前，可先將3個(gè)Map函數(shù)的對應(yīng)的輸出結(jié)果合并，然后再開始計(jì)算。

4.2 降維

（1）計(jì)算公共特征空間。這一步是在求得共同核矩陣之后進(jìn)行的，僅需按常規(guī)的方法調(diào)用API函數(shù)編寫一個(gè)函數(shù)即可，這里僅給出具體執(zhí)行步驟：

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

②根據(jù)式（12），得到p個(gè)基于角度優(yōu)化的共同核主成分。

③由式（13）對基于角度優(yōu)化的共同核主成分進(jìn)行規(guī)范化，構(gòu)成共同特征空間F。

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

因Map函數(shù)降維后的輸出結(jié)果為獨(dú)立的向量，有利于進(jìn)一步的并行分析，所以本文未將輸出結(jié)果用Reduce函數(shù)合并成MTS。

5 實(shí)驗(yàn)結(jié)果與分析

5.1 實(shí)驗(yàn)環(huán)境

本實(shí)驗(yàn)選取5臺(tái)計(jì)算機(jī)搭建Hadoop平臺(tái)，每臺(tái)機(jī)器都是雙核IntelCorei3處理器，1 GB內(nèi)存，操作系統(tǒng)采用國產(chǎn)麒麟操作系統(tǒng)，Hadoop版本選用Hadoop 0.20.2；jdk的版本是1.6.30；一臺(tái)機(jī)器作為Master和JobTracker服務(wù)節(jié)點(diǎn)，其他機(jī)器作為Slave和TaskTracker服務(wù)節(jié)點(diǎn)。

5.2 實(shí)驗(yàn)數(shù)據(jù)集選取

Australian Sign Language（記為 ASL）數(shù)據(jù)集包含95種語意（95個(gè)類），每種語意都有27組序列，每個(gè)序列包含22個(gè)變量。實(shí)驗(yàn)中σ=80%，選用RBF核函數(shù)。

5.3 實(shí)驗(yàn)

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

從實(shí)驗(yàn)結(jié)果可以看出，在對于含有噪聲的MTS進(jìn)行模式識別中，AOCKPCA方法較CKPCA方法有效，有較強(qiáng)的抗噪聲能力。

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

（2）集群加速比實(shí)驗(yàn)：加速比是同一個(gè)任務(wù)在單處理器系統(tǒng)和并行處理器系統(tǒng)中運(yùn)行消耗時(shí)間的比率，是衡量并行系統(tǒng)或并行化程序性能和效果的主要指標(biāo)。本文對比了AOCKPCA算法和基于Hadoop平臺(tái)的AOCKPCA算法處理不同大小數(shù)據(jù)集時(shí)的加速比，實(shí)驗(yàn)選取了10種語意，30種語意，60種語意，90種語意的4組數(shù)據(jù)進(jìn)行測試，結(jié)果如表1所示。

表1 AOCKPCA方法的加速比

從實(shí)驗(yàn)結(jié)果可以看出，當(dāng)集群節(jié)點(diǎn)數(shù)增加時(shí)，系統(tǒng)在處理相同數(shù)量的作業(yè)時(shí)性能逐漸提升，當(dāng)集群節(jié)點(diǎn)不變時(shí)，系統(tǒng)的性能也穩(wěn)步提升。但實(shí)驗(yàn)中發(fā)現(xiàn)，當(dāng)實(shí)驗(yàn)數(shù)據(jù)較少（10中語意）時(shí)，系統(tǒng)性能不如單機(jī)性能，這是因?yàn)镸apReduce編程模型本身存在在大量的I/O操作，影響了運(yùn)算速度。

（3）集群伸縮性實(shí)驗(yàn)：為了驗(yàn)證算法的可伸縮性，分別在由5臺(tái)和3臺(tái)機(jī)器組成的集群上進(jìn)行測試并與傳統(tǒng)AOCKPCA算法進(jìn)行對比，實(shí)驗(yàn)選取了10種語意，30種語意，60種語意，90種語意的4組數(shù)據(jù)進(jìn)行測試，結(jié)果如表2所示。

表2 運(yùn)算時(shí)間對比 ms

從表中的數(shù)據(jù)可以看出，隨著數(shù)據(jù)量的增加，AOCKPCA方法的效率優(yōu)勢會(huì)逐步體現(xiàn)，數(shù)據(jù)集規(guī)模越大，優(yōu)勢越明顯。當(dāng)集群中節(jié)點(diǎn)數(shù)的增多時(shí)，算法效率隨之增加，節(jié)點(diǎn)越多，對數(shù)據(jù)集規(guī)模的增大越不敏感。

5.4 結(jié)果分析

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

6 結(jié)束語

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

[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]劉勝藍(lán)，閆德勤.一種新的全局嵌入降維算法[J].自動(dòng)化學(xué)報(bào)，2011，37（7）：828-835.

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

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

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

[9]周大鐲，姜文波，李敏強(qiáng).一個(gè)高效的多變量時(shí)間序列聚類算法[J].計(jì)算機(jī)工程與應(yīng)用，2010，46（1）：137-139.

[10]李正欣，張鳳鳴，張曉豐，等.多元時(shí)間序列特征降維方法研究[J].小型微型計(jì)算機(jī)系統(tǒng)，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實(shí)戰(zhàn)[M].北京：人民郵電出版社，2011：85-112.

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