王 臻，楊 敏

(南京郵電大學(xué) 自動化學(xué)院，江蘇 南京 210023)

核范數(shù)隨機矩陣求解新方法及其RPCA應(yīng)用

王 臻，楊 敏

(南京郵電大學(xué) 自動化學(xué)院，江蘇 南京 210023)

RPCA(穩(wěn)健主成分分析)從原始觀測數(shù)據(jù)中恢復(fù)低秩成分和稀疏成分。RPCA常用交替方向法迭代求解，算法的效率取決于核范數(shù)優(yōu)化求解，即SVD分解。而RPCA在計算機視覺應(yīng)用中，圖像和視頻等巨大的數(shù)據(jù)量為大規(guī)模數(shù)據(jù)SVD分解帶來了很大困難。采用隨機矩陣算法對SVD分解進行改進，分別為計數(shù)縮略算法、標準隨機k-SVD算法和快速隨機k-SVD算法。主要是對原有大規(guī)模數(shù)據(jù)矩陣進行降維隨機采樣，使用隨機投影算法得到原數(shù)據(jù)矩陣的一個近似，對這個近似矩陣進行QR分解，得到對應(yīng)的酉矩陣。對酉矩陣進行相關(guān)操作，得到與原矩陣SVD相似的結(jié)果。算法的時間效率和存儲空間得到極大改善?；趩螐垐D像和視頻前景檢測等仿真實驗，表明所提方法大大提高了RPCA迭代優(yōu)化求解的效率。

穩(wěn)健主成分分析；交替方向法；標準隨機k-SVD；快速隨機k-SVD

0 引 言

PCA在高維數(shù)據(jù)樣本中尋找和挖掘低維特征空間。在較小的高斯隨機噪聲時，可通過奇異值分解準確求解。當不滿足上述條件時，所得結(jié)果會有很大偏差，于是用RPCA來解決此種情況。RPCA即低秩矩陣恢復(fù)[1]，又稱為稀疏與低秩矩陣分解[2]。

RPCA模型(穩(wěn)健主成分分析)[3]在視頻前景提取、人臉識別圖像預(yù)處理等計算機視覺領(lǐng)域中應(yīng)用廣泛。求解RPCA常用交替方向法[4]，充分利用RPCA模型具有的很好的可分離結(jié)構(gòu)特性，是對帶有線性約束條件凸規(guī)劃問題的增廣拉格朗日乘子法的一種改進。算法每次迭代的主要計算量在于SVD分解[5]。像矩陣求逆、特征值分解、奇異值分解等這些矩陣計算，不但非常耗時，存儲所需的空間也很大，因此，這些缺點限制了它的拓展和應(yīng)用范圍。為了對一個大規(guī)模的矩陣進行計算，隨機矩陣近似技術(shù)[6]應(yīng)運而生。隨機矩陣計算已經(jīng)在現(xiàn)代數(shù)據(jù)科學(xué)領(lǐng)域占有舉足輕重的地位。特別在過去的十年里，隨機數(shù)字線性代數(shù)取得了卓越進步，現(xiàn)在大規(guī)模的矩陣計算不再是一個不可完成的任務(wù)了。

文中將改進的SVD分解，即隨機矩陣算法應(yīng)用到視頻圖像的前景檢測中。對于交替方向法中每次迭代時，對核范數(shù)的優(yōu)化求解，即大規(guī)模SVD分解，進行優(yōu)化改進，對于大規(guī)模的數(shù)據(jù)矩陣進行隨機近似求解，大大提高了算法效率。數(shù)據(jù)表明，改進算法具有更快的計算速度和存儲優(yōu)勢。

1 RPCA模型及交替方向法求解算法

1.1 RPCA模型

RPCA問題可以抽象描述[7]為：已知觀測數(shù)據(jù)矩陣M，且M=L+S，而L和S是未知的，但是已知L為低秩，S為稀疏且非零元素可以任意大，要求恢復(fù)L。基于問題的描述，首先想到的解決方法是尋求觀測數(shù)據(jù)中主成分L的最小秩，且干擾誤差矩陣S是稀疏的，即非零元素個數(shù)較少。于是形成了如下優(yōu)化問題：

