魯棒概率矩陣三分解*

史加榮，陳姣姣

1.省部共建西部綠色建筑國家重點(diǎn)實驗室/西安建筑科技大學(xué)，西安710055

2.西安建筑科技大學(xué) 理學(xué)院，西安710055

1 引言

在人工智能與大數(shù)據(jù)時代，推薦系統(tǒng)、圖像處理和運(yùn)動結(jié)構(gòu)重建等諸多應(yīng)用領(lǐng)域的科學(xué)研究都是在處理數(shù)據(jù)矩陣。這些矩陣不僅規(guī)模龐大，還可能伴隨數(shù)據(jù)缺失、被污染和存在異常值等情況。解決上述問題的常用方法是低秩分解，例如主成分分析（principal component analysis，PCA）[1]、奇異值分解（singular value decomposition，SVD）[2]和非負(fù)矩陣分解（non-negative matrix factorization，NMF）[3]等。傳統(tǒng)的矩陣分解方法一般采用l2范數(shù)來度量逼近誤差，但它對數(shù)據(jù)中的異常值具有很高的敏感性。

為了克服PCA 對稀疏噪聲的敏感性，文獻(xiàn)[4]提出了魯棒主成分分析（robust PCA，RPCA），該模型假設(shè)噪聲是稀疏的，通過求解非凸的l1范數(shù)最優(yōu)化問題來獲得低秩逼近。隨后，一些學(xué)者研究了RPCA的凸優(yōu)化模型[5-8]，并提出了求解模型的主成分追蹤方法（principal component pursuit，PCP）[9]。在魯棒主成分分析等魯棒矩陣分解模型中，仍假設(shè)低秩矩陣的元素是確定性的，這可能會導(dǎo)致過擬合現(xiàn)象，也不利于研究數(shù)據(jù)的生成方式。為了避免上述弊端，概率低秩矩陣分解模型假設(shè)低秩成分和噪聲矩陣均服從某種隨機(jī)分布[10]。對于概率低秩矩陣分解模型，可以通過期望最大化（expectation maximization，EM）算法、Gibbs 采樣或變分貝葉斯來求解模型參數(shù)[11-20]。為了增強(qiáng)模型的魯棒性，可假設(shè)噪聲矩陣的元素服從拉普拉斯分布[21-22]。

概率矩陣分解（probabilistic matrix factorization，PMF）[23]和魯棒概率矩陣分解（robust PMF，RPMF）[24]是兩種重要的概率分解模型。在PMF 中，均假設(shè)低秩矩陣和噪聲矩陣服從高斯分布。RPMF 基于馬爾科夫鏈蒙特卡洛（Markov chain Monte Carlo，MCMC）技術(shù)，能達(dá)到比PMF 更高的預(yù)測精度。但從本質(zhì)上講，這兩種方法都等價于l2范數(shù)優(yōu)化模型，它們對數(shù)據(jù)的異常值依然比較敏感。而魯棒貝葉斯矩陣分解（robust Bayesian matrix factorization，RBMF）[25]和貝葉斯魯棒矩陣分解（Bayesian robust matrix factorization，BRMF）[26]在非確定性框架下解決模型的魯棒性問題。RBMF 將貝葉斯概率矩陣分解中的高斯噪聲替換成高斯混合噪聲，并且采用重尾的學(xué)生t分布作為低秩矩陣的先驗分布。為了減輕異常值對模型結(jié)果的不利影響，BRMF綜合運(yùn)用了高斯先驗和拉普拉斯噪聲，比其他算法具有更好的魯棒性。前述概率低秩矩陣分解模型均將數(shù)據(jù)矩陣分解為兩個低秩矩陣的乘積，使得模型缺乏靈活性和實用性。最近，文獻(xiàn)[27]提出的矩陣三分解模型（tri-decomposition model，Tri-Decom）將數(shù)據(jù)矩陣分解為低秩、稀疏噪聲和稠密噪聲三個分量矩陣之和，但該模型沒有考慮低秩矩陣的分解，也沒有使用概率模型。為此，本文建立了魯棒概率矩陣三分解模型（robust probabilistic matrix tri-factorization，RPMTF），并提出了求解該模型的條件EM算法。

