最大期望模擬退火的貝葉斯變分推理算法

劉浩然 張力悅 蘇昭玉 張 赟 張 磊

①(燕山大學(xué)信息科學(xué)與工程學(xué)院 秦皇島 066004)

②(河北省特種光纖與光纖傳感重點實驗室 秦皇島 066004)

③(北京市機電研究院 北京 100027)

1 引言

在機器學(xué)習(xí)中，貝葉斯推理已經(jīng)成為求解不可觀測變量后驗概率的重要方法，它類似最大期望算法。當模型更為復(fù)雜時，貝葉斯精確推理求解具體的參數(shù)分布時間開銷巨大，而許多時候?qū)τ趶?fù)雜模型的參數(shù)分布精確求解不是必要的，而推理的真正目的是得到參數(shù)的期望或者近似分布，以了解數(shù)據(jù)遵循分布的規(guī)律[1,2]。精確的貝葉斯推理表現(xiàn)出計算量隨變量增加呈指數(shù)增長以及計算時間也迅速增加的問題，近似貝葉斯推理針對于此類問題有較好的方法，其中基于蒙特卡羅抽樣(Markov Chain Monte Carlo, MCMC)[3]的隨機近似算法通過采樣的方式計算參數(shù)后驗概率，MCMC算法提供了從目標分布漸進地生成精確樣本的精確率保證，然而MCMC方法在數(shù)據(jù)量較小時，無法保證采集過程的準確率，算法的結(jié)果較差，算法會隨著數(shù)據(jù)量和采樣量的增大，其結(jié)果越來越好，但與此同時計算消耗和時間消耗也越來越大，而另一種確定性近似算法即變分推理算法[4]將求解隱參數(shù)后驗問題轉(zhuǎn)化成變分優(yōu)化問題，通過迭代找到較優(yōu)的后驗分布解，這種方法在應(yīng)對更為復(fù)雜的模型下，表現(xiàn)出良好的計算效率[5]。

變分推理算法廣泛應(yīng)用于計算機科學(xué)、模式識別、圖像處理、卡爾曼濾波、戰(zhàn)場決策等領(lǐng)域[6–9]。近年來，許多學(xué)者深入研究變分推理算法并提出了有效的改進方案，以解決其存在的相關(guān)問題。Gianniotis等人[10]使用梯度下降優(yōu)化變分參數(shù)實現(xiàn)模型似然下界重新構(gòu)建，通過循環(huán)迭代得到模型參數(shù)后驗分布，算法具有較強的泛化能力，可應(yīng)用于多種場景，如貝葉斯線性回歸、貝葉斯多目標分類、概率圖像去噪等，但算法在引入梯度下降優(yōu)化似然下界時加重了原算法的局部最優(yōu)問題。Shekaramiz等人[11]提出使用貪婪策略篩選出支持度子集進行變分推理，該算法能有效地防止算法過擬合，但算法主要針對稀疏的貝葉斯學(xué)習(xí)，在數(shù)據(jù)樣本屬性較多時，該算法無法得到較好的模型參數(shù)后驗結(jié)果。Katahira等人[12]將最大熵引入確定性退火算法控制退火過程，算法得到相對原變分算法更優(yōu)的結(jié)果，但其隨機初始先驗對最終結(jié)果影響較大。Tabushi等人[13]使用非線性最大化非廣延統(tǒng)計力學(xué)的Tsallis熵提出了廣義確定性退火最大期望(Generalize Deterministic Annealing Expectation-Maximization, GDAEM)算法，通過控制參數(shù)得到全局較優(yōu)解，由于控制參數(shù)設(shè)置的問題，收斂效率有待提升。Salimans等人[14]提出了基于馬爾科夫鏈的變分近似，該算法將兩種近似方式有效結(jié)合，實現(xiàn)了速度和精度的相對平衡，但算法針對復(fù)雜的模型時，依然表現(xiàn)出采樣時間消耗大的問題。

