基于非負矩陣分解的低秩矩陣恢復模型

徐夢珂許道云魏明俊

（1.貴州大學數學與統計學院貴陽550025）（2.貴州大學計算機科學與計算學院貴陽550025）

基于非負矩陣分解的低秩矩陣恢復模型

徐夢珂1許道云2魏明俊2

（1.貴州大學數學與統計學院貴陽550025）（2.貴州大學計算機科學與計算學院貴陽550025）

針對低秩矩陣恢復需要求解大規模矩陣核范數奇異值分解，計算復雜度高的缺陷，提出基于非負矩陣分解的低秩矩陣恢復模型。新模型通過對傳統低秩矩陣恢復模型中的低秩矩陣進行非負因子分解，不但可以保持原始數據的局部特征，而且其低秩性可以快速求解矩陣低秩分解，從而避免了矩陣核范數求解大規模奇異值分解問題。在算法上采用多乘子交替迭代法（ADMM），將全局問題分解為多個易求解的局部子問題，對每個子問題利用拉格朗日乘子法分別對低秩矩陣和稀疏矩陣進行迭代求解。在ORL，AL_Gore和Windows三個圖像數據庫中Matlab仿真實驗結果表明，新模型求解算法比傳統低秩矩陣恢復模型識別率高，降秩效果明顯，算法的時間復雜度低，從而提高算法運行速度。

非負矩陣分解；低秩矩陣恢復；多乘子交替迭代法；奇異值分解；圖像識別

Class NumberTP391.41

1 引言

在大數據時代，如何高效率地處理大規模數據越來越受到人們重視。在現實生活中，巨大的數據量往往給數據分析帶來麻煩，即所謂的“維度災難”［1］。雖然許多數據都可以用矩陣形式表達，但是直接處理矩陣沒有意義，矩陣分解就是將矩陣“分而治之”為數個矩陣的運算。非負矩陣分解是1999年由D.D.Lee和H.S.Seung首次提出，其基本思想是尋找一組基矩陣和線性組合參數來近似原始數據，同時約束表達基和參數都是非負的［2］。非負矩陣分解形式描述如下：

給定非負矩陣A∈Rm×n，k是所降的維度，求解最優化問題：其中W為基矩陣，W各列表示原始數據的潛在結構，具有稀疏性。H為系數矩陣，H各列表示各個結構在該原始數據中所占的比重?！珹F是Frobenius范數：

一方面，矩陣進行非負矩陣分解后，其低秩和稀疏性使得算法速度快；另一方面，現實中的很多物理特性都是非負的，負數沒有實際意義，非負約束使得模型具有很強的可解釋性。

對非負矩陣分解的算法有很多，大體上分為兩類，一類是公式逼近法，優點是迭代公式構造簡單，主要的運算是矩陣之間加法、乘法和求逆等，但收斂性難以得到證明［3］。另一類是優化逼近法，如梯度投影法、二次規劃序列逼近法等，優點是算法都比較成熟，但是只能保證解的局部最優性，因此仍然需要對算法進行改進［4～5］。

更具挑戰的是一些大規模數據可能含有空缺或者毀損元素。主成分分析法是利用K-L變換將原始數據投影到一個低維特征空間，提取主要特征［6］，該方法可以用在人臉識別中。然而實際中，人臉圖像是一個二維矩陣，主成分分析把矩陣變成向量，破壞了矩陣的空間結構，其性能對于有噪聲的情況缺乏魯棒性。2011年John Wright等提出來的低秩矩陣恢復解決了上述問題［7］。低秩矩陣恢復把向量的稀疏表示推廣到矩陣的低秩情形，將含有噪音的數據矩陣表示為低秩矩陣與稀疏噪聲矩陣之和，再通過求解矩陣核范數優化問題從含“噪音”原始矩陣中恢復出“干凈”的低秩矩陣。同時，求解低秩矩陣近似的一系列算法研究也在積極展開。

2009年，Beck，Teboulle用迭代閾值算法來求解矩陣核范數最小化問題［8］，2010年Cai、Candes、and Shen提出了奇異值閾值算法，在求矩陣奇異值時用軟閾值算子代替，該算法僅對于秩比較低的矩陣有效［9］。Lin、Chen and Wu提出用增廣拉格朗日乘子法求解低秩矩陣恢復，增廣拉格朗日乘子法比前面的算法速度快［10］。但是以上算法每一步迭代都要求解核范數奇異值分解問題，如果矩陣大小是n×n時，它的奇異值計算時間復雜度是O(n3)，所以當矩陣規模比較大時，這些算法計算復雜度會比較高，因此如何有效處理矩陣核范數奇異值分解是個難點。