2 預(yù)備知識

定義1（矩陣范數(shù)）矩陣A=(aij)m×n∈?m×n的l1范數(shù)為，F(xiàn)robenius范數(shù)為。

定義2（Kronecker 積）設(shè)A=(aij)m×n∈?m×n，B=(bst)p×q∈?p×q，稱如下的分塊矩陣

為A與B的Kronecker積。

定義3（矩陣向量化）設(shè)A=(aij)m×n∈?m×n，稱mn維列向量

vec(A)=(a11,a21,…,am1,a12,a22,…,am2,…,a1n,a2n,…,amn)T為矩陣A的逐列向量化。

定理1（奇異值分解）對于秩為r的矩陣A∈?m×n，它的奇異值分解為：

其中，U∈?m×m和V∈?n×n均為正交矩陣，對角矩陣Σr=diag(σ1,σ2,…,σr)的對角線元素滿足σ1≥σ2≥…σr＞0。

本文用到的概率分布如表1所示，其中廣義逆高斯分布中的Kp(·)表示二階修正的貝塞爾函數(shù)，c=，d=b/a。

3 魯棒概率矩陣三分解

3.1 模型建立

假設(shè)數(shù)據(jù)矩陣D=(dij)m×n∈?m×n具有近似低秩結(jié)構(gòu)，通常使用矩陣分解方法進(jìn)行維數(shù)約減和噪聲移除。在加性高斯噪聲腐蝕下，可使用奇異值分解來得到D的最優(yōu)秩k逼近矩陣L，即：

其中，k＜min{m,n}，S∈?k×k為對角矩陣，U∈?m×k與V∈?n×k滿足UTU=VTV=Ik，Ik為k階單位矩陣。

將矩陣U和V按列分塊，即U=(u·1,u·2,…,u·k)，V=(v·1,v·2,…,v·k)，易知：

為了更加靈活地獲得矩陣的低秩逼近，本文考慮將數(shù)據(jù)矩陣D的低秩成分重新分解為三個矩陣的乘積，即有如下公式：

與式（1）相比，上式中E∈?m×n是噪聲矩陣，且不要求矩陣S為對角方陣。下面建立魯棒概率矩陣三分解模型。

假設(shè)矩陣U、S、V和E均是連續(xù)型隨機(jī)變量，對應(yīng)的概率密度函數(shù)分別記為p(U)、p(S)、p(V)和p(E)。根據(jù)極大后驗估計、期望最大化或變分貝葉斯等方法可求得D的低秩近似?？紤]數(shù)據(jù)矩陣D受到大的稀疏噪聲的腐蝕，為增強(qiáng)矩陣分解模型的魯棒性，可假設(shè)噪聲矩陣E服從拉普拉斯分布。本文給出一類簡單的概率矩陣三分解模型：

矩陣集合{U,S,V}的后驗概率密度為：

對上式取對數(shù)，得：

其中，Const是與U、S、V無關(guān)的常數(shù)項。根據(jù)極大后驗估計法，最優(yōu)的低秩逼近矩陣可通過求解下列最優(yōu)化問題來獲得：

Table 1 Several commonly-used probability distributions表1 幾種常用的概率分布

最優(yōu)化問題（5）的目標(biāo)函數(shù)是非光滑的，故可使用分層模型來處理拉普拉斯分布[28]。因為：

故建立針對拉普拉斯分布的分層模型：

隨機(jī)變量矩陣T=(tij)m×n為隱變量，{D,U,S,V,T}為完整數(shù)據(jù)。結(jié)合式（3）和式（7），使用條件期望最大化方法估計D的最優(yōu)低秩逼近。

