一種貝葉斯對數(shù)正態(tài)分布的張量分解插補算法

李小沛 李凡長 梁合蘭

(蘇州大學計算機科學與技術學院 江蘇 蘇州 215006)

0 引 言

智能交通系統(tǒng)已經(jīng)被用于在各種場景中獲得大量的城市交通數(shù)據(jù)。這些數(shù)據(jù)可以作為各種研究的數(shù)據(jù)集為交通管理者和研究者使用。通常可以將數(shù)據(jù)提取成矩陣存儲，考慮將矩陣擴展成高維度的數(shù)據(jù)，可以保存數(shù)據(jù)本身的信息結構，高維度的數(shù)據(jù)存儲形式被稱作張量。所以可以將取得的數(shù)據(jù)變成一個多維數(shù)組，如按照路段×天×一天的時間間隔，或者用戶×路段×某天。

雖然可以利用取得的數(shù)據(jù)進行整理分析預測，但是總會有各種外界因素，例如某地的網(wǎng)絡信號有所浮動等故障，這會對數(shù)據(jù)的獲取帶來影響，最直接的結果就是使取得的時空數(shù)據(jù)有所缺陷，即數(shù)據(jù)缺失問題。如果該道路在給定的時空范圍內(nèi)沒有任何數(shù)據(jù)信息，會使得數(shù)據(jù)集變得稀疏，導致無法準確地進行預測和估算，這會給基于這些數(shù)據(jù)集的各種工作帶來負擔。為了解決上述問題，就產(chǎn)生了對不完整數(shù)據(jù)集的缺失條目提供可靠估計的研究。

由于每個時空數(shù)據(jù)由不同含義的維度來定義，選擇用張量來處理這類數(shù)據(jù)。處理稀疏張量數(shù)據(jù)的方式就是張量分解。當前較為流行的張量分解為CANDECOMP/PARAFAC模型和Tucker模型[1]。在插補任務中，考慮對未知條目的可靠預測，引入了一種貝葉斯矩陣分解算法[2]，并將其拓展為高階張量來使用，即貝葉斯張量分解[3]。

這種完全的貝葉斯模型表征了數(shù)據(jù)生成過程，因此能夠有效地估算缺失的條目。利用貝葉斯分析，在推導過程中，將共軛先驗使用在超參數(shù)上，分析并得出其后驗分布，得到馬爾可夫鏈蒙特卡洛算法的采樣模型[2]。MCMC算法為缺失值提供了多種估計，避免單點估計帶來的較大誤差[4]。由于數(shù)據(jù)集可以很好地擬合成高斯分布而沒有極端值，考慮到速度數(shù)據(jù)集非負這一特性，選取一組服從對數(shù)正態(tài)分布的隨機數(shù)據(jù)作為初始值，并選用了如下三種數(shù)據(jù)的表示形式[10]：(1) 矩陣(路段×時間間隔)；(2) 三階張量(路段×天×時間間隔)；(3) 四階張量(路段×周×天×時間間隔)。本文在這三種表示上評估了不同的模型的插補性能，發(fā)現(xiàn)本模型在三階張量結構上的性能較好。

1 相關工作

張量已經(jīng)是統(tǒng)計學問題中最為常見的數(shù)據(jù)存儲形式。前人已經(jīng)對張量的知識及其分解應用做了非常全面的綜述，例如特征提取、數(shù)據(jù)填補、軌跡預測、數(shù)據(jù)降噪、信號處理及推薦系統(tǒng)方面等[6-8,15-16,24-25]。

已知的對缺失數(shù)據(jù)插補的方法中最常用的就是張量分解方法。利用觀察到的數(shù)據(jù)估計低秩近似的模型，并且利用估計的低秩近似模型計算出缺失的條目[9]。

