基于雙重正則矩陣分解的缺失數據恢復

劉 歌, 芮國勝, 田文飚

(海軍航空大學, 山東 煙臺 264001)

0 引 言

現實場景中,在某一監測區域內通過部署多個傳感器對同一對象進行持續的感知,獲取豐富的信息以支撐不同類型的感知應用[1]。這些多傳感器網絡中收集到的數據通常被稱為多源時間序列。例如,海上浮標多個傳感器監測海洋環境(maritime environment, ME)數據(溫度、濕度、壓強、風速、風向等)以獲取蒸發波導整體態勢感知[2];個人醫療系統中通過穿戴設備布設多個傳感器監測血壓、脈搏、心電等數據獲知病人的整體健康狀況[3]。由于惡劣的工作條件或不可控制的因素,導致傳感器網絡短暫甚至長時間無法收集數據從而出現信息盲區,影響系統的感知應用。將上述問題歸結為多源時間序列數據的缺失信息恢復問題。該問題廣泛存在于海上氣象要素監測傳感器網絡[4]、電網系統[5]、物聯網系統[6]、森林防火監測系統[7]等各種傳感器網絡中,因此解決缺失數據的恢復問題具有重要的實際應用價值。

缺失數據恢復方法中最簡單的是插值法,文獻[8]提出了基于M分量的概率主成分分析模型和期望最大化(M-component probabilitic principal component analysis-expectation maximization,MPPCA-EM)算法的多重插補算法來預測歷史降雨時間序列數據,該方法簡單易操作,但是僅適用于丟失少量數據并且時間序列變化非常穩定的情況[9]。常用的方法還有建模法,這是一種通過掌握數據內部的隱含規律,建立模型以預測缺失數據的方法。Frasconi[10]等使用季節性核來測量時間序列實例之間的相似性,并提出將季節性自回歸綜合移動平均模型與卡爾曼濾波器結合使用,可實現丟失數據預測性能。文獻[11]使用了一種基于廣義空間動態自回歸模型(generalized spatial dynamic autoregression model, GSDAM)的方法對缺失數據進行預測。但是基于模型的方法局限性比較大,一旦脫離對應的數據類型,模型將失效。近年來,基于矩陣填充和矩陣分解(matrix factorization, MF)的方法逐漸興起。矩陣填充方法中各種低秩約束的算法如奇異值閾值(singular value thresholding, SVT)算法[12]、奇異值投影(singular value projection, SVP)算法[13]、基于近似SVD的FPC算法(approximate SVD based FPC algorithm, FPCA)[14]等能夠實現缺失數據的恢復,但是需要滿足矩陣低秩性和非相干性以及R-RIP條件,約束條件嚴格。Song等用矩陣分解法預測流量矩陣,其方法比傳統方法表現出更有效的性能[15],但是未能充分利用數據內部的先驗信息,預測性能有待提升。文獻[16]時域動態矩陣分解(time-domain dynamic matrix factorization, TDMF)的方法,通過施加不同的正則化項來約束矩陣分解的目標函數,并建立了5個相應的模型來預測多變量時間序列中的缺失數據,但是對于傳感器數據的相關性計算忽略了數據的變化趨勢相關性,可能導致相似傳感器的判斷不夠準確。

為解決上述問題,本文提出一種基于雙重正則矩陣分解(double regularization matrix factorization, DRMF)的方法,該方法明確地考慮了多源時間序列在時間和傳感器兩個角度下的數據特性,充分挖掘先驗信息。對于時間序列,通過定義時間序列的穩定性,利用二階差分正則化進行約束;對于多傳感器數據,引入圖論的基本原理,設計了一種基于雙重皮爾遜系數的相關性度量策略,獲取表示各傳感器數據關系的拉普拉斯矩陣;最終將圖拉普拉斯正則化與二階差分正則化融入到的矩陣分解框架下,利用梯度下降法實現目標函數的優化。

1 問題描述

傳感器網絡在其網絡分布區域收集到的各種ME監測數據通常以矩陣形式表示,在本文中將該矩陣稱為多源時間序列矩陣,記作X(X∈Rm×n),其元素xij表示第i(i=1,2,…,m)個氣象傳感器源在第j(j=1,2,…,n)個時刻獲得的某一氣象要素數據。缺失信息矩陣X具體可表示為