3.2 條件EM算法

記θ={U,S,V}為待估參數(shù)矩陣集合。采用條件期望最大化（conditional expectation maximization，CEM）[29]方法求解最優(yōu)的參數(shù)θ。具體而言，在算法迭代過程中先求隱變量tij的倒數(shù)的期望，再使用最大化模型更新U、S和V。令上次迭代過程得到的參數(shù)為，據(jù)此得到參數(shù)的更新公式。

（1）更新V

①E步

根據(jù)貝葉斯規(guī)則，V的后驗分布可表示為：

對上式取對數(shù)，得:

其中，Const是不依賴于V的常數(shù)。

結(jié)合式（7），得到tij的條件概率分布密度：

②M步

通過最大化Q函數(shù)來更新矩陣V。記：

（2）更新U

（3）更新S

為了最大化H(S)，將矩陣S逐列向量化，得到k2維列向量vec(S)。于是：

再將vec(S)矩陣化，得到S。

下面列出求解概率矩陣三分解的條件EM算法。

算法1求解RPMTF的條件EM算法

在上述算法中，可設(shè)置終止條件為達(dá)到最大迭代次數(shù)或：

其中，ε為一個比較小的正常數(shù)，{Uiter,Siter,Viter}表示第iter次迭代的結(jié)果。

3.3 計算復(fù)雜度

下面討論算法1 在一次迭代過程中的計算復(fù)雜度。不失一般性，假設(shè)m≤n。E 步的更新過程主要涉及殘差矩陣R的計算，因為USVT的計算復(fù)雜度為O(mnk+mk2)，所以E步計算復(fù)雜度均為O(mnk+mk2)。對于U、V和S更新過程的M步，對應(yīng)的計算復(fù)雜度分 別為O(mk2+k3)、O(nk2+k3) 和O(mnk4+k6) 。因此，一次迭代中總計算復(fù)雜度為O(mnk4+k6)。可以看出：算法1 的計算復(fù)雜度主要集中在S的更新上。在實際應(yīng)用中，k的取值不宜太大。當(dāng)k的取值存在上界時，算法1 在一次迭代過程中的計算復(fù)雜度為O(mn)，即它線性地依賴于數(shù)據(jù)矩陣的元素數(shù)目。

3.4 與其他模型的聯(lián)系

本文提出的概率矩陣三分解算法本質(zhì)上等價于根據(jù)極大后驗法求解模型（5）。若令λ→0+，則算法1可用來求解基于l1范數(shù)的奇異值分解或主成分分析。與傳統(tǒng)的概率矩陣分解模型相比，概率矩陣三分解多了矩陣S。當(dāng)S為對角矩陣、且λs→+∞時，矩陣三分解就變成了二分解；當(dāng)S為對角矩陣且sii可表示為伯努利分布與高斯分布之積時，矩陣三分解等價于穩(wěn)定魯棒主成分分析[30]。與奇異值分解相比，概率矩陣三分解模型不要求矩陣U和V滿足列正交性，而且S也不是對角的。因此，概率矩陣三分解使得模型更加靈活實用。

還可以使用張量Tucker分解[31]來表示式（2），即：

其中，“×1”表示核心張量S與模式矩陣U的1-模式積，“×2”表示核心張量S與模式矩陣V的2-模式積。將矩陣三分解向量化，得到：

文獻(xiàn)[32]提出了指數(shù)族張量分解，即假設(shè)vec(D)服從參數(shù)(V?U)vec(S)下的指數(shù)族分布，但此模型并未考慮U、S和V的先驗分布。

4 數(shù)值實驗

將魯棒概率矩陣三分解應(yīng)用到圖像去噪及視頻背景建模中。實驗采用個人計算機(jī)，其處理器為Intel?CoreTMi5-7400@3 GHz，使用Matlab R2012a 進(jìn)行編程。

4.1 圖像去噪