本文針對部分算法對求解模型參數(shù)后驗分布時間消耗長、收斂精度低的問題，提出一種基于最大期望模擬退火的貝葉斯變分推理算法(Expectationmaximization and Simulated annealing for Variational Bayesian Inference, ES-VBI)，將雙重最大期望(Expectation-Maximization, EM)策略應(yīng)用于初始先驗的生成，將模擬退火策略的逆溫度參數(shù)用于似然下界的優(yōu)化，最終返回算法迭代最優(yōu)解，利用收斂性準則理論分析算法的收斂性，將該算法應(yīng)用于混合高斯分布實例的實驗仿真，說明本算法的優(yōu)勢。

2 問題描述

變分推理經(jīng)常被用于計算模型的參數(shù)后驗分布問題，通過變分法將求解模型參數(shù)后驗問題轉(zhuǎn)換為變分自由能最大化迭代尋優(yōu)問題，最終返回全局較優(yōu)解，即為最合理的近似模型參數(shù)后驗分布。本文提出的算法針對初始先驗和變分自由能在尋優(yōu)過程中進行優(yōu)化：(1)在保證算法速率的情況下，滿足模型參數(shù)的后驗分布的近似度最高；(2) 盡量降低算法對初始先驗的敏感性，保證算法具有對全局有效的初始先驗，防止算法過早地陷入局部最優(yōu)。針對變分推理過程中的上述特點，給出如下推理模型。

2.1 參數(shù)描述

本文給出3個定義描述推理模型。

(1)觀測變量和不可觀測變量。貝葉斯網(wǎng)絡(luò)中，變量類型分為觀測變量和不可觀測變量，其中觀測變量是直接可觀察或可采集的變量；不可觀測變量是不可直接觀察或不可直接采集的變量，不可觀測變量包括潛變量和模型參數(shù)，潛變量用于解釋觀測變量，可以看作其對應(yīng)觀測變量的抽象和概括。模型參數(shù)是描述模型數(shù)據(jù)自身存在規(guī)律的參數(shù)。

(2)庫爾貝克距離(Kullback-Leibler, KL)。衡量假設(shè)的近似分布q與真實分布p之間差異的量稱為KL散度，也稱KL距離[15]。其表達式如式(1)所示

其中，KL(q‖p)表示近似分布q與真實分布p的KL距離；W表示模型參數(shù)集合(隱變量集合)；X表示觀測變量數(shù)據(jù)。

在變分推理中，通過隨機賦值或者其他賦值方法給出一個初始先驗分布q，假設(shè)q是精確條件的候選假設(shè)分布，算法將KL(q‖p)不斷優(yōu)化達到最小，達到假設(shè)分布與真實分布的差異無限接近于0，得到對應(yīng)的q為

(3)證據(jù)下界(Evidence Lower BOund, ELBO)。觀測變量數(shù)據(jù)似然對數(shù)lgp(X)與KL距離的差值稱為證據(jù)下界(ELBO)[16]。在貝葉斯變分推理過程中，初始計算參數(shù)后驗分布時，根據(jù)貝葉斯公式，將數(shù)據(jù)分布表示為似然對數(shù)的形式lgp(X)，經(jīng)過化簡得到ELBO與lgp(X), KL距離的關(guān)系表達式

其中，L(q)稱為證據(jù)下界，即ELBO。此時lgp(X)是常數(shù)并保持不變，要使KL距離最小，則使ELBO達到最大，KL距離與ELBO的關(guān)系如圖1所示。由于在變分推理過程中，最小化KL距離無法通過式(1)直接計算，故通過最大化ELBO(即最大化L(q))，最終得到最優(yōu)的近似分布q?＝arg maxL(q)。

圖1 KL距離與ELBO的關(guān)系

2.2 模型后驗分布描述