張量已經(jīng)被應用在智能交通數(shù)據(jù)的研究中，例如在時空交通速度數(shù)據(jù)和軌跡預測數(shù)據(jù)集中應用張量來處理數(shù)據(jù)：文獻[9]提出了一種SVD組合張量分解恢復不完全數(shù)據(jù)的方法；文獻[10]提出了一種低秩張量分解模型(HaLRTC)來恢復從城市道路網(wǎng)收集的缺失交通速度數(shù)據(jù)集；文獻[11]提出了一種基于上下文感知的張量Tucker分解的方法對收集的軌跡數(shù)據(jù)進行預測；文獻[12]提出了一種基于上下文感知的將CP和Tucker分解結合的BTD方法對數(shù)據(jù)進行預測；文獻[13]提出了一種Tucker分解方法來促進核心張量的簡約性；文獻[14]提出了一種基于動態(tài)張量進行短期交通流量預測的方法；文獻[3]提出了一種貝葉斯高斯CP張量方法，基于貝葉斯推斷，利用服從標準正態(tài)分布的隨機變量迭代，對缺失數(shù)據(jù)集進行CP分解來進行數(shù)據(jù)集的插補。除此之外，還有很多運用張量的算法研究[15-17,21-22]，也有將張量和深度學習聯(lián)合應用的研究[18-19]。

先前的研究大多數(shù)采用單點計算，在逼近過程中就會產(chǎn)生較大的誤差，不能得出較為準確的結果。由于貝葉斯推斷只需要較少的觀察來進行估計，為了解決這個問題，將貝葉斯推理方法加入到張量分解問題當中，這可以有效地解決上述問題[16,20]。

2 貝葉斯對數(shù)正態(tài)分布張量分解

2.1 張量模型

幾何代數(shù)中定義的張量是基于向量和矩陣的推廣。簡單而言，可以將標量視為零階張量，多個標量組成的矢量視為一階張量，多個矢量組成的矩陣視為二階張量，多個矩陣的排列可以視為三階張量，之后增加了維度的數(shù)據(jù)結構統(tǒng)稱為高階張量，空間上更為立體。為了便于張量運算，介紹以下張量分解中的數(shù)學知識。

(1) Kronecker積。Kronecker積在張量計算中十分常見，它是銜接矩陣計算和張量計算的橋梁。Kronecker積的定義為：給定一個大小為m1×m2的矩陣A，一個大小為n1×n2的矩陣B，矩陣A和B的Kronecker積為：

(1)

式中：矩陣A?B的大小為(m1n1)×(m2n2)；符號?表示Kronecker積。

(2) Khatri-Rao積。在矩陣分解過程中，利用一種特殊的矩陣運算加快采樣過程，這個運算就是Khatri-Rao積(KR乘積)。其定義為：給定如下兩個矩陣：

A=(a1,a2,…,ar)∈Rm×r

B=(b1,b2,…,br)∈Rn×r

(2)

矩陣A、B有相同的列數(shù)，則其KR積為：

A⊙B=(a1?b1,a2?b2,…,ar?br)∈R(mn)×r

(3)

式中：符號⊙表示KR積；符號?表示Kronecker積。

(3) 張量的秩。張量的秩與矩陣的秩的定義類似，但是有一些不同。張量的秩在實數(shù)域和復數(shù)域可以不同，并且沒有直接的算法可以計算給定的張量的秩。通常情況下都是選擇多個不同的值。本文分別選擇了三個不同的值進行比較：r=50,80,110。

(4) CANDECOMP/PARAFAC(CP)張量分解。CP分解是將高維度的張量分解成為秩一張量的和，如圖1所示。

圖1 張量CP分解模型

k=1,2,…,K

(4)

定義因子矩陣表示向量的組合，例如A=[a1,a2,…,aR]，然后對三維情況，可以寫成如下形式：

X(1)≈A(C⊙B)T

X(2)≈B(C⊙A)T

X(3)≈C(B⊙A)T

(5)