X=

式中,符號“”表示該位置出現缺失數據;xi∈Rn表示第i個氣象傳感器源獲取的數據向量;tj表示第j個采樣時刻。

2 矩陣分解模型

矩陣分解是矩陣恢復中常用的方法,其優點在于發現數據中的潛在結構,并且可以用概率解釋;容易擴展到一些指定特定先驗信息的領域;可以用許多優化方法例如梯度下降法來找到一個最優解。奇異值分解是其中最為普遍的一種分解方法,給定X∈Rm×n,X的奇異值分解為

X=UΣVT

(1)

式中,U∈Rm×r表示傳感器源隱含因子矩陣;V∈Rn×r表示時間序列隱含因子矩陣;Σ∈Rr×r為對角陣,對角元素是矩陣X奇異值;r表示隱含因子的維度。

式(1)可以轉換為兩個矩陣相乘的表示形式,即

X=SQT

(2)

為了后續計算方便,定義一個矩陣I∈Rm×n,用來表示矩陣X的缺失位置,即其中元素值為0表示此處信息缺失,元素值為1則相反。將缺失矩陣看作是完整矩陣與缺失位置矩陣的哈達瑪乘積,即

X′=I°X

(3)

式中,X′表示缺失信息的矩陣。為實現矩陣X的信息恢復,需要完成最優化問題:

(4)

為了防止過擬合,分別對兩個隱含因子進行約束:

(5)

式中,λS和λQ表示正則化參數。

3 算法實現

3.1 時間序列的正則化設計

現實場景中,各傳感器數據通常呈現為有限時間內的緩慢變化,把這種一定時間段內的緩變現象稱為時間序列的穩定性。

多源時間序列矩陣X的行向量為某一傳感器的時間序列,計算xij前后兩個相鄰位置的差值并作歸一化處理,得

(6)

以海洋氣象數據的獲取為例,海上傳感器網絡通過在指定海域部署大量海上傳感器節點,采集網絡分布區域內的各種環境要素,然后通過各類通信手段將采集到的數據傳送到監測中心。TAO/TRITON和PIRATA的浮標測量海洋和海表面氣象參數多達20余種,本文限于篇幅,僅對幾種典型數據進行分析,即蒸發波導特性研究所需要的氣溫、風速、海表溫度、壓強、相對濕度等環境要素。氣溫、風速、海表溫度、壓強和相對濕度5種海洋氣象環境要素數據的統計累積分布如圖1所示。

圖1 5種ME數據時間序列穩定性分析

從圖1中曲線可以看出,5種傳感器的時間序列在一定的時間范圍內變化緩慢,其中溫度和海表溫度以及濕度3種傳感器數據的r(i)≤0.1的比重占90%以上,風速和壓強兩種傳感器r(i)≤0.2的比重占90%以上。所以,多源傳感器時間序列的穩定性可以作為先驗信息加入矩陣分解的框架中。

基于上述分析,將時間序列隱含因子二階差分的正則化約束引入到目標函數:

(7)

式中,H表示二階差分矩陣;α表示正則化參數，有

(8)

3.2 多源傳感器隱含因子的正則化設計

在多源傳感器網絡中,盡管各傳感器的任務不同,獲取的數據單位不同,但是其目標相同,因此某些傳感器之間可能存在很強的相關性。例如,環境監測傳感器網絡中,溫度和濕度之間可能存在著較強的相關性;個人醫療監護系統中,心跳和血壓之間也可能存在高度相關性。這些相關性的信息都是多傳感器元數據的先驗信息。

為了更好地將多源傳感器時間序列的先驗信息融入到矩陣分解的框架中,受到圖論被廣泛用于探索數據內部的幾何結構的啟發,本文引入圖正則化對多傳感器數據之間的相互關系進行約束。從文獻[17]和文獻[18]中獲知,如果兩組數據在數據分布的內在幾何結構中很接近,那么這兩組數據相對于基函數的表示也很接近。這種假設通常被稱為局部不變假設。因此對多源時間序列,式(2)還可以表示為

(9)