采用人工方式合成一幅256×256的低秩圖像（對應(yīng)矩陣L），并在其上添加若干簡單的英文單詞（對應(yīng)稀疏噪聲矩陣E），從而合成含噪圖像（對應(yīng)矩陣D）[33-34]，如圖1 所示。根據(jù)魯棒概率矩陣三分解（RPMTF）來恢復(fù)低秩成分，從而得到噪聲矩陣。在進(jìn)行實驗時，先將矩陣L、E和D的元素按相同的公式標(biāo)準(zhǔn)化到區(qū)間[-1,1]，再將英文單詞對應(yīng)的元素用[-0.5,1.5]區(qū)間上的均勻分布替換。

Fig.1 Synthesis of input image圖1 輸入圖像的合成

使用低秩恢復(fù)誤差和F1測度兩種指標(biāo)來度量算法的性能。低秩成分的恢復(fù)誤差定義為，其值越小，低秩成分的逼近性能越好。記的絕對值最小元素和最大元素分別為Emin和Emax。對于給定的實數(shù)x，定義如下兩個函數(shù)：

使用RPMTF 方法對前述合成的圖像進(jìn)行低秩與噪聲成分的分離，并與PCP、RPMF、BRMF 和Tri-Decom 的結(jié)果進(jìn)行比較。在RPMTF 參數(shù)設(shè)置中，取k=10，λu=λv=λs=λ=1，最大迭代次數(shù)為100。按照標(biāo)準(zhǔn)高斯分布隨機(jī)初始化矩陣U和V，S的初始化矩陣為U?D(VT)?，其中“?”表示矩陣的偽逆。在RPMF 和BRMF 的實驗中，仍取k=10，最大迭代次數(shù)為100，其他參數(shù)按照默認(rèn)值設(shè)置。由于RPMF、BRMF、Tri-Decom 和RPMTF 的結(jié)果具有隨機(jī)性，將實驗重復(fù)20 次，最終報告實驗結(jié)果的平均值與標(biāo)準(zhǔn)差。5種低秩矩陣恢復(fù)方法的實驗結(jié)果如表2所示。

Table 2 Comparison of low rank matrix recovery results on synthetic images表2 合成圖像的低秩矩陣恢復(fù)結(jié)果比較

從表2 可以看出：RPMTF 得到了最小的低秩恢復(fù)誤差，Tri-Decom的恢復(fù)誤差最大；BRMF取得了最大的F1 測度，RPMF 和RPMTF 得到了次優(yōu)的F1 測度，Tri-Decom 的F1 測度最?。籖PMF 的運(yùn)行時間最短，RPMTF具有最長的運(yùn)行時間。此外，BRMF低秩恢復(fù)誤差的標(biāo)準(zhǔn)差較大，穩(wěn)定性能較差。圖2繪出了5 種方法得到的低秩圖像與噪聲圖像，其中第1 行為低秩圖像，第2行為噪聲圖像。從該圖可以看出：Tri-Decom在低秩圖像和噪聲圖像方面均具有最差的恢復(fù)性能，而BRMF 和RPMTF 都較好地恢復(fù)了低秩與噪聲成分。

結(jié)合表2和圖2可以得出如下結(jié)論：PCP和RPMF的低秩恢復(fù)誤差較大，F(xiàn)1 測度較小，實驗效果較差，這可能是由較小的k值所致；Tri-Decom 的恢復(fù)誤差最大，F(xiàn)1 測度最小，實驗效果最差，這可能是因為圖像合成過程中沒考慮稠密噪聲；BRMF 和RPMTF 的恢復(fù)誤差小，F(xiàn)1測度大，實驗效果較好。

Fig.2 Low rank and noise images recovered by different methods圖2 不同方法恢復(fù)的低秩與噪聲圖像

4.2 靈敏度分析