(1)

其中，‖·‖0表示S中非零元的個數(shù)。

通過對式(1)做松弛變化，可以把問題轉(zhuǎn)化為一個易于解決的凸優(yōu)化問題，即用1范數(shù)代替0范數(shù)，用L的核范數(shù)代替L的秩。式(1)轉(zhuǎn)化為：

(2)

這個凸優(yōu)化問題，在一定條件下，可以求解并能得到理想的恢復(fù)結(jié)果(L,S)，且解是唯一的。

1.2 交替方向法

交替方向法面向目標函數(shù)可分解[8](即為兩個不同變量凸函數(shù)的和)、帶線性約束和凸集約束的問題，它是Lagrange乘子法的一種變形。先構(gòu)造問題的增廣Lagrangian函數(shù)，然后通過交替極小化增廣Lagrangian函數(shù)來更新兩個變量，最后更新Lagrange乘子，如此迭代。其好處是更新變量子問題要比原問題簡單，甚至有閉解。

對于穩(wěn)健PCA問題(式(2))，易得其增廣拉格朗日函數(shù)為：

Γ(L,S,Y)=‖L‖*+λ‖S‖1+〈Y,L+S-

(3)

其中，Y∈Rm×n為線性約束乘子；μ>0為標量，是違背線性約束的懲罰參數(shù)；〈·〉表示標準內(nèi)積；‖·‖F(xiàn)表示Frobenius范數(shù)。

很明顯，可以直接運用經(jīng)典增廣拉格朗日乘子法求解，其迭代主題是：

(4)

每個迭代是一個大型的非光滑優(yōu)化問題。經(jīng)典的增廣拉格朗日乘子法把式(2)看作一個一般的最小化問題，而忽視了它在目標函數(shù)和約束條件中都體現(xiàn)出的很好的可分離結(jié)構(gòu)特性，即式(4)是同時最小化L和S的。

而對于S的極小化是非常容易得到的：

(5)

對于L的極小化也很容易得到：

(6)

交替更新，可以得到結(jié)果：

(7)

其中，prox是在如下空間的歐幾里得投影：

(8)

懲罰參數(shù)比較適合于動態(tài)調(diào)整，從相關(guān)文獻[9-11]可以了解具有動態(tài)變化參數(shù)的交替方向法的收斂性以及一些具體有效的調(diào)整參數(shù)的方法。

交替方向法求解RPCA問題(式(4))時，算法效率取決于每一次迭代的核范數(shù)優(yōu)化求解，即大規(guī)模的SVD分解。因此，SVD算法的好壞對于交替方向法處理大規(guī)模矩陣的稀疏和低秩矩陣恢復(fù)問題有很重要的影響。

2 隨機矩陣算法

在對矩陣進行操作前，通常會對大規(guī)模矩陣進行降維操作。矩陣的降維方式一般有隨機投影和列選擇兩種方式[12]。假設(shè)給定一個矩陣A∈Rm×n，而S∈Rn×s是一個縮略矩陣，比如是一個隨機投影或者是列選擇矩陣[13]，C=AS∈Rm×s是A的一個縮略。矩陣C的大小是遠小于矩陣A的，但是C卻保存了矩陣A的一些重要性質(zhì)。文中使用的降維方式主要是使用隨機投影中的計數(shù)縮略算法[14]。

算法1：計數(shù)縮略算法。

輸入矩陣A∈Rm×n及縮略矩陣的列數(shù)s；

初始化C為一個m×s的全零矩陣；

Fori=1,2,…,ndo；

(1)隨機均勻采樣矩陣C的l個列向量；

(2)隨機均勻采樣g個+1或者-1；

(3)通過將矩陣C的第l個列的各元素加上g×矩陣A的第i個列向量來更新矩陣C的第l個列；

End For

得到縮略矩陣C∈Rm×s。

