巴麗偉,童常青
(杭州電子科技大學(xué)理學(xué)院,浙江 杭州 310018)
因子分析(Factor Analysis,F(xiàn)A)是對(duì)高維數(shù)據(jù)建模分析的重要方法,通過(guò)對(duì)高維觀測(cè)數(shù)據(jù)降維,用少數(shù)幾個(gè)不相關(guān)的綜合指標(biāo)來(lái)表示多個(gè)實(shí)測(cè)變量,通常用于行為科學(xué)、社會(huì)科學(xué)、市場(chǎng)營(yíng)銷、產(chǎn)品管理等領(lǐng)域,在機(jī)器學(xué)習(xí)領(lǐng)域應(yīng)用廣泛[1]。岳博等[2]使用期望最大化算法(Expectation Maximization,EM)進(jìn)行因子分析,運(yùn)用最大似然估計(jì)得到模型參數(shù)計(jì)量,但未考慮模型的復(fù)雜性,容易陷入局部最大值。為了克服這些問(wèn)題,Mayekaw[3]引入貝葉斯理論進(jìn)行因子分析,使用EM算法找到參數(shù)的點(diǎn)估計(jì)并作為后驗(yàn)分布的邊際分布,但無(wú)法證明所使用的EM算法是收斂的。實(shí)際模型中很難通過(guò)貝葉斯理論估計(jì)模型的參數(shù),無(wú)法得到精確解,一般采用近似解去近似精確解。Blei等[4]將近似解的求解方式分為隨機(jī)近似方法和確定性近似方法,一般情況下,確定性近似方法比隨機(jī)近似方法更容易收斂。變分貝葉斯推斷算法是一種確定性近似方法,類似于EM算法,在求解概率模型中有著廣泛的應(yīng)用。
Ghahramani等[5]給出了混合因子分析的變分貝葉斯推斷算法,從貝葉斯角度估計(jì)模型參數(shù),解決了聚類問(wèn)題和局部維數(shù)減少問(wèn)題。Wei等[1]從貝葉斯角度進(jìn)行混合公因子分析,提出公因子分析的貝葉斯混合模型,并結(jié)合變分貝葉斯推斷算法,對(duì)高維合成數(shù)據(jù)聚類,最后使用貝葉斯信息準(zhǔn)則(Bayes Information Criterion,BIC)進(jìn)行模型選擇。變分貝葉斯信息準(zhǔn)則(Variational Bayes Information Criterion,VBIC)[6]是BIC的變分貝葉斯近似,即使沒(méi)有明確定義數(shù)據(jù)集維度和個(gè)數(shù),仍然可以用VBIC進(jìn)行模型的選擇。但是,上述文獻(xiàn)中,直接為先驗(yàn)分布的參數(shù)設(shè)置一個(gè)超先驗(yàn)分布,并沒(méi)有說(shuō)明設(shè)置超先驗(yàn)分布時(shí)模型的參數(shù)估計(jì)效果。為了解決這個(gè)問(wèn)題,本文為數(shù)據(jù)樣本建立FA模型,并提出一種新的變分貝葉斯推斷算法,用于FA模型的參數(shù)估計(jì),重點(diǎn)討論因子分析模型中超先驗(yàn)分布對(duì)變分貝葉斯估計(jì)結(jié)果的影響。
yi=Axi+ξi,i=1,…,n
(1)
式中,xi為k(k
假設(shè)觀測(cè)值yi∈Rp是連續(xù)的,就像線性回歸一樣,設(shè)均值是(隱變量)輸入的線性函數(shù),從而似然函數(shù)為:
p(yi|A,D,xi)=N(yi|Axi,D)
(2)
在貝葉斯框架中,模型的參數(shù)被視為隨機(jī)變量,因此需要設(shè)置參數(shù)的先驗(yàn)分布。在共軛指數(shù)域內(nèi),為了保證先驗(yàn)分布和后驗(yàn)分布具有一致性,在變分貝葉斯推斷算法中將參數(shù)的先驗(yàn)分布設(shè)置為共軛先驗(yàn)分布,否則,很難求出變分后驗(yàn)分布的具體形式。詳細(xì)的選擇過(guò)程如下。
參考文獻(xiàn)[1]中對(duì)于因子載荷矩陣A的先驗(yàn)分布設(shè)置,假設(shè)A的每一列服從相同均值和不同精度的高斯分布,
(3)
式中,A.j表示矩陣A的第j列;v={v1,…,vk}是一個(gè)隨機(jī)變量,考慮到對(duì)A分布的影響,選擇在精度v的每一個(gè)元素上放置一個(gè)超先驗(yàn)分布:
(4)
式中,a*和b*分別為形狀參數(shù)和逆尺度參數(shù)。A的分布從行向量的角度進(jìn)行描述:
(5)
式中,Aq.表示矩陣A的第q行;diag(v)為對(duì)角線上元素為{v1,…,vk}的對(duì)角矩陣。
因?yàn)镈的自由度數(shù)不隨組件的維數(shù)增加而增加,所以,進(jìn)行參數(shù)的先驗(yàn)分布假設(shè)時(shí),沒(méi)有假設(shè)D的先驗(yàn)分布,而是把D當(dāng)成超參數(shù)直接進(jìn)行優(yōu)化。用概率圖形來(lái)表示FA的模型結(jié)構(gòu),如圖1所示。
圖1 FA概率圖形模型
(6)
對(duì)于采樣后的數(shù)據(jù)集,根據(jù)貝葉斯定理,參數(shù)和隱變量的后驗(yàn)分布表示為:
(7)
事實(shí)上,可以采用最大似然估計(jì)或最大后驗(yàn)估計(jì)等點(diǎn)估計(jì)的方式來(lái)求解式(7)的后驗(yàn)分布,但在此模型中,聯(lián)合概率分布函數(shù)比較復(fù)雜,采用上述方法很難求解,因此本文提出一種變分貝葉斯推斷算法,用一種簡(jiǎn)單的分布q(θ)來(lái)近似后驗(yàn)分布,假定Z是不可觀測(cè)變量的概率密度分布族,q(θ)∈Z。為了描述q(θ)與p(θ|y)的逼近程度,引入相對(duì)熵(Kullback-Leibler,KL),即:
(8)
依據(jù)泛函理論[8],證據(jù)下界L(q)關(guān)于q(θ)的泛函求導(dǎo),得到各參數(shù)的通解表達(dá)式:
lnq(θj)=Ei≠jlnp(y,θ)+c
(9)
式中,Ei≠j[lnp(y,θ)]為聯(lián)合概率密度函數(shù)的期望不涉及第j個(gè)變分因素,c表示不依賴于參數(shù)θj的常數(shù)。
假設(shè)因子載荷矩陣A的共軛先驗(yàn)分布中的參數(shù)v有超先驗(yàn)分布,此時(shí)對(duì)應(yīng)的模型稱為有超先驗(yàn)分布的FA模型。根據(jù)式(9)計(jì)算出有超先驗(yàn)分布的FA模型中各組參數(shù)的變分后驗(yàn)分布。
2.1.1 因子載荷矩陣A的第j列的精度vj的變分后驗(yàn)分布
由圖1可以看出,精度vj的更新可以利用因子載荷矩陣和參數(shù)的先驗(yàn)信息推導(dǎo)得出,由共軛性知,vj,?j∈[1,k]的變分后驗(yàn)是Gamma分布,即
lnGa(vj|a,bj)
具有形狀參數(shù)a和逆尺度參數(shù)bj,后驗(yàn)參數(shù)可利用下式更新:
(10)
其中Aqj表示因子載荷矩陣的第q行第j列,〈·〉表示對(duì)相關(guān)后驗(yàn)分布的期望。
2.1.2 q(A)的變分后驗(yàn)分布
A的更新可以利用參數(shù)的先驗(yàn)信息和可觀測(cè)數(shù)據(jù)的聯(lián)合估計(jì)信息推導(dǎo)得出,
其中tr(·)為跡函數(shù),q(A)可以寫(xiě)成行向量的后驗(yàn)分布的乘積,即:
參數(shù)的更新公式為:
(11)
2.1.3 隱變量xi的變分后驗(yàn)分布
隱變量xi的后驗(yàn)推斷利用(2)式似然函數(shù)的信息和xi的先驗(yàn)信息進(jìn)行推導(dǎo),
lnq(xi)=EAlnp(yi|A,xi,D)+lnp(xi)+c=
(12)
根據(jù)式(10)—式(12),計(jì)算得到參數(shù)更新公式中的期望值為:
2.1.4 模型的證據(jù)下界
根據(jù)參數(shù)的先驗(yàn)分布和變分后驗(yàn)分布信息,計(jì)算式(8)的證據(jù)下界L(q),
(13)
另外,相對(duì)于參數(shù)D,通過(guò)L(q)對(duì)D的求導(dǎo)來(lái)更新D,此對(duì)角矩陣的第q個(gè)元素Dqq為:
(14)
將式(10)—式(12)、式(14)的更新關(guān)系代入式(13),不斷迭代更新模型的參數(shù)直到證據(jù)下界收斂,當(dāng)L(qt)-L(qt-1)<ε時(shí),說(shuō)明L(q)達(dá)到收斂狀態(tài),迭代停止,否則,再進(jìn)行下一次迭代。其中L(qt)表示第t次迭代證據(jù)下界的值,ε表示閾值。L(q)的具體表達(dá)式在此不再展開(kāi)。
綜上所述,整個(gè)算法流程如下。
輸入:數(shù)據(jù)集y1︰n以及因子的個(gè)數(shù)k輸出:變分密度q(θ),變分下界L(q)初始化:隱變量xi的均值和方差,以及D,超參數(shù)a*,b*(1)通過(guò)式(10)—式(12)、式(14)更新超參數(shù);(2)通過(guò)計(jì)算得到證據(jù)下界L(q);(3)當(dāng)L(qt)-L(qt-1)<ε,L(q)達(dá)到收斂狀態(tài),迭代停止,否則,重復(fù)步驟1至步驟3;(4)輸出變分密度q(θ),變分下界L(q)。
參數(shù)的更新公式為:
(15)
變分下界變?yōu)椋?/p>
(16)
參數(shù)D只與似然函數(shù)有關(guān),因此D的更新形式?jīng)]有發(fā)生改變,如式(14)。
為了比較有超先驗(yàn)分布的FA模型和無(wú)超先驗(yàn)分布的FA模型中,變分貝葉斯推斷算法對(duì)參數(shù)估計(jì)的效果,可以采用一種標(biāo)準(zhǔn)進(jìn)行衡量。在早期的研究中,主要利用赤池信息準(zhǔn)則(Akaike Information Criterion,AIC)、貝葉斯信息準(zhǔn)則BIC、偏差信息準(zhǔn)則(Deviance Information Criterion,DIC)等比較模型的好壞[9]。本文中,主要采用變分貝葉斯信息準(zhǔn)則VBIC對(duì)模型進(jìn)行選擇,在不知道數(shù)據(jù)集維度和個(gè)數(shù)情況下,VBIC仍可以進(jìn)行模型的選擇。VBIC是BIC的變分貝葉斯近似,表達(dá)式定義為[6]:
VBIC=-2L(q)+2Eq[lnp(θ)]
根據(jù)式(13)、式(16)和參數(shù)的先驗(yàn)信息,得到有超先驗(yàn)分布的FA模型和無(wú)超先驗(yàn)分布的FA模型的VBIC分別為:
(17)
本文中,運(yùn)用VBIC比較有超先驗(yàn)分布的FA模型和無(wú)超先驗(yàn)分布的FA模型的優(yōu)劣,VBIC的值越小,說(shuō)明模型的參數(shù)估計(jì)效果越佳。
運(yùn)用式(1)生成一個(gè)高維數(shù)據(jù)集,且樣本容量是可變的,為了方便起見(jiàn),選擇特征個(gè)數(shù)p=10,公共因子的個(gè)數(shù)k=2,模型參數(shù)的真實(shí)值為:
D=diag(0.6,0.1,0.2,0.2,0.3,0.4,0.3,0.5,0.7,0.8)
為了避免在更新迭代過(guò)程中出現(xiàn)跟蹤停滯問(wèn)題,通常為共軛先驗(yàn)分布中的超參數(shù)設(shè)置一個(gè)初始值[10]。變分貝葉斯推斷算法對(duì)初始值不敏感,本文設(shè)定的初始值與真實(shí)值較近,目的是讓迭代步驟更少。將模型的參數(shù)初始化設(shè)置為:a*=10,b*=1,xi的均值設(shè)定為[-0.040,0.011],xi的方差設(shè)置為I,噪聲方差的對(duì)角線元素均設(shè)置為0.01,當(dāng)參數(shù)無(wú)超先驗(yàn)分布時(shí),v設(shè)置為[50,100]。采用本文提出的變分貝葉斯推斷算法進(jìn)行FA模型參數(shù)的估計(jì),并運(yùn)用VBIC進(jìn)行模型的選擇。
本文在樣本容量n=80和n=300下,分別對(duì)有超先驗(yàn)分布的FA模型和無(wú)超先驗(yàn)分布的FA模型的參數(shù)進(jìn)行變分貝葉斯估計(jì),由于參數(shù)個(gè)數(shù)較多,估算結(jié)果只給出參數(shù)D的估計(jì)及真實(shí)值和估計(jì)值之間的偏差絕對(duì)值,如表1所示。
表1 參數(shù)D的估計(jì)值和真實(shí)值的偏差
從表1可以看出,即使是在小樣本的情形下,運(yùn)用有超先驗(yàn)分布的變分貝葉斯推斷算法和無(wú)超先驗(yàn)分布的變分貝葉斯推斷算法對(duì)不同F(xiàn)A模型的參數(shù)進(jìn)行估計(jì),參數(shù)的估計(jì)值和真實(shí)值的偏差絕對(duì)值都很小,例如,超過(guò)50%的偏差絕對(duì)值小于0.050,大約有90%的偏差絕對(duì)值是小于0.100,這證明了變分貝葉斯推斷算法在FA模型的參數(shù)估計(jì)上的有效性。在小樣本和參數(shù)無(wú)超先驗(yàn)分布的情況下,參數(shù)估計(jì)量的精度不夠高,不同模型和樣本容量下,參數(shù)估計(jì)均存在一定誤差,這是因?yàn)樽兎重惾~斯算法是近似推斷而非精確推斷。而當(dāng)樣本容量較大以及參數(shù)有超先驗(yàn)分布時(shí),參數(shù)的估計(jì)值和真實(shí)值更逼近,說(shuō)明本文提出的變分貝葉斯推斷算法更適合較大的樣本容量和有超先驗(yàn)分布的FA模型,可以通過(guò)擴(kuò)大樣本的容量和設(shè)置參數(shù)的超先驗(yàn)分布來(lái)減少參數(shù)估計(jì)的誤差,提高參數(shù)估計(jì)精度。
變分下界的迭代收斂過(guò)程如圖2所示,迭代設(shè)置的閾值ε=10-3,迭代的最大次數(shù)設(shè)置為1 000。有超先驗(yàn)分布的變分貝葉斯推斷算法和無(wú)超先驗(yàn)分布的變分貝葉斯推斷算法的目標(biāo)函數(shù)的變化趨勢(shì)如圖2所示。從圖2可以看出,隨著迭代次數(shù)的增加,變分下界不斷增加,直到收斂。
圖2 變分下界的迭代收斂過(guò)程
在不同數(shù)據(jù)集下,對(duì)有超先驗(yàn)分布的FA模型和無(wú)超先驗(yàn)分布的FA模型的迭代次數(shù)、ELBO,VBIC進(jìn)行比較,結(jié)果如表2所示。
表2 不同數(shù)據(jù)集下,不同F(xiàn)A模型的衡量指標(biāo)比較
ELBO是邊際似然的下界,當(dāng)ELBO的值越高,說(shuō)明模型對(duì)數(shù)據(jù)擬合程度越好,模型產(chǎn)生數(shù)據(jù)的概率也越高。由表2可以看出,在不同樣本容量下,有超先驗(yàn)分布的FA模型的ELBO值更大一些,并且有超先驗(yàn)分布的FA模型的VBIC小于無(wú)超先驗(yàn)分布的FA模型。綜合來(lái)說(shuō),有超先驗(yàn)分布的FA模型對(duì)數(shù)據(jù)的擬合程度更好,因此,有超先驗(yàn)分布的FA模型是最好的選擇。
本文提出一種有關(guān)FA模型的變分貝葉斯推斷算法,重點(diǎn)討論FA模型中超先驗(yàn)分布對(duì)變分貝葉斯估計(jì)結(jié)果的影響,并采用變分貝葉斯信息準(zhǔn)則進(jìn)行模型的選擇。實(shí)驗(yàn)研究發(fā)現(xiàn),當(dāng)參數(shù)設(shè)置超先驗(yàn)分布時(shí),F(xiàn)A模型的參數(shù)估計(jì)效果更佳。后續(xù)計(jì)劃將本文算法應(yīng)用于結(jié)構(gòu)方程模型(Structural Equation Modeling,SEM)的參數(shù)估計(jì)和模型選擇中。