算法1涉及4個超參數(shù){λu,λs,λv,λ}，其取值對實驗結(jié)果產(chǎn)生一定的影響。下面以前述人工圖像為例，探討RPMTF方法中超參數(shù)不同取值對實驗結(jié)果的影響。固定λ=1，取λu∈{1,7}，λv∈{1,7}，λs∈{1,7}。仍令迭代次數(shù)為100，重復(fù)實驗20次。表3列出了平均實驗結(jié)果，圖3 顯示了部分低秩圖像與噪聲圖像。從表3和圖3可以看出，RPMTF在一定程度上具有穩(wěn)定性。

Table 3 Comparison of experimental results under different combinations of hyperparameters表3 超參數(shù)不同組合下的實驗結(jié)果比較

Fig.3 Partial low rank and noise images under different hyperparameters(λu,λs,λv)圖3 不同超參數(shù)(λu,λs,λv)下的部分低秩與噪聲圖像

下面考慮更多超參數(shù)取值對實驗結(jié)果的影響。為簡便起見，固定λu=λv=1，取λs=i/2,λ=j/2，其中i∈[10],j∈[10]。實驗結(jié)果如圖4所示。從該圖可以看出：當(dāng)λs≤4,λ≤4 時，低秩恢復(fù)誤差與F1測度相對穩(wěn)定。當(dāng)λ比較大時，稀疏噪聲的方差比較大，對恢復(fù)性能產(chǎn)生不利的影響。此外，當(dāng)λs給定時，λ的取值對實驗結(jié)果影響不顯著。

4.3 視頻背景建模

在Hall視頻數(shù)據(jù)集上進(jìn)行視頻背景建模實驗[34]。選取此視頻序列的前200幀圖像，其中每幀圖像的大小均為144×176，部分圖像如圖5 所示。因此得到維數(shù)為25 344×200的數(shù)據(jù)矩陣。使用峰值信噪比（peak signal-to-noise ratio，PSNR）來度量算法的性能，其定義為：

Fig.4 Effect of different(λs,λ)on experimental results圖4 (λs,λ)的不同取值對實驗結(jié)果的影響

Fig.5 Partial images of Hall video圖5 Hall視頻的部分圖像

Fig.6 Separation of background and foreground under different methods圖6 不同方法下背景與前景的分離

其中，MSE是真實背景和重建背景的均方誤差。峰值信噪比越大，說明失真越少，實驗效果越好。

在實驗中，取k=5，最大迭代次數(shù)仍為100。對應(yīng)于某幀圖像，5種方法的背景與前景分離結(jié)果如圖6 所示，其中第1 行為背景圖像，第2 行為前景圖像。從圖6可以看出：前兩種方法PCP和RPMF背景中還存在陰影，分離效果不理想；后3 種方法BRMF、Tri-Decom 和RPMTF 的實驗結(jié)果比前兩種方法好很多。5種方法PCP、RPMF、BRMF、Tri-Decom和RPMTF的運(yùn)行時間（單位：s）分別為144.341 9、129.054 7、409.786 7、56.981 2、525.747 4，其中BRMF和RPMTF的運(yùn)行時間較長，可以考慮采用并行計算策略。5種方法的峰值信噪比如表4 所示。從表4 可以看出：Tri-Decom的峰值信噪比最小，RPMTF的峰值信噪比最大，這說明RPMTF的實驗效果最好。

Table 4 PSNR under different methods表4 不同方法的峰值信噪比

5 結(jié)束語

本文提出了概率矩陣三分解的一個魯棒模型。為了簡化模型求解，將拉普拉斯分布表示為高斯分布與指數(shù)分布的混合?；跇O大后驗策略，提出期望最大化算法。實驗結(jié)果驗證了所提模型的可行性與有效性。在建立模型時，僅考慮了U、S和V的元素服從均值為0的相互獨(dú)立的高斯分布，對于其他的概率分布仍值得進(jìn)一步研究。此外，魯棒概率張量分解和基于變分貝葉斯推斷的魯棒概率矩陣三分解也將是今后的兩個研究方向。