該算法有以下一些性質(zhì)：

(1)時間復(fù)雜度為O(nnz(A))

(2)空間復(fù)雜度為O(ms)。當矩陣A不適合存儲時，該算法將矩陣C保存在內(nèi)存中，并且只經(jīng)歷一次矩陣A的各列向量。

(3)理論保證

(4)計數(shù)矩陣縮略算法在矩陣A是稀疏矩陣時特別有效。

3 改進的SVD算法

3.1 標準隨機k-SVD算法

本節(jié)描述標準隨機k-SVD算法[15]，它計算矩陣A的k-SVD分解并且到達了1+ε的F范數(shù)相對誤差。算法描述如下：

算法2：標準隨機k-SVD算法。

輸入：目標秩k和矩陣A∈Rm×n；

S∈Rn×s，C=AS∈Rm×s

(9)

算法推導(dǎo)：

由于C的列空間和QC的列空間是一樣的，式(9)的極小化問題就可以被等價轉(zhuǎn)換為：

(10)

(11)

該算法有以下幾個特性：

(1)該算法只有2次經(jīng)過A；

(3)時間復(fù)雜度為O(nnz(A)k/ε)。

3.2 快速隨機k-SVD算法

標準算法將大多數(shù)時間花費在解決式(10)上。假如式(10)可以更高效地解決，那么標準隨機k-SVD可以變得更加迅速?？梢宰⒁獾?，可以通過不精確的最小二乘回歸算法[17]解決式(10)的問題。

現(xiàn)在構(gòu)造一個m×p大小的計數(shù)縮略(或者是高斯投影+計數(shù)縮略)矩陣P來解決這個問題：

(12)

(13)

算法描述如下：

算法3 ：快速隨機k-SVD算法。

輸入：目標秩k和矩陣A∈Rm×n；

構(gòu)造一個n×s計數(shù)縮略矩陣S，然后令C=AS；

構(gòu)造一個m×pcs的計數(shù)縮略矩陣Pcs，還有一個pcs×p的矩陣Psrht；

算法推導(dǎo)：

(14)

基于此，通過式(15)進行分解。

(15)

這里需要注意：

(1)該算法只經(jīng)過2次矩陣A；

(2)算法的時間復(fù)雜度為O(nnz(A)+(m+n)poly(k/ε))；

(3)參數(shù)應(yīng)該滿足k

(4)算法中的縮略算法也可以使用其他的矩陣縮略方式；

(5)矩陣A，矩陣sketch，矩陣L這些都是整個系統(tǒng)中存儲花銷最大的，但是幸運的是，它們只需要1次或2次掃描。假如矩陣A，縮略矩陣，矩陣L不適合隨機存儲，那么它們應(yīng)該被存儲到硬盤中，再逐個部分被加載到內(nèi)存中進行計算。

4 仿真實驗

4.1 單張圖像分解仿真

測試圖像為MATLAB自帶圖片，其大小為512*512。如圖1所示，前10個奇異值中，奇異值的衰減速度相當快，并且逐漸趨向于0，圖像表現(xiàn)為低秩矩陣?，F(xiàn)分別使用SVD算法、標準隨機k-SVD算法以及快速隨機k-SVD算法進行低秩矩陣重建仿真，比較算法之間的運行效率、恢復(fù)精度，以及與目標秩之間的關(guān)系。

圖1 奇異值衰減圖解

重建結(jié)果與原圖的比較如圖2所示。

圖2 三種SVD分解效果圖

這是分別使用三種SVD算法，利用主成分累加和逼近原圖，得到原圖的低秩逼近效果圖?？梢钥闯?，三種算法都很好地實現(xiàn)了對于原矩陣的處理。

圖3為標準隨機k-SVD算法和快速隨機k-SVD算法的計算精度隨目標秩變化的曲線。

圖3 兩種改進算法的精度比較

由圖3可以看出，標準隨機k-SVD算法的分解精度要比快速隨機k-SVD算法高。隨著目標秩的升高，兩種隨機算法的精度都逐漸提高并趨向平穩(wěn)。