非負矩陣分解通過對矩陣進行非負因子分解，能夠快速求解矩陣的低秩分解，從而可以避免矩陣核范數求解大規模奇異值分解問題。因此，受到非負矩陣分解的啟發，本文研究如何提高低秩矩陣恢復的魯棒性，采用低秩矩陣的非負矩陣分解技術，提出基于非負矩陣分解的低秩矩陣恢復模型，新模型不但可以保持原始數據的局部特征，而且避免求解矩陣核范數大規模奇異值分解問題，從而降低時間復雜度。在算法方面，由于基于非負矩陣分解的低秩矩陣恢復模型涉及到多個變量，所以采用多乘子交替迭代算法，將變量求解分成多個子問題，固定其他子問題，交替最小化求解每一個子問題的增廣拉格朗日函數。基于非負矩陣分解的低秩矩陣恢復模型求解與傳統的低秩矩陣恢復模型在ORL，AL_Gore和Windows三個圖像數據庫的應用結果相比，提高了圖像的識別率和低秩矩陣的降秩能力，降低了時間復雜度，加快了算法的運行速度，進而表明了低秩矩陣恢復的非負矩陣分解技術的有效性。

2 基于非負矩陣分解的低秩矩陣恢復模型

矩陣低秩近似把壓縮感知理論通過一維向量拓展到二維矩陣上。矩陣的稀疏性分為兩類：矩陣元素的稀疏性和矩陣奇異值的稀疏性。元素的稀疏性指矩陣非零元素個數比較少；矩陣奇異值的稀疏性指矩陣奇異值非零元素個數比較少。若同時考慮矩陣元素和奇異值的稀疏性，可以得到低秩矩陣恢復模型。低秩矩陣恢復模型描述如下［2］

其中X∈Rm×n是數據矩陣，Z∈Rm×n是低秩矩陣，Z的秩rank(Z)=r，r＜＜min(m，n)，E∈Rm×n是稀疏噪音矩陣。低秩矩陣恢復解決了這樣一個問題：假設原始數據含有噪音，就將含有噪音的數據矩陣表示為低秩矩陣與稀疏噪聲矩陣之和，再通過求解矩陣核范數優化問題從含有“缺損”數據矩陣中恢復“干凈”的低秩矩陣。

然而矩陣秩的估計問題是NP難的，關鍵問題是，如何處理矩陣的秩？Candes和Recht從理論上解決了這個問題：如果一個矩陣的秩足夠小，類比壓縮感知中的l0范數和l1范數，自然而然把矩陣的秩用矩陣核范數來代替，從而轉化為核范數的最小化問題［11］。設矩陣A∈Rm×n，則矩陣核范數‖‖A*：，即矩陣奇異值之和。

核范數是凸的，最小化核范數就是一個凸優化問題。2011年Emmanuel和Candes在文獻［7］給出了矩陣秩問題如何轉化為矩陣核范數凸優化問題的證明，即下面定理1和2。

定理1［7］假設矩陣A0∈Rm×n（m≥n），如果存在參數μ＞0，有：

則稱A0以參數μ滿足非相干性條件。

其中ei是單位向量，矩陣A0的奇異值分解為是矩陣A0的秩，范數

定理2［7］在非相干性條件下，從一個秩為r的矩陣A0中均勻抽取m個元素，滿足：，則核范數最小化問題的最優解Aˉ依至少1-cn-3的概率接近原來矩陣A0，其中C，c是常數。

但是很多求解低秩矩陣恢復模型的算法每一步迭代都要求解核范數奇異值分解問題，當矩陣大小是m×n時，它的奇異值分解時間復雜度是O(n3)，所以考慮到當矩陣規模比較大時，以上這些算法時間復雜度會比較高。因此，基于非負矩陣分解是通過對矩陣進行非負因子分解，能夠快速求解矩陣低秩分解的想法，本文采用低秩矩陣的非負分解技術，提出基于非負矩陣分解的低秩矩陣恢復模型，其數學表達式如下：

X∈Rm×n是原始數據矩陣，Z∈Rm×n是低秩矩陣，E∈Rm×n是稀疏噪音矩陣，是Z的非負因子。通過把非負矩陣分解運用到低秩矩陣恢復中，不但可以保持原始數據的局部特征，而且避免求解大規模矩陣核范數的奇異值分解問題，從而對原始數據矩陣進行降秩，降低算法的時間復雜度，提高算法運行速度。