設(shè)推理模型的觀測變量集合為X＝{x1,x2,···,xn}，其中xi,(i＝1,2,···,n)表示第i個觀測變量，n表示觀測變量的總數(shù)；設(shè)不可觀測變量集合為模型參數(shù)ω和潛變量集合z＝{z1,z2,···,zm}，其中模型參數(shù)ω表示描述模型數(shù)據(jù)在統(tǒng)計模型下的自身存在的規(guī)律，理論上最優(yōu)的模型參數(shù)與數(shù)據(jù)無限接近吻合(擬合度無限接近1 0 0%)，zj,(j＝1,2,···,m)表示第j個潛變量，m表示潛變量的總數(shù)。通過對聯(lián)合分布p(z,ω,X)積分得到觀測數(shù)據(jù)的邊緣密度分布(邊際似然)p(X)，如式(4)所示

根據(jù)貝葉斯定理，z和ω在X上的后驗概率表示為

將式(4)代入式(5)，式(5)分母部分表示成邊緣密度分布形式

將邊際似然等式(式(4))的兩邊進行對數(shù)運算，式(4)進一步分解為

根據(jù)Jeason不等式的性質(zhì)[11]，式(8)等式右邊部分存在如式(9)的不等關(guān)系

將式(9)中不等式右邊部分進一步分解為

2.3 雙重EM模型

當前多數(shù)的推理算法模型參數(shù)初始化都是在一定范圍內(nèi)進行隨機賦值操作，但較差的初始化結(jié)果，導(dǎo)致算法存在無法收斂到較優(yōu)值的情況，因此模型參數(shù)初始化對初始先驗敏感，在增加模型參數(shù)的數(shù)量時，這個問題顯得更加突出。針對先驗初始化問題，本算法采用雙重最大期望(EM)[19,20]方法，算法在第1次EM算法結(jié)果的基礎(chǔ)上運行第2次EM算法。在每次循環(huán)迭代，在第1階段，對一組隨機初始化數(shù)據(jù)使用EM算法，得到了第1階段的模型參數(shù)。在第2階段，針對第1階段結(jié)果再進行第2階段EM操作，完成ES-VBI的初始化。

求解含觀測變量、模型參數(shù)的條件概率最大時的模型參數(shù)ω,z，以求解ω為例，如式(12)所示

根據(jù)對數(shù)函數(shù)的單調(diào)性可得式(12)的等價形式為

算法將最大似然估計準則應(yīng)用到第1階段所提供的模型分布。在執(zhí)行該步驟過程中，選取來源于第1階段模型分布的數(shù)據(jù)樣本作為第2階段的初始先驗，該階段的參數(shù)分布后驗概率為

2.4 模擬退火模型

統(tǒng)計物理學(xué)的基本公式有L＝E-TS,(L是變分自由能，E是總能量，T是溫度，S是熵)[21]，根據(jù)此公式原理，本文在原目標函數(shù)(式(11))的基礎(chǔ)上引入逆溫度參數(shù)?＝1/T構(gòu)建新的目標函數(shù)，則此時變分目標函數(shù)表示為

逆溫度參數(shù)有兩種特殊情況：如果?→0，則第2項成為整個式子主導(dǎo)，使整個式子最大化相當于使第2個式子最大化，因均勻分布的熵最大，于是整個式子趨于一個均勻分布。如果?＝1，整個式子與式(11)相同，可知得到原始的后驗分布，模擬退火對算法的控制逐漸降低。

獨立隨機變量的聯(lián)合熵等于各變量熵的和，因此式(19)最右邊一項變?yōu)槭?20)最右邊兩項，式(20)為

3 算法實現(xiàn)

3.1 算法執(zhí)行過程

貝葉斯變分(VB)推理的核心是通過構(gòu)造一個簡單分布q去近似待求解的復(fù)雜分布p，通過不斷縮小它們之間的差異即KL距離，不斷增大變分自由能(ELBO)，使分布q與p無限接近，最終求得分布q即為最終的參數(shù)分布。

對以式(22)和式(23)構(gòu)造的拉格朗日式求導(dǎo)并令其等于0(其中F1(z,ω,X,λ1)對qz(z)求導(dǎo)，F(xiàn)2(z,ω,X,λ2)對qω(ω)求導(dǎo))，得

在原變分推理算法流程中加入雙重EM和逆溫度參數(shù)構(gòu)建退火循環(huán)，可以推導(dǎo)出如表1所示的ESVBI算法。