圖4為三種SVD算法的計算時間隨目標秩的變化曲線。

圖4 三種算法的計算時間比較

由圖4可以看出，兩種改進算法的計算時間都隨目標秩的增加而增加。由于快速隨機k-SVD每次計算都需要構(gòu)造好幾個縮略矩陣，所以在目標秩一定的情況下，標準隨機k-SVD比快速隨機k-SVD用時短。而MATLAB自帶的SVD每次計算的時間幾乎是不變的，比兩種改進算法都要耗時。

4.2 視頻前景提取仿真

待操作數(shù)據(jù)矩陣大小為20 800*200，取視頻中的4幀。取目標秩k為50，并且同時構(gòu)造出原矩陣的縮略矩陣C，其中s取60。仿真結(jié)果如圖5所示。

圖5 視頻前景提取效果

對于快速隨機k-SVD算法，同樣先構(gòu)造縮略矩陣C，p和pcs分別取100和80，取值條件需滿足k

如圖5所示，第一行是視頻中的4幀數(shù)據(jù)。第二行是文中改進算法對每幀數(shù)據(jù)分解得出的低秩部分，即背景。最后一行是分解得出的稀疏部分，即視頻的運動前景。可以看出，改進算法都較好地實現(xiàn)了視頻前景與背景的分離。

三種SVD算法的比較如表1所示。

表1 三種SVD算法的比較

從表1的運行結(jié)果可以看出，兩種改進的隨機k-SVD分解與普通的SVD分解一樣也實現(xiàn)了前景和背景的分離。標準隨機k-SVD算法雖然迭代次數(shù)比普通的SVD算法要多一些，但是運行時間卻大幅縮短，并且由于采用縮略矩陣算法，程序?qū)τ趦?nèi)存的占用也相對較少，很好地實現(xiàn)了算法的改進?？焖匐S機k-SVD算法雖然每次迭代的時間比標準隨機k-SVD算法要長一些，因為快速隨機k-SVD算法中進行了多次的矩陣縮略降維，但是其迭代次數(shù)比普通的SVD算法和標準隨機k-SVD算法都要少，運行時間也更短。雖然該算法中使用了多個縮略矩陣，比較占內(nèi)存，但幸運的是，它們只需1次或2次的掃描就可以實現(xiàn)算法的功能。實驗結(jié)果證明，該算法也有良好的收斂性，提高了核范數(shù)求解的效率。

5 結(jié)束語

針對RPCA模型中用來實現(xiàn)圖像前景和背景分離的交替方向法，對矩陣算法進行了研究和改進。通過分析，當用交替方向法求解RPCA問題時，每次計算量消耗最大的就在于核范數(shù)的優(yōu)化求解，即大規(guī)模SVD的計算。于是介紹了隨機矩陣算法，并提出了兩種新的隨機SVD算法。兩種算法都大大提高了代碼的運行效率，并且經(jīng)過仿真實驗表明，兩種算法達到的實際分離效果和普通的SVD分解也相差無幾?？傊?，兩種改進的隨機奇異值分解算法都加快了算法的收斂速度，很好地提高了算法效率，并且都不需要占用大量的內(nèi)存空間。

[1] Candès E J,Li X,Ma Y,et al.Robust principal component analysis[J].Journal of the ACM,2011,58(3):11-12.

[2] Lin Z,Chen M,Wu L,et al.The augmented Lagrange multiplier method for exact recovery of a corrupted low-rank matrices[R].Illinois:University of Illinois at Urbana-Champaign,2010.

[3] 江明陽,封舉富.基于魯棒主成分分析的人臉子空間重構(gòu)方法[J].計算機輔助設(shè)計與圖形學(xué)學(xué)報,2012,24(6):761-765.

[4] Yuan X,Yang J.Sparse and low-rank matrix decomposition via alternating direction methods[R].Hong Kong:Hong Kong Baptist University,2009.