3 基于非負矩陣分解的低秩矩陣恢復模型求解

由于基于非負矩陣分解的低秩矩陣恢復模型涉及到多個變量，所以本文用多乘子交替迭代算法來求解目標函數［12］，其增廣拉格朗日函數表示為［10］：

其中X∈Rm×n是原始數據矩陣，Z∈Rm×n是低秩矩陣，E∈Rm×n是稀疏噪音矩陣，是Z的非負因子，為兩個拉格朗日乘子，μ＞0為懲罰因子，A，B表示矩陣的內積，增廣拉格朗日函數消除了等式約束，在函數里同時增加一次懲罰項和二次懲罰項，采用迭代方法可求得函數解和拉格朗日乘子。

交替迭代求解各個子問題：

1）首先固定其他變量，求解子問題Z，Z的增廣拉格朗日函數表示如下

對Z的求解會用到核范數的次微分，通過文獻［13］，Z的最優解Z*的求解可以通過奇異值收縮得到如下的閉式解：

那么可以得到：

其中(σ-δ)+=max(0，σ-δ)。

所以Z的最優解Z*的求解式（7）可以轉化成：2）固定其他變量，求解子問題E，求解E有下面的優化模型：

為了方便計算以及簡化符號，設A=X-WH-E+Y1/μ+Y2/μ，δ=1/μ。

即Z*=SVTδ(A)。

由于矩陣A的SVD分解為

根據文獻［14］，該模型有如下閉式解：

其中[，]:，i矩陣的第i列且bi指矩陣的第i列。

3）固定其他變量，求解子問題W和H：

對W和H求解是最小二乘優化問題［15］，讓D=Z+Y2/μ，當固定變量H時，W=DHT，根據文獻［15］，本文設計如下迭代：

其中QR(DHT)表示通過QR分解得到矩陣DHT的一個正交基的操作。

當固定變量W，H=WTD，同理，對H的迭代如下：

4）對拉格朗日乘子的更新：

綜上所述，基于非負矩陣分解的低秩矩陣恢復模型求解算法步驟如下：

輸入原始數據：X∈Rm×n。

初始化：Z=0，E=0，W=0，H=eye(m，n)，Y1=0，Y2=0，μ。

輸出：Z，E

步驟：（1）根據式（8）更新Z；

（2）根據式（10）更新E；

（3）根據式（13）更新W；

（4）根據式（14）更新H；

（5）根據式（15）（16）更新拉格朗日乘子Y1，Y2；

（6）更新懲罰因子μ；

（7）迭代停止條件：

重復1～7，直到停止迭代，返回Z和E。

4 數值實驗及分析

4.1 算法復雜度分析

首先對基于非負矩陣分解的低秩矩陣恢復模型求解算法進行時間復雜度分析：因為該模型求解的主要算法包括：低秩矩陣Z∈Rm×n的奇異值分解，求解非負因子W、H的QR分解和一些矩陣乘法計算這三個主要部分。其中Z的秩rank(Z)=r，r＜＜min(m，n)。

低秩矩陣Z奇異值分解時間復雜度是：O(r3)；QR分解時間復雜度為O(mr2+nr2)；其余一些矩陣乘法所需時間復雜度為O(rmn)；所以基于非負矩陣分解的低秩矩陣恢復模型算法總的時間復雜度為O(r3+r2(m+n)++rmn)。

低秩矩陣恢復算法主要是求解原始矩陣的奇異值分解，時間復雜度是O(m3)或者是O(n3)，都遠遠大于基于非負矩陣分解的低秩矩陣恢復模型算法的時間復雜度O(r3+r2(m+n)++rmn)，因此基于非負矩陣分解的低秩矩陣恢復模型算法求解的時間復雜度比傳統低秩矩陣恢復模型要低。

4.2 數值實驗及分析

本文從ORL，AL_Gore和Windows三個圖像數據庫中分別隨機選取10張圖像作為樣本進行Matlab數值實驗，分別把傳統模型的交替迭代方法（ADM）和基于非負矩陣分解的低秩矩陣恢復模型的多乘子交替迭代法（ADMM）來作為實驗對比。交替迭代方法是解決凸優化問題的增廣拉格朗日函數法的進一步改進。在偏微分方程中，對于結構型變分不等式，增廣拉格朗日求解沒有充分利用結構型變分不等式可分離特性，為此Gabay和Mercier［16］首先提出了交替迭代方向法，其思想是先求解一個子變分不等式問題得到一個解x，再利用新的迭代點x去求解另一個子變分不等式得到迭代點y，最后用x、y去迭代拉格朗日乘子。