式中,sj=[sj1,sj2,…,sjr]表示傳感器隱含因子矩陣S的行向量;qi表示矩陣Q的列向量;xj表示矩陣X的行向量。式(9)可以近似看作由S分量加權的Q列向量的線性組合。換個角度來看,將qi看作是基向量,[sj1,sj2,…,sjr]就是此基向量下的表示。因此,若兩個傳感器時間序列數據本身是相似的,那么傳感器隱含因子也是相似的。針對給定的大小為m的傳感器數據集{X1,X2,…,Xm},為這個數據集構造一個具有m個結點的無向圖,節點之間邊的權重值則表示兩節點之間的相似程度,這樣就為數據集構建了一個最近鄰圖。首先構造近鄰圖,如果xi,xj是近鄰關系,就在i節點和j節點之間添加一條邊,然后確定近鄰圖邊上的權重值,得到權重矩陣W。基于此,兩個傳感器數據之間的相似關系可以表示為

(10)

式中,si和sj表示第i個傳感器隱含因子和第j個傳感器隱含因子;Wij表示權重,具體表示為

(11)

式中,si和sj滿足的相似條件將在下面進行詳細探討。

2tr(ST(D-W)S)=2tr(STLS)

(12)

將傳感器隱含因子的拉普拉斯正則化加入到式(4)中,可以得到

(13)

式中,γ表示正則化參數。

基于此,為更好地說明不同傳感器數據之間的相互關系,本文從兩個傳感器數據之間的相關性和時間序列變化趨勢之間的關系入手,利用統計學中的皮爾遜相關系數(Pearson correlation coefficient, PCC)來度量兩個變量之間的相關程度,其表達式為

(14)

而僅采用 PCC計算兩個傳感器數據樣本之間相似性,忽略了變化趨勢的相關性,因此增加對數據變化趨勢的相似程度對比對于相似傳感器的判定,會更有說服力。基于上述分析,本文提出聯合一階擬合系數的PCC和數據本身的PCC的相似傳感器判定方法,為表示方便,將兩次求PCC的方法在本文中稱為D-PCC。

一階擬合系數是指通過對各傳感器時間序列進行分段擬合得到的系數組成的一組向量,這組向量就是傳感器數據變化趨勢的真實反應。因此,利用任意兩個傳感器的一階擬合系數來感知數據變化趨勢的相似性。假設兩個傳感器時間序列將Xi=[xi1,xi2,…,xin]和Xj=[xj1,xj2,…,xjn]均分為z段,分別對每一段進行一階擬合,得到的系數組分別為a=[a1,a2,…,az]和b=[b1,b2,…,bz],因此一階擬合系數的相似度可以計算為

(15)

因此,定義當且僅當任意兩個傳感器感知數據的綜合相關系數Ci, j

相似傳感器的定義建立在對兩個傳感器數據本身大小和變化趨勢兩者比較的基礎之上,該定義為衡量任意兩個傳感器之間是否具有穩定的數據相關性提供了依據。兩個傳感器的數據達到該衡量標準,無向圖中就可以在相似傳感器之間加一條邊,鄰接矩陣W就是對圖中的邊進行統計。

3.3 多源時間序列缺失信息恢復模型確立

第3.1節和第3.2節的分析,旨在從傳感器和時間序列兩個角度充分挖掘多源時間序列的先驗,對于提高缺失信息的恢復性能具有重要作用,因此本小節將兩種先驗約束統一于矩陣分解的框架下,提出基于時序差分正則化和圖拉普拉斯正則化的矩陣分解方法用于解決缺失數據的恢復問題。目標函數的最優化問題可以重新歸結為

(16)

通過梯度下降來優化式(16)所述的目標函數,即

(17)

(18)

結合上述分析和計算,將算法步驟總結如算法1所示來解決式(16)中所示的問題。在給定多變量時間序列X以及各種參數的條件下,設計該算法以獲得隱含因子特征矩陣S和Q更精確的解。算法更新的停止條件可以選擇預先設置的迭代次數,也可以是前后兩次迭代的誤差小于某一閾值,每次迭代都根據回溯線搜索策略更新步長τS和τQ,最終缺失數據可以從恢復矩陣中得到。

算法 1 DRMF算法輸入 多源時間序列X,指示矩陣I,參數r,λS,λQ,α,β,步長τ,閾值c,迭代次數iters輸出 恢復后的多源時間序列^X步驟1 獲取拉普拉斯矩陣(1) 綜合相關系數判別法確定權重矩陣W。根據式(14)、式(15)計算傳感器i和j的PCC_ab和PCC_X,從而得到Ci,j=12·(PCC_ab+PCC_X)當且僅當Ci,j