[5] Wang Ruoxi, Li Yingzhou, Mahoney M W,et al.Structured block basis factorization for scalable kernel matrix evaluation[J].Numerische Mathematik,2015,83:313-323.

[6] Mahoney M W.Randomized algorithms for matrices and data[J].Foundations and Trends in Machine Learning,2013,3(2):647-672.

[7] 戴瓊海,付長軍,季向陽.壓縮感知研究[J].計算機學(xué)報,2011,34(3):425-434.

[8] 彭義剛,索津莉,戴瓊海,等.從壓縮傳感到低秩矩陣恢復(fù):理論與應(yīng)用[J].自動化學(xué)報,2013,39(7):981-994.

[9] He B S,Liao L Z,Han D,et al.A new inexact alternating directions method for monontone variational inequalities[J].Mathematical Programming,2002,92(1):103-118.

[10] He B S,Yang H.Some convergence properties of a method of multipliers for linearly constrained monotone variational inequalities[J].Operations Research Letters,1998,23(1):151-161.

[11] Kontogiorgis S,Meyer R R.A variable-penalty alternating directions method for convex optimization[J].Mathematical Programming,1998,83(1):29-53.

[12] Clarkson K,Woodruff D P.Low rank approximation and regression in input sparsity time[C]//Proceedings of annual ACM symposium on theory of computing.[s.l.]:ACM,2013:121-128.

[13] Woodruff D P.Sketching as a tool for numerical linear algebra[J].Foundations and Trends in Theoretical Computer Science,2014,10(1-2):1-157.

[14] Wang Shusen,Zhang Zhihua,Zhang Tong.Towards more efficient symmetric matrix sketching and cur matrix decomposition[J].Journal of Machine Learning Research,2015,13:2729-2769.

[15] Halko N,Martinsson P G,Tropp J A.Finding structure with randomness:probabilistic algorithms for constructing approximate matrix decompositions[J].SIAM Review,2011,53(2):217-288.

[16] Gorinov S A,Tyrtyshnikov E E.The maximum-volume concept in approximation by low-rank matrices[J].Contemp. Math.,2001,280:47-51.

[17] Halko N.Randomized methods for computing low-rank approximations of matrices[D].Colorado:University of Colorado,2012.

ANewMethodforSolvingNuclearNormwithRandomMatrixandItsApplicationinRobustPrincipalComponentAnalysis

WANG Zhen,YANG Min

(College of Automation，Nanjing University of Posts and Telecommunications，Nanjing 210023，China)

RPCA (Robust Principal Component Analysis) recovers sparse and low rank components from the original observation data.It commonly uses ADM (Alternate Direction Method) for iterative solving,the efficiency of which depends on the nuclear norm optimization solution,that is SVD.The application of RPCA in computer vision,large amounts of data from images and video make it difficult for large-scale data SVD.Therefore,a random matrix algorithm is adopted to improve the SVD,respectively the algorithm of count sketch,the prototype randomizedk-SVD and the faster randomizedk-SVD.Its main idea is to reduce the size of the original large-scale data matrix and sample randomly.Using the random projection algorithm to obtain an approximation of the original matrix,and operating QR decomposition of this approximate matrix,the unitary matrix corresponding to it is obtained,and then the results which is similar to the SVD can be achieved through correlated operation of unitary matrix.The time and space of the algorithm have been greatly optimized.Simulation based on single image and video foreground detection shows that the proposed method can greatly improve the efficiency of RPCA iterative optimization.

RPCA;ADM;prototype randomizedk-SVD;faster randomizedk-SVD

TP391

A

1673-629X(2017)12-0071-06

10.3969/j.issn.1673-629X.2017.12.016

2016-10-21

2017-02-23 < class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2017-08-01

國家自然科學(xué)基金資助項目(61271234)

王 臻(1992-)，男，碩士，研究方向為核范數(shù)極小化隨機優(yōu)化求解；楊 敏，博士，副教授，研究方向為計算機視覺與圖像理解。

http://kns.cnki.net/kcms/detail/61.1450.TP.20170801.1550.022.html