魯棒概率矩陣三分解*

史加榮，陳姣姣

1.省部共建西部綠色建筑國家重點實驗室/西安建筑科技大學，西安710055

2.西安建筑科技大學 理學院，西安710055

1 引言

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

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

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

2 預備知識

定義1（矩陣范數）矩陣A=(aij)m×n∈?m×n的l1范數為，Frobenius范數為。

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

為A與B的Kronecker積。

定義3（矩陣向量化）設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(·)表示二階修正的貝塞爾函數，c=，d=b/a。

3 魯棒概率矩陣三分解

3.1 模型建立

假設數據矩陣D=(dij)m×n∈?m×n具有近似低秩結構，通常使用矩陣分解方法進行維數約減和噪聲移除。在加性高斯噪聲腐蝕下，可使用奇異值分解來得到D的最優秩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)，易知：

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

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

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

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

對上式取對數，得：

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

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

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

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

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

3.2 條件EM算法

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

（1）更新V

①E步

根據貝葉斯規則，V的后驗分布可表示為：

對上式取對數，得:

其中，Const是不依賴于V的常數。

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

②M步

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

（2）更新U

（3）更新S

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

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

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

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

在上述算法中，可設置終止條件為達到最大迭代次數或：

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

3.3 計算復雜度

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

3.4 與其他模型的聯系

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

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

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

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

4 數值實驗

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

4.1 圖像去噪

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

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

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

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

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

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

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

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

4.2 靈敏度分析

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

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

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

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

4.3 視頻背景建模

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

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

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

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

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

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

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

5 結束語

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