在基于非負矩陣分解的低秩矩陣恢復模型實驗中，本文設正則參數λ=0.03，懲罰因子原始數據X的所有元素

絕對值的平均值。采用算法恢復矩陣與原始矩陣的誤差作為識別率，其定義如下

其中，X∈Rm×n是原始數據矩陣，Z∈Rm×n是低秩矩陣，E∈Rm×n是稀疏噪音矩陣。

首先在400張ORL人臉數據庫中隨機選取10幅56×56原始數據圖像X，分別從原始數據矩陣X得到低秩矩陣Z的降秩效果，算法運行時間和所恢復矩陣的識別率這三個方面，對傳統低秩矩陣恢復模型的交替迭代法（ADM）和基于非負矩陣分解的低秩矩陣恢復模型的多乘子交替迭代法（ADMM）的實驗結果進行對比，結果如表1所示。

表1 ORL數據庫ADM和ADMM實驗結果比較

為了更好地驗證本文所提算法的有效性，再次在140張AL_Gore人臉數據庫中隨機選取10幅107×107原始數據圖像X進行實驗，結果如表2所示。

為了驗證新模型的ADMM算法對大數據的降秩效果，本文在18張AL_Gore數據庫中隨機選取10幅大小為1200×1600的圖像，為了方便實驗，把圖像規模預處理為1200×1200的原始圖像X進行實驗，結果如表3所示。

從表1，表2和表3比較結果可以得出：本文提出的基于非負矩陣分解的低秩矩陣恢復模型比傳統的低秩矩陣恢復模型求解算法不僅僅識別率高，時間復雜度低，運行時間短，而且可以明顯地降低原始矩陣的秩。

表2 AL-Gore數據庫ADM和ADMM實驗結果比較

表3 Windows數據庫ADM和ADMM實驗結果比較

為了比較出算法對不同秩的原始數據X矩陣的降秩效果，任意選取了6幅秩不等的圖像矩陣X，把基于非負矩陣分解的低秩矩陣恢復模型的ADMM算法，傳統低秩矩陣恢復模型的ADM算法和利用梯度下降來求低秩和稀疏矩陣方法的近端梯度算法（APG）進行比較，結果如圖1，從圖中可以得出：隨著原始數據矩陣秩的增加，基于非負矩陣分解的低秩矩陣恢復模型的ADMM算法降秩效果最明顯，而ADM算法降秩效果不佳。

圖1 ADM，ADMM，APG算法對不同維數的原始數據降維效果

最后，本文把基于非負矩陣分解的低秩矩陣恢復模型的ADMM算法和傳統低秩矩陣恢復模型的ADM算法的識別率進行比較，結果如圖2。

圖2 ADM，ADMM算法的識別率比較

從圖2可以得出：隨著原始數據矩陣的秩的增加，基于非負矩陣分解的低秩矩陣恢復模型的ADMM算法識別率高，而ADM算法識別率逐漸降低。

5 結語

低秩矩陣恢復是重要的數據分析工具，低秩特性解決了大規模數據效率問題，稀疏捕捉了大數據中的關鍵信息，兩者抓住了主要矛盾問題，是處理大數據的核心技術，在人臉識別、推薦系統等領域有了不少的應用。本論文觀察到矩陣分解的低秩性可以降低算法的時間復雜度，因此把非負矩陣分解的思想加入到低秩矩陣恢復模型中，采用把低秩矩陣部分進行非負矩陣分解技術，不但可以得到原始數據的局部特征，避免了求解大規模矩陣核范數的奇異值分解，從而達到降秩效果在算法上采用了多乘子交替迭代算法框架，降低了時間復雜度，從而提高算法的運行速度，而且矩陣分解結果不含負值元素，具有實際意義。本論文進一步改進是：1）進一步挖掘大規模數據的潛在結構如正定性，結構稀疏性等，從而提高矩陣恢復的識別率。2）算法迭代求解非負因子W，H過程中，考慮如何更加有效地處理負元素，進一步去提高算法精度。

［1］David L，Donoho.High-dimensional data analysis：the cures and blessings of dimensionality［C］//Proc.Amer. Math Soc.Conf.Math Challenges Lecture，2000.