式中:X(i)表示張量X在模式i上的展開。CP分解可以簡明地表述為：

(6)

如果將A、B和C標準化，然后提出權重λ∈RR，可以得到如下表述：

(7)

對于N維度的情況：

x≈[[λ;A(1),A(2),…,A(N)]]≡

(8)

(9)

(10)

2.2 對數(shù)正態(tài)分布

正態(tài)分布，又稱為高斯分布，若隨機變量X服從一個位置參數(shù)為μ、尺度參數(shù)為σ的概率分布，且其概率密度函數(shù)為：

(11)

這個隨機變量就可稱為正態(tài)隨機變量，正態(tài)隨機變量服從的分布就稱為正態(tài)分布，記作：

X～N(μ,σ2)

(12)

當μ=0,σ=1時，正態(tài)分布就成為標準正態(tài)分布，記為：

(13)

但是由于有些量的值必須是正數(shù)，由于高斯分布中X中的值可正可負，所以它不能很好地描述這樣的變量分布。在很多應用中，可能數(shù)據(jù)不符合正態(tài)分布，但是隨機變量的對數(shù)可能符合正態(tài)分布，此情況就稱之為對數(shù)正態(tài)分布，它是指一個隨機變量的對數(shù)服從正態(tài)分布，則該隨機變量服從對數(shù)正態(tài)分布。

若X是取值為正數(shù)的隨機變量，若lnX～N(μ,σ2)，且其概率密度函數(shù)為：

(14)

則稱隨機變量X服從對數(shù)正態(tài)分布，記為：

lnX～N(μ,σ2)

(15)

2.3 貝葉斯對數(shù)正態(tài)分布CP分解

對多維交通數(shù)據(jù)變成的張量進行插補缺失條目的操作，利用了貝葉斯張量分解，貝葉斯分析方法是將未知參數(shù)的先驗信息，與樣本信息聯(lián)合，根據(jù)貝葉斯公式，得出后驗分布公式，而后依照后驗信息推斷未知參數(shù)。由于本文算法選用的數(shù)據(jù)集是速度數(shù)據(jù)集，這要求每一個值都是非負的，所以選用服從對數(shù)正態(tài)分布的隨機變量作為初始值進行算法循環(huán)。

首先，假設觀察到的項遵循獨立的對數(shù)正態(tài)分布：

lnxi～N(μ,λ-1)

(16)

式中：μ代表是均值；λ代表精度，λ-1代表方差。假設所有的因子矩陣中的行向量的先驗分布為多變量的正態(tài)分布：

(17)

多變量正態(tài)分布是單維正態(tài)分布向多維的推廣。在貝葉斯環(huán)境中，使用多變量共軛Gaussian-Wishart先驗來增強模型的魯棒性，在μ和Λ都是未知的情況下，Gaussian-Wishart分布公式為：

p(μ,Λ|μ0,β0,W0,v0)=N(μ|μ0,(β0Λ)-1)W(Λ|W0,v0)

將超參數(shù)μ(k)和Λ(k)代入該公式，結果為：

(μ(k),Λ(k))～Gaussian-Wishart(μ0,β0,W0,v0)

p(μ(k),Λ(k)|μ0,β0,W0,v0)=

N(μ(k)|μ0,(β0Λ(k))-1)×Wishart(Λ(k)|W0,v0)

(18)

在式(16)的隨機變量服從對數(shù)正態(tài)分布的假設下，方差用來捕獲數(shù)據(jù)中偏離預期值的數(shù)據(jù)。然而，方差在已知的數(shù)據(jù)集中對我們而言也是未知參數(shù)。因此對精度參數(shù)使用一個單變量高斯分布，其似然函數(shù)為：

(19)

其共軛Gamma先驗分布公式為：

其后驗分布為：

記為：

λ～Gamma(a0,b0)

(20)