4 算法性能分析及仿真實驗

4.1 算法可行性分析

仿真實驗中，在完整數據的基礎上剔除一部分數據，剔除數據的數目占原始數據總數目的比例稱為信息缺失率。數據缺失的類型根據產生原因不同大體上可以分為隨機型、均勻型和連續型3類。在進行D-PCC相關系數計算時，對隨機型和均勻型缺失來說，缺失數據或隨機或均勻地分布在數據中，將缺失數據剔除，剔除前后的曲線變化趨勢沒有受到很大影響。將剔除缺失值后的數據曲線與原始數據曲線進行對比，如圖2所示。數據缺失類型為隨機缺失，缺失率為0.7。

圖2 原始時間序列與剔除缺失數據序列的曲線對比

從兩幅圖的對比中可以看出，剔除缺失數據前后的曲線變化趨勢極其相似，說明剔除缺失數據不但不會對整體的相關計算產生影響，反而因為缺失位置上的數據點的剔除，實現了數據的整體縮減，使得后續計算的數據量減少。從另一個角度來看，這樣做也提高了處理較長時間跨度內時間序列的能力。對于均勻缺失也是同樣的情況。不同的是連續型缺失，這種類型的缺失情況可能導致一段時間內的數據缺失，而其他時間數據沒有發生缺失，這種情況下將缺失數據的位置剔除，獲得的曲線與前兩種類型不同，因此不能直接用來進行相關性計算。針對這種情況，本文采用分階段相關計算的方法，以缺失位置為界，前后分別對相應位置上的數據進行相關性計算和曲線變化趨勢相關性計算，其余計算方法不變。

4.2 算法收斂條件

若要保證梯度下降法始終有效收斂，需要在正確范圍內選擇迭代過程中的步長。下面通過Lipschitz條件對梯度下降過程中的步長上界進行推導。

定理 1(β平滑) 若函數f(x)滿足Lipschitz條件，則存在一個常量β(β>0)，使f(x)在定義域上的任意2個值滿足:

(19)

定理 2(收斂條件) 假設凸函數f(x)滿足β平滑條件，以xn為當前點，以x*為最優點，如果τn≤1/β，則‖xn-x*‖隨n的減小而減小。

現在目標函數如式(16)所示，那么首先根據定理1推導變量S和Q的β平滑條件。對于變量S，將式(17)代入式(19)中得

‖(S1QT-X)Q+λSS1+γ(L+LT)S1-

(S2QT-X)Q+λSS2+γ(L+LT)S2‖=

‖(S1-S2)QTQ+(λSI+γ(L+LT))(S1-S2)‖≤

‖(S1-S2)QTQ‖+‖(λSI+γ(L+LT))(S1-S2)‖≤

‖S1-S2‖‖QTQ‖+‖λSI+γ(L+LT)‖‖S1-S2‖=

(‖QTQ‖+‖λSI+γ(L+LT)‖)‖S1-S2‖

(20)

由于此時‖QTQ‖+‖λSI+γ(L+LT)‖是已知的數據，所以令β=‖QTQ‖+‖λSI+γ(L+LT)‖且β>0時，E(S)滿足L-Lipschitz條件，若想保證式(17)的整個迭代過程收斂，根據定理2，第t次迭代過程中的步長滿足:

同理，對于變量Q，將式(18)代入式(19)中得

‖(Q1-Q2)STS+λQ(Q1-Q2)+α(Q1-Q2)HTH‖=

‖(Q1-Q2)(STS+λQI-αHTH)‖≤

‖STS+λQI-αHTH‖‖Q1-Q2‖

(21)

同樣，由于此時‖STS+λQI-αHTH‖是已知的數據，所以β=‖STS+λQI-αHTH‖且β>0時，E(Q)滿足L-Lipschitz條件。若要保證式(18)的整個迭代過程收斂，根據定理2，第t次迭代過程中的步長應滿足:

4.3 不同算法性能對比

為了評估該方法的性能，利用均方根誤差(root mean square error，RMSE)來衡量信息恢復質量。RMSE計算方式為

(22)