表1 ES-VBI算法

其中，關(guān)于模擬退火的迭代參數(shù)const 設(shè)置，退火過程是根據(jù)控制該參數(shù)的取值使得迭代的逆溫度參數(shù)?逐漸從一個小正數(shù)(正向趨于0小于1的數(shù))增加至1，以達到最大終止迭代條件，而且逆溫度參數(shù)初值也是正向趨于0小于1，所以const的初值應(yīng)設(shè)為大于1，這樣? ←?×const才能逐漸從0趨于1。本算法作為DA-VB和SA-VB算法的變體，對于該參數(shù)的設(shè)置結(jié)合了退火算法的整體性能以及大體參照DA-VB的參數(shù)const設(shè)置 1.1和1.2[12]，通過實驗發(fā)現(xiàn)兩種設(shè)置情況下均能達到最終的退火性能，因此本文選擇使用const＝1.1。

3.2 算法收斂性分析

定理1 ES-VBI算法中每次迭代增加L(q)，即L(qt+1)≥L(qt)，當且僅當?(qt+1|qt)＝?(qt)時，等號成立。

4 實驗仿真

4.1 高斯混合實例和實驗設(shè)置

高斯模型混合模型是比較經(jīng)典的包含隱參數(shù)的數(shù)據(jù)模型，我們利用提出的算法從混合高斯模型的數(shù)據(jù)近似擬合出混合高斯模型中的相關(guān)參數(shù)的估計值，而且對于在高斯混合模型可以較為準確地衡量算法對模型的魯棒性和準確性，具有較強的說服力，文獻[13]使用了高斯混合模型，文獻[14]雙變量高斯模型，而且提到可以推廣到多變量的混合模型。除此之外在關(guān)于變分推理算法對比文獻[10–12]中還使用了隱馬爾可夫模型、多目標分類、線性回歸、邏輯回歸等模型進行實驗仿真，對比算法的性能，提升算法的可信度。本算法與對比算法屬于同類型的推理算法可以針對上述模型進行實驗對比，本文由于文章內(nèi)容限制，僅在高斯混合模型下進行仿真實驗，具體如下：本文給出一個高斯混合模型對提出的ES-VBI算法性能進行實驗驗證，其中該混合模型包含K個高斯分布，以及n個可觀測變量xi,i＝1,2,···,n，潛在變量包含K個高斯分布的均值μk,k＝1,2,···,K和n個類別參數(shù)ci,i＝1,2,···,n，該模型表示為

為檢驗ES-VBI算法的迭代效率、收斂精度，將本文提出的ES-VBI算法與A-VI算法(原文中未給出算法具體名稱，本文記作A-VI(Advanced-Variational Inference)算法)[10]、OSBL-VB(Ordinary Sparse Bayesiam Learning-Variational Bayesian)算法[11]、DA-VB(Deterministic Annealing-Variational Bayesian)算法[12]、GDAEM(Generalized Deterministic Annealing Expectation-Max-imization)算法[13]、MCVI(Markov Chain Variational Inference)算法[14]在混合模型中進行仿真對比，另外本文使用了一些常用的變分推理的評價量：ELOB和時間t(s)。每個算法在相同數(shù)據(jù)集下獨立運行100次求取平均值為最后的統(tǒng)計結(jié)果。實驗環(huán)境：處理器Intel(R) Core(TM),CPU i7-7700HQ，主頻2.80 GHz，內(nèi)存為16 GB，Windows10 64 bit操作系統(tǒng)，python3.8版本。

4.2 各個算法針對高斯混合模型的實驗對比

在混合高斯模型優(yōu)化中生成器數(shù)據(jù)為1000時各個算法的ELOB(表示變分自由能，越大越好)和時間t(s)的對比結(jié)果如表2所示，它們的迭代過程如圖2所示。

圖2 各算法迭代過程對比圖

