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

王 臻，楊 敏

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

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

王 臻，楊 敏

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

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

穩(wěn)健主成分分析；交替方向法；標(biāo)準(zhǔn)隨機(jī)k-SVD；快速隨機(jī)k-SVD

0 引 言

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

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

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

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

1.1 RPCA模型

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

(1)

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

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

(2)

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

1.2 交替方向法

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

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

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

(3)

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

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

(4)

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

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

(5)

對(duì)于L的極小化也很容易得到：

(6)

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

(7)

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

(8)

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

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

2 隨機(jī)矩陣算法

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

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

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

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

Fori=1,2,…,ndo；

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

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

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

End For

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

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

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

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

(3)理論保證

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

3 改進(jìn)的SVD算法

3.1 標(biāo)準(zhǔn)隨機(jī)k-SVD算法

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

算法2：標(biāo)準(zhǔn)隨機(jī)k-SVD算法。

輸入：目標(biāo)秩k和矩陣A∈Rm×n；

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

(9)

算法推導(dǎo)：

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

(10)

(11)

該算法有以下幾個(gè)特性：

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

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

3.2 快速隨機(jī)k-SVD算法

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

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

(12)

(13)

算法描述如下：

算法3 ：快速隨機(jī)k-SVD算法。

輸入：目標(biāo)秩k和矩陣A∈Rm×n；

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

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

算法推導(dǎo)：

(14)

基于此，通過(guò)式(15)進(jìn)行分解。

(15)

這里需要注意：

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

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

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

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

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

4 仿真實(shí)驗(yàn)

4.1 單張圖像分解仿真

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

圖1 奇異值衰減圖解

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

圖2 三種SVD分解效果圖

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

圖3為標(biāo)準(zhǔn)隨機(jī)k-SVD算法和快速隨機(jī)k-SVD算法的計(jì)算精度隨目標(biāo)秩變化的曲線。

圖3 兩種改進(jìn)算法的精度比較

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

圖4為三種SVD算法的計(jì)算時(shí)間隨目標(biāo)秩的變化曲線。

圖4 三種算法的計(jì)算時(shí)間比較

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

4.2 視頻前景提取仿真

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

圖5 視頻前景提取效果

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

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

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

表1 三種SVD算法的比較

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

5 結(jié)束語(yǔ)

針對(duì)RPCA模型中用來(lái)實(shí)現(xiàn)圖像前景和背景分離的交替方向法，對(duì)矩陣算法進(jìn)行了研究和改進(jìn)。通過(guò)分析，當(dāng)用交替方向法求解RPCA問(wèn)題時(shí)，每次計(jì)算量消耗最大的就在于核范數(shù)的優(yōu)化求解，即大規(guī)模SVD的計(jì)算。于是介紹了隨機(jī)矩陣算法，并提出了兩種新的隨機(jī)SVD算法。兩種算法都大大提高了代碼的運(yùn)行效率，并且經(jīng)過(guò)仿真實(shí)驗(yàn)表明，兩種算法達(dá)到的實(shí)際分離效果和普通的SVD分解也相差無(wú)幾?？傊?，兩種改進(jìn)的隨機(jī)奇異值分解算法都加快了算法的收斂速度，很好地提高了算法效率，并且都不需要占用大量的內(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] 江明陽(yáng),封舉富.基于魯棒主成分分析的人臉子空間重構(gòu)方法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),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] 戴瓊海,付長(zhǎng)軍,季向陽(yáng).壓縮感知研究[J].計(jì)算機(jī)學(xué)報(bào),2011,34(3):425-434.

[8] 彭義剛,索津莉,戴瓊海,等.從壓縮傳感到低秩矩陣恢復(fù):理論與應(yīng)用[J].自動(dòng)化學(xué)報(bào),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ò)出版時(shí)間

時(shí)間：2017-08-01

國(guó)家自然科學(xué)基金資助項(xiàng)目(61271234)

王 臻(1992-)，男，碩士，研究方向?yàn)楹朔稊?shù)極小化隨機(jī)優(yōu)化求解；楊 敏，博士，副教授，研究方向?yàn)橛?jì)算機(jī)視覺(jué)與圖像理解。

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