本文以3種數據集合作為實驗對象，以RMSE為評價標準，對算法性能進行測試，每次試驗得到的結果都是50次結果的平均值。

實驗部分采用一個合成的MTS數據集和兩個真實的MTS數據集。

(1) 合成(synthetic, SYN)數據集：由Asin(ωt+φ)+n產生，其中振幅A取值范圍為[1,2,3]，角頻率ω取值范圍為[1,π,2π]，初相φ的取值范圍為[0,π/2,π]，n為加性高斯白噪聲，時間序列長度為104。

(2) 海洋環境數據集：來自國內山東煙臺芝罘島氣象站的ME數據集，含氣溫、風速、海表溫度、壓強和相對濕度的每小時測量值，時間范圍為2017年1月1日20時-2018年9月4日23時，時間分辨率為1 h，取時間序列長度為104。

(3) Motes數據集：來自英特爾-伯克利實驗室54個MieahDot傳感器獲得的數據，時間從2004年2月28日-2004年4月5日。選取其中23個傳感器數據進行分析，時間序列長度為104。

利用算法對剔除的數據進行預測，然后通過預測值與原實際值的對比來評價預測算法的性能。依據相關文獻中通常采用的數據劃分方式，將數據集劃分為90%的訓練集和10%的測試集。

為方便計算，令正則化參數λS=λQ=0.1。對于SYN數據集、ME數據集以及Motes數據集，3組數據梯度下降法的步長在每次迭代中都根據回溯線搜索策略更新。下面對其他參數設置進行具體討論。實驗中具體的參數設置如表1所示。

表1 DRMF算法在不同數據集下的參數

4.3.1 算法恢復性能對比

將DRMF算法與其他4種算法進行仿真比較，以對算法性能進行有效評估，對比算法包括：基于GSDAM的算法、基于SVP的方法、基于TDMF的算法以及不添加正則項的MF算法。

DRMF算法與對比算法在缺失率不同的條件下的缺失信息恢復性能對比如圖3所示，通過對圖中數據的分析可以得到以下結論。

圖3 不同數據集下算法性能對比

從曲線的總體趨勢來看，DRMF算法顯示出更低的RMSE，證明了基于DRMF的方法是恢復多變量時間序列缺失值的有效方法。隨著缺失率的升高，誤差逐漸增大，但是在缺失率為90%時，誤差仍處于合理范圍內。具體來說，不同數據集中，DRMF算法的恢復性能不同，在ME數據集中，恢復誤差更小，這可能與該數據集的傳感器較少且相關性較強有關。對比算法中，SVP算法的重構精度較低，這是因為缺失率的提高會導致采樣數目的降低，從而使得矩陣填充算法的精確重構條件難以達到，因此誤差較大。未加正則項的MF算法RMSE顯然要高于本文算法，充分證明了融合先驗信息之后，對MF算法恢復數據的準確性提高起到了極大的作用。

4.3.2 算法運算時間對比

在3種數據集中的缺失率都設置為0.7的前提下，表2給出了本文算法與4種對比算法在上述3種數據集中100次計算過程的平均運行時間。仿真實驗環境為Inter Core i7, 2.30 GHz處理器,Windows10操作系統，仿真軟件為Matlab R2016 b。

表2 不同算法平均運行時間的比較

從表2中數據可以看出,未加正則項的MF算法的運行時間最短,DRMF算法的運行時間比MF算法長,但是明顯比GSDAM算法短,總體來說運算時間在一個合理的范圍內。

5 結 論

本文針對多源時間序列缺失數據恢復精度不高的問題,提出了一種基于雙重正則矩陣分解的方法,該方法在挖掘數據內部先驗信息的基礎上,利用時間序列的平滑性構建時間序列隱含因子的二階差分正則項,引入反映數據內部結構的圖拉普拉斯正則項對傳感器隱含因子進行約束,并在圖拉普拉斯矩陣獲取過程中設計了一種聯合歐式距離和D-PCC的相關性度量策略,最后將雙正則項統一于矩陣分解的框架中,從而獲得較高的恢復性能。但是復雜的先驗融入矩陣分解的結構之后,會導致計算復雜度隨著數據量的增加而上升,因此DRMF算法框架下的大規模數據的計算效率問題是下一步需要解決的問題。