圖2引用自文獻[3]，展示了上述貝葉斯概率CP分解的數(shù)據(jù)生成過程的結構圖。三個模型圖分別對應矩陣、三階張量和四階張量。在此模型圖中，節(jié)點xi代表觀測到的數(shù)據(jù)；λ是參數(shù)；μ(k)、Λ(k)、a0和b0都是需要預先設置的超參數(shù)。當初始化采樣時，應設置這些超參數(shù)。這樣可以基于此模型推導出用于貝葉斯推斷的Gibbs采樣算法。

(a) 矩陣

(b) 三階張量

(c) 四階張量圖2 貝葉斯概率CP分解的數(shù)據(jù)生成過程的結構圖

2.4 Gibbs采樣

不斷循環(huán)使得隨機變量的概率分布將收斂于x的聯(lián)合概率分布p(x)。

Gibbs采樣的思想是在每次迭代中順序更新所有變量。在所有其他變量的當前值的條件下，從其分布中采樣每個變量。Gibbs采樣算法的關鍵是為所有變量定義這種分布，這些條件分布通常稱為完全條件。由于已經(jīng)將共軛先驗用于參數(shù)和超參數(shù)，因此圖中模型的后驗分布可以通過分析得出，以圖2中三階張量為例進行推導。

(21)

式中：?代表Hadamard乘積。結合式(17)和式(21)，給出后驗分布多變量高斯的形式：

(22)

其中：

(23)

2) 采樣超參數(shù)μ(k)和Λ(k)(k=1,2,3)。因子矩陣U(1)∈Rn1×R的似然函數(shù)為：

p(U(1)|μ(1),Λ(1))∝

(24)

結合式(24)和式(18)中的Gaussian-Wishart超先驗項，超參數(shù)μ(1)和Λ(1)的聯(lián)合后驗分布為：

p(μ(1),Λ(1)|U(1),μ0,W0,v0,β0)∝

p(U(1)|μ(1),Λ(1))×N(μ(1)|μ0,(β0Λ(1))-1)×

(25)

這兩個分布中的參數(shù)可以通過以下方式計算：

(26)

(27)

3) 采樣方差的倒數(shù)λ。

似然函數(shù)為：

(28)

結合式(28)中的似然項和式(20)的先驗項，可以給出λ的后驗分布公式：

(29)

到此已經(jīng)推導出三階張量的Gibbs采樣中所需要的所有變量的公式，可同樣拓展至高維數(shù)據(jù)。Gibbs采樣算法達到收斂后，就可以通過Monte Carlo近似估計缺失值。

通過上述對各項參數(shù)的推導和后驗分布迭代公式，得出貝葉斯對數(shù)正態(tài)分布張量分解插補算法BLCP，如算法1所示。

算法1貝葉斯對數(shù)正態(tài)分布張量分解插補算法

初始化：μ0=0,β0=1,v0=r,W0=I其中I是單位矩陣，r是設置好的張量秩low rank。

設置隨機變量U(k)使其服從一個對數(shù)正態(tài)分布作為初始值循環(huán)：

Foriter=1,2，…,T

Fork=1,2，…,d

Forik=1,2,…,nk

End for

End for

End for

//T為迭代次數(shù)

3 實 驗

本文使用的交通速度數(shù)據(jù)集是通過智能手機上廣泛應用的導航程序生成的選取了文獻[3]中給出的在中國廣州收集的大規(guī)模交通速度數(shù)據(jù)集。數(shù)據(jù)集包含61天(2016年8月1日—2016年9月31日)，一天被劃分為144個間隔(以10分鐘為間隔)，并且記錄了214個路段上行駛速度的觀測值[3]。將這個數(shù)據(jù)集分別變成二階、三階和四階張量進行實驗。其中觀測到的數(shù)據(jù)大約有1 855 527個，此數(shù)據(jù)集有約1.29%的數(shù)據(jù)缺失值。圖3顯示了觀測到的速度數(shù)據(jù)的直方圖，可以看出，數(shù)據(jù)的分布可以擬合成一個最大似然的高斯分布而沒有極端值。