［2］Daniel D，Lee H，Sebastians.Learning the parts of objects by non-negative matrix factorization［J］.Nature，1999，

（401）：788-791.

［3］JIANG J.Improvement of non-negative matrix decomposition methods and its applications［D］.Beijing：Beijing University of Techology，2011：1-60.

［4］張可村，李換琴.工程優化方法及其應用［M］.西安：西安交通大學出版社，2007：83-156.

ZHANG Kecun，LI Huanqin.Engineering optimization method and its application［M］.Xi'an：Xi'an Jiaotong University Press，2007：83-156.

［5］陽明盛，羅長童.最優原理、方法及求解軟件［M］.北京：科學出版社，2006：46-132.

YANG Mingsheng，LUO Changtong.Best principles，methods and software［M］.Beijing：Science Press，2006：46-132.

［6］Hotelling H.Analysis of a complex of statistical variables into principal components［J］.Journal of Educational Psychology，1933，24：417-441.

［7］Candes E，Li X，Ma Y，Wright J.Robust principal component analysis？［J］.Journal of ACM，2011，58（3）：1-37.

［8］Beck A，Teboulle M.A fast iterative shrinkage thresholding algorithm for linear inverse problems［J］.SIAM Journal on Imaging Sciences，2009，2（1）：183-202.

［9］Cai J，Candes E，Shen Z.A singular value thresholding algorithm for matrix completion［J］.SIAM Journal on Optimization 2010，20（4）：1956-1982.

［10］Lin Z，Chen M，Wu L.The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices［C］//In Technical Report UILU ENG，2009.

［11］Candes E，Recht B.Exact matrix competion via convex optimization［J］.Foundation of Computational Math，2009，9（6）：717-772.

［12］XIAO M，JUN F.Sparse and low-rank matrix decomposition via alternation direction methods［J］.Pacific Journal of Optimization，2013，9（1）：167-180.

［13］Cai J，Emmanuel J，Candes.A singular value thresholding algorithm for matrix completion［J］.SIAM，2010，1956-1982.

［14］Liu Y，Jiao L，Sheng F.A fast tri-factorization method for low-rank matrix recovery and completion［J］.Pattern Recognition，2013，46（1）：163-173.

［15］Yuan S，Zai W，Yin Z.Augmenteed lagrangian alternating direction method for matrix separation based on low-rank factorization［J］.Optimization method and Software，2014，29（2）：239-263.

［16］Gabay D，Mercier B.A dual algorithm for the solution of nonlinear variational problems via finite elements approximations［J］.Compute Math Appl，1976（2）：178-40.

Low-rank Matrix Recovery Model Based on Non-negative Matrix Factorization

XU Mengke1XU Daoyun2WEI Mingjun2
（1.College of Mathematics and Statistics，Guizhou University，Guiyang550025）（2.College of Computer Science and Technology，Guizhou University，Guiyang550025）

To overcome the shortage of large-scale nuclear matrix singular value decomposition existing in low-rank matrix recovery model，the paper proposed low-rank matrix recovery model based on non-negative matrix factorization.Non-negative matrix factorization（NMF）applied to the low-rank matrix，which could quickly deal with the problem of the decomposition matrix of low-rank and avoid large-scale nuclear matrix singular value decomposition.Then the algorithm used alternarting directions method of multipliers（ADMM）.ADMM divided the global problem into partial sub-problems.Each sub-problem used Lagrange multipliers to solve low rank matrix and sparse matrix.Experimental results in ORL，AL_Gore and Windows databases showed that low-rank recovery model based NMF has higher recognition rate，better reduction rank and lower the complexity of the algorithm than other traditional low-rank recovery model.

non-negative matrix factorization（NMF），low-rank matrix recovery，alternating directions method of multipliers，singular value decomposition（SVD），image recognition

TP391.41

10.3969/j.issn.1672-9722.2017.06.002

2016年12月18日，

2017年1月25日

國家自然科學基金項目（編號：61262006，61540050）；貴州省重大應用基礎研究項目（編號：黔科合JZ字［2014］2001）；貴州省科技廳聯合基金（編號：黔科合LH字［2014］7636）資助。

徐夢珂，女，碩士研究生，研究方向：密碼學理論與工程。許道云，男，博士，教授，研究方向：可計算性與計算復雜性、算法設計與分析。魏明俊，男，碩士研究生，研究方向：可計算性與計算復雜性。