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

李小沛 李凡長(zhǎng) 梁合蘭

(蘇州大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院 江蘇 蘇州 215006)

0 引 言

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

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

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

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

1 相關(guān)工作

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

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

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

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

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

2.1 張量模型

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

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

(1)

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

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

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)

式中：符號(hào)⊙表示KR積；符號(hào)?表示Kronecker積。

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

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

圖1 張量CP分解模型

k=1,2,…,K

(4)

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

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

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

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

(5)

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

(6)

如果將A、B和C標(biāo)準(zhǔn)化，然后提出權(quán)重λ∈RR，可以得到如下表述：

(7)

對(duì)于N維度的情況：

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

(8)

(9)

(10)

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

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

(11)

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

X～N(μ,σ2)

(12)

當(dāng)μ=0,σ=1時(shí)，正態(tài)分布就成為標(biāo)準(zhǔn)正態(tài)分布，記為：

(13)

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

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

(14)

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

lnX～N(μ,σ2)

(15)

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

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

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

lnxi～N(μ,λ-1)

(16)

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

(17)

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

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

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

(μ(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)的隨機(jī)變量服從對(duì)數(shù)正態(tài)分布的假設(shè)下，方差用來捕獲數(shù)據(jù)中偏離預(yù)期值的數(shù)據(jù)。然而，方差在已知的數(shù)據(jù)集中對(duì)我們而言也是未知參數(shù)。因此對(duì)精度參數(shù)使用一個(gè)單變量高斯分布，其似然函數(shù)為：

(19)

其共軛Gamma先驗(yàn)分布公式為：

其后驗(yàn)分布為：

記為：

λ～Gamma(a0,b0)

(20)

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

(a) 矩陣

(b) 三階張量

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

2.4 Gibbs采樣

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

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

(21)

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

(22)

其中：

(23)

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

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

(24)

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

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

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

(25)

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

(26)

(27)

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

似然函數(shù)為：

(28)

結(jié)合式(28)中的似然項(xiàng)和式(20)的先驗(yàn)項(xiàng)，可以給出λ的后驗(yàn)分布公式：

(29)

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

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

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

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

設(shè)置隨機(jī)變量U(k)使其服從一個(gè)對(duì)數(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í) 驗(yàn)

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

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

3.1 對(duì)比實(shí)驗(yàn)

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

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

使用平均絕對(duì)百分比誤差(MAPE)和均方根誤差(RMSE)作為模型的評(píng)價(jià)指標(biāo)，其計(jì)算公式分別為：

(30)

(31)

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

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

3.3 實(shí)驗(yàn)結(jié)果

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

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

表1 實(shí)驗(yàn)結(jié)果

4 結(jié) 語

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

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