由表2數(shù)據(jù)可得，在6種算法中，ES-VBI算法的變分自由能E L O B 最大，且時間消耗僅比GDAEM算法多一些，但ES-VBI算法的收斂精度高于GDAEM算法，這是由于ES-VBI算法使用了雙重EM相對于GDAEM算法的單層EM時間消耗有所增加，但ES-VBI通過雙重EM的時間增加代價換來了算法的收斂精度提升。按照ELOB的收斂結(jié)果可以將6種算法的收斂精度從大到小排序為ESVBI算法、OSBL-VB算法、A-VI算法、GDAEM算法、MCVI算法、DA-VB算法。而且可以看出ES-VBI算法對比于確定性退火算法DA-VB收斂精度提升較大，說明基于逆溫度參數(shù)對算法的收斂進行有效的控制，使得算法在迭代尋優(yōu)過程的目標函數(shù)值ELOB更大，在每次迭代中更能優(yōu)化出較好的模型參數(shù)。

表2 各算法ELOB和時間對比

圖2顯示各個算法在1至100代的目標函數(shù)優(yōu)化圖，可以看出ES-VBI算法的初始化過程的結(jié)果顯然優(yōu)于其他的隨機初始過程結(jié)果，通過雙重EM模型初始化模型參數(shù)使得算法最終的收斂結(jié)果較優(yōu)。其中A-VI算法和OSBL-VB算法由于分別使用了梯度優(yōu)化和貪婪策略使得其相對另外3種算法最終的結(jié)果較好，通過圖2中A-VI的迭代線可以看出其在35代以后很難擺脫梯度優(yōu)化對局部最優(yōu)的加強影響，使得后期的迭代線較平，算法后期的優(yōu)化效率較低，且結(jié)果精度有待提升。

將高斯混合模型的K設(shè)置為5，通過數(shù)據(jù)生成器生成5000數(shù)據(jù)樣本，對比各個算法對高斯混合模型各個分量的擬合程度，其對比圖如圖3所示，圖3中每種顏色代表一個高斯分量。

圖3 各算法針對高斯混合模型擬合圖

由各個算法對應(yīng)的高斯混合模型擬合圖，可以看出本文提出的ES-VBI算法擬合程度最好，很且明顯優(yōu)于DA-VB算法和MCVI算法，DA-VB算法通過最大熵控制退火過程，其并不能很好地適應(yīng)算法迭代的過程，且隨機初始先驗也影響了最終的收斂結(jié)果。而本算法通過模擬退火的逆溫度參數(shù)有效地控制迭代的過程，相比于DA-VB算法的確定性退火原理在前期和后期都有對目標函數(shù)相對較優(yōu)的尋優(yōu)效率，使得算法最終的收斂精度較高，擬合混合高斯分布的各個分量較優(yōu)。

雖然本算法設(shè)計了一個改進的變分推理算法，算法對初始化的敏感性降低和迭代過程的有效控制，但所提出的方案仍然沒有得到全局最優(yōu)，只是在模型參數(shù)求解方面提供了一個優(yōu)化的方法，且上述實驗驗證了提出算法有效性和較優(yōu)性能。

5 結(jié)束語

針對貝葉斯變分推理的初始先驗和變分自由能優(yōu)化兩個問題，本文基于模擬退火理論和最大期望理論提出了一種針對降低初始先驗影響和變分自由能優(yōu)化的變分推理方法。通過雙重EM模型的構(gòu)建，使得算法保證近似精度的情況下，進一步提升算法在前期的參數(shù)分布質(zhì)量，通過模擬退火的逆溫度參數(shù)控制迭代優(yōu)化過程，進一步提升求解后驗分布的精度。通過收斂性理論證明了該推理模型收斂于全局較優(yōu)解。經(jīng)過實驗仿真證明本文的雙重EM模型和模擬退火的逆溫度參數(shù)的有效性，針對高斯混合模型本文提出算法表現(xiàn)出更好的性能，擬合模型參數(shù)的效果較優(yōu)，同時該算法能夠為卡爾曼濾波、貝葉斯神經(jīng)網(wǎng)絡(luò)參數(shù)求解、圖像去噪等提供理論支持。