圖3 觀測到的速度數(shù)據(jù)直方圖和正態(tài)分布的擬合

3.1 對比實驗

通過隨機刪除一定數(shù)量的已知值來創(chuàng)建實驗數(shù)據(jù)集，然后將原始數(shù)據(jù)集作為缺失數(shù)據(jù)的基礎事實，可以利用這兩個數(shù)據(jù)集評估本文算法插補的性能。

在對比實驗中，對比了以下三個基于張量分解的插補方法和一個傳統(tǒng)的方法：(1) K最近鄰分類算法(KNN)；(2) 高精度低秩張量算法(HaLRTC)[10]；(3) SVD組合張量分解算法(STD)[9]；(4) 基于標準正態(tài)分布的貝葉斯張量分解方法(BGCP)[3]。

使用平均絕對百分比誤差(MAPE)和均方根誤差(RMSE)作為模型的評價指標，其計算公式分別為：

(30)

(31)

3.2 數(shù)據(jù)表示

張量的結構表示也是需要考慮的關鍵因素。基于交通速度數(shù)據(jù)集，本文使用了三種不同類型的數(shù)據(jù)表示形式[3]。

3.3 實驗結果

對于隨機丟失的情況，隨機刪除矩陣表示形式(1)中的某些條目，并將矩陣中的數(shù)據(jù)重塑為三階張量和四階張量的存儲形式。

隨機刪除數(shù)據(jù)在不同模型下的性能如表1所示。由于STD模型主要解決非凸性問題，因此它在矩陣表示的數(shù)據(jù)集和四階張量表示的數(shù)據(jù)集中沒有提供太多的改進，所以只在三階張量表示的數(shù)據(jù)集中有所應用。分別刪除了原始數(shù)據(jù)集的10%、20%、30%、40%、50%，整理出五個數(shù)據(jù)集。本文模型對Gibbs采樣算法進行了1 000次迭代估算缺失條目(T=1 000)。分別選取了張量秩為r=50、80、110來運行所有算法模型，最優(yōu)結果已在表1中標粗。可以看出，BLCP在處理三階張量的缺失數(shù)據(jù)集上有較好的性能。

表1 實驗結果

4 結 語

本文提供了一種貝葉斯對數(shù)正態(tài)分布的張量CP分解模型，對時空數(shù)據(jù)集的缺失位置進行插補任務。并利用了收集到的時空交通速度數(shù)據(jù)集作為實驗數(shù)據(jù)集[3]，將數(shù)據(jù)集變成二階張量(矩陣)、三階張量和四階張量三種不同的維度作為實驗數(shù)據(jù)存儲形式。通過隨機刪除二階張量中的10%到50%的數(shù)據(jù)，并依次將其重塑成三階張量和四階張量作為實驗數(shù)據(jù)，以原始未刪除數(shù)據(jù)的數(shù)據(jù)集作為事實依據(jù)，并與HaLRTC、傳統(tǒng)的KNN算法、STD算法及基于標準正態(tài)分布的貝葉斯高斯張量分解算法(BGCP)進行了比較，結果表明，在隨機刪除數(shù)據(jù)情況下，本文算法在處理三階張量存儲形式的數(shù)據(jù)上能獲得較好的結果。

由于張量分解方法都是基于對低秩張量的逼近假設，并且時空數(shù)據(jù)集可能存在著強烈的時空關聯(lián)性，在之后的研究可以考慮將時空關聯(lián)的影響加入到算法中，并且由于CP分解張量的秩難以求出并且取得一個較好的值，可以考慮非參數(shù)的貝葉斯模型應用于CP分解上。除此之外，由于張量CP分解可能存在大量的冗余因子，可以考慮開發(fā)基于貝葉斯的Tucker算法解決上述問題。