非負(fù)矩陣分解的快速收斂算法

王桂榮， 林海鵬， 靳慶貴

(1.黑龍江科技大學(xué) 機(jī)械工程學(xué)院，哈爾濱 150022;2.哈爾濱工程大學(xué) 信息與通信工程學(xué)院，哈爾濱 150001)

0 引言

信號(hào)和圖形處理中的許多問(wèn)題可以通過(guò)矩陣分解的方式來(lái)表達(dá)。非負(fù)矩陣分解(NMF)最近已被廣泛用于稀疏信號(hào)和圖形分析，甚至是普通數(shù)據(jù)處理的有效代表。NMF的主要目的是獲取數(shù)據(jù)或信號(hào)中潛在的重要結(jié)構(gòu)，進(jìn)而從這些結(jié)構(gòu)中挖掘出更多生物內(nèi)在信息。

較系統(tǒng)的 NMF理論1999年由D.D.Lee和H.S.Seung在《Nature》上提出，同時(shí)用于人臉識(shí)別[1]，引起了科學(xué)界的廣泛關(guān)注。近十幾年來(lái)，有關(guān)NMF方面的算法和應(yīng)用的理論研究也在積極的進(jìn)行中。在這些研究中，主要解決不同的代價(jià)(損失)函數(shù)和附加限制條件可以用來(lái)獲得不同類(lèi)型的矩陣分解。在額外的條件限制如非負(fù)、稀疏、光滑、低復(fù)雜性或較好可預(yù)見(jiàn)性等被強(qiáng)加到相互并不獨(dú)立的信號(hào)源上時(shí)，非負(fù)矩陣分解算法可以用來(lái)在盲信號(hào)分離中進(jìn)行有效取出或分離稀疏光滑信源。NMF實(shí)際上是對(duì)于代價(jià)函數(shù)進(jìn)行近似非線(xiàn)性?xún)?yōu)化問(wèn)題。

自L(fǎng)ee和Seung提出這兩個(gè)經(jīng)典乘性算法后，大量NMF算法被提出，其應(yīng)用過(guò)程中都存在一些限制，如一些NMF算法為了在實(shí)際應(yīng)用中增強(qiáng)某些因子特性，在代價(jià)函數(shù)中增加一些限制條件，從而形成限制NMF算法。在這些分離中主要包括稀疏性[2－3]、正交性[4]、空間位置[5]、時(shí)間連續(xù)性[6]和去相關(guān)性[7]等等限制條件。

在NMF問(wèn)題的研究文獻(xiàn)中很多使用倍乘更新MU算法對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化，顯然MU算法已經(jīng)成為NMF研究中最受歡迎的。但是對(duì)其收斂性的證明卻很少涉及。Lee和Seung[8]以歐幾里德距離EUC和散度或稱(chēng)為相對(duì)熵KL為目標(biāo)函數(shù)來(lái)證明，通過(guò)MU優(yōu)化算法目標(biāo)函數(shù)迭代減少。文獻(xiàn)[9－10]證明目標(biāo)函數(shù)的限制值既不是局部最小解，也不能連續(xù)迭代收斂到固有點(diǎn)。MU算法由于結(jié)構(gòu)簡(jiǎn)單及易于使用，成為NMF研究的基本算法，因此，對(duì)MU算法的收斂性研究顯得非常重要。

1 加速算法IILSMU-EUC

阻礙LSMU(Seung的MU算法)算法收斂主要有兩個(gè):一是這些乘性公式不能使目標(biāo)函數(shù)足夠降低。這一點(diǎn)可以從對(duì)LSMU算法形成過(guò)程的分析中得出。二是在考慮KKT時(shí)W和H可能出現(xiàn)零值。這使得算法不能保證收斂到一階固有點(diǎn)。

文獻(xiàn)[4]提出另一種簡(jiǎn)單方法解決這種問(wèn)題，即一個(gè)稍微修改的算法(MU被重寫(xiě)為變比例的梯度降低方法并且步進(jìn)長(zhǎng)度相應(yīng)修改)。主要理論如下:

對(duì)于任何常數(shù)ε≥0，V≥0和任何(W，H)≥ε，用式(1)更新‖V－WH‖:

文獻(xiàn)[4]提出的改進(jìn)算法雖然有效，但收斂速度有待提高，在上述改進(jìn)算法的基礎(chǔ)上，重新提出一個(gè)加速算法。

1.1 計(jì)算量

在MU算法中更新W所需要的計(jì)算量要首先確定(對(duì)于VT=HTWT，H因子的更新完全與W對(duì)稱(chēng))。原MU是在每次迭代時(shí)交替更新W和H。帶來(lái)的問(wèn)題是應(yīng)該更新多少次W，進(jìn)而更新多少次H都無(wú)法確定。即應(yīng)該執(zhí)行多少次MU是最好的。

為了確定計(jì)算量，先對(duì)稀疏和密集矩陣進(jìn)行有效分析，引入?yún)?shù)K表示在V中非零項(xiàng)的數(shù)量。如果V是稀疏的，那么K=mn，假設(shè)NMF已經(jīng)完成壓縮(實(shí)際中是通常要求的)，這意味著存儲(chǔ)W和H比存儲(chǔ)V要便宜，粗略地說(shuō)，在V中非零項(xiàng)的數(shù)量一定比W和H的數(shù)量大，既K≥r(m+n)。

表1給出在MU更新W矩陣相乘的浮點(diǎn)運(yùn)算數(shù)量(flops)，從中可以查到不同矩陣運(yùn)算的計(jì)算花費(fèi)，尤其矩陣相乘次序變化時(shí)計(jì)算花費(fèi)有很大不同。

函數(shù)max使成分智能化，另外，相應(yīng)算法的每一個(gè)限制點(diǎn)是以下優(yōu)化算法的固有點(diǎn):

表1 MU更新W運(yùn)算量Table 1 Flop of W in MU

表1中可以看出:在MU算法中若假設(shè) K≥r(m+n)，第一步計(jì)算(A=VHT)計(jì)算量是最大的。因此花費(fèi)時(shí)間多的運(yùn)算步驟應(yīng)盡可能少的執(zhí)行，即應(yīng)該充分利用已經(jīng)計(jì)算出的VHT和HHT矩陣積，在下一次更新H之前，完成更新W多次，即在MU中完成步驟3和4多次。

1.2 內(nèi)部迭代

從LSMU-EUC算法可以看出，依據(jù)浮點(diǎn)運(yùn)算量，估計(jì)出W的第一次更新對(duì)下一次是多么費(fèi)時(shí)(H固定時(shí))。由下列迭代因子pW給出(對(duì)于H的相應(yīng)數(shù)值由pH給出):

對(duì)于不同圖像數(shù)據(jù)庫(kù)降維的pW和pH的值，如表2所示。由于主要目的是用降維矩陣來(lái)提取原矩陣的主要數(shù)據(jù)特征，為了減少貯存或運(yùn)算時(shí)的數(shù)據(jù)量，所以降到表示維數(shù)?？勺⒁獾剑瑢?duì)于K≥r(m+n)，看到pW≥2說(shuō)明W的第一次更新比下一次花費(fèi)計(jì)算多2倍。對(duì)于一個(gè)密集矩陣K=mn，并且確切，它的確很大，因?yàn)閚是比r更大的量。例如，當(dāng)IILSMU-EUC算法對(duì)W能更新次的時(shí)間花費(fèi)相當(dāng)于LSMU-EUC算法兩次獨(dú)立更新。

表2 r=30和 r=60時(shí)的[pW]、[pH]Table 2 pWand pHwhen r=30 and 60

一個(gè)較自然的選擇是依賴(lài)pW和pH的數(shù)值更新W和H一個(gè)固定倍數(shù)。為了更柔性化可以引入?yún)?shù)a≥0，H的下次更新之前W被更新(1+apW)次;W的下次更新之前H被更新(1+apH)次。

不管怎樣，矩陣V的行的數(shù)量m和列數(shù)量n不在同一個(gè)數(shù)量級(jí)上，例如當(dāng)n＞＞m時(shí)，有，pW＞＞pH。因此，在一方面，矩陣W比H有顯著少項(xiàng)(mr＜＜nr)，并且相應(yīng)NNLS(凸?fàn)罘秦?fù)最小二乘)子問(wèn)題把很少數(shù)量變量作為特征;在另一方面，pW＞＞pH和更多W的更新被執(zhí)行。更明確地說(shuō)，很多迭代執(zhí)行在處理簡(jiǎn)單的問(wèn)題上，似乎是很不合理的。例如對(duì)于CBCL 人臉數(shù)據(jù)庫(kù)，有 m=361，n=2 429，r=20，可得到pW≈123、pH=18，用W的上百次更新才能得到NNLS子問(wèn)題近似優(yōu)化解顯得非常不必要。

1.3 內(nèi)部迭代停止標(biāo)準(zhǔn)

另一個(gè)考慮是什么時(shí)候更新W到H(更新H到W)應(yīng)該使用一個(gè)適當(dāng)?shù)臉?biāo)準(zhǔn)，因此采用下列停止標(biāo)準(zhǔn)。W(l)是W的第l次更新的迭代值(保持V固定)。只要滿(mǎn)足以下公式停止優(yōu)化迭代:

δ是一個(gè)小的正數(shù)，通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)δ=0.01時(shí)取得很好的效果。

1.4 IILSMU-EUC算法

在上述基礎(chǔ)上形成新的快速M(fèi)U算法(IILSMUEUC算法)。要求:數(shù)據(jù)矩陣V∈和初始化迭代(W，H×Rr×n

(1)當(dāng)停止標(biāo)準(zhǔn)不滿(mǎn)足時(shí)計(jì)算A=VHT，B=HHT，W(0)=W;

(2)循環(huán)(1+apW)次;

(3)使用MU計(jì)算W(l);

(4)若‖W(l+1)－W(l)‖F(xiàn)≤0.01‖W(1)－W(0)‖F(xiàn)不滿(mǎn)足到(6)步;

(5)到(7)步;

(6)結(jié)束如果判斷;

(7)結(jié)束循環(huán);

(8)W←W(l);

(9)從W更新H，V使用與 W更新相對(duì)應(yīng)(2)～(9)步;

(10)結(jié)束。

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

數(shù)值仿真實(shí)驗(yàn)比較以下算法性能:

(1)LSMU-EUC:Li和 Seung的倍乘算法[8];

(2)IILSMU-EUC:文中提出的快速M(fèi)U算法;

(3)PG:Lin[11]推薦的映射梯度算法;

(4)HALS:Cichocki和 Zdunek[12]提出的 HALS算法。

所有的測(cè)試都在MATLAB7.1，intel Core Dual CPU、內(nèi)存為1 G的 PC上運(yùn)行程序。實(shí)驗(yàn)取ε=10－16，δ=0.01，降維秩數(shù)取 r=60。為了驗(yàn)證算法的真實(shí)和有效性，分別針對(duì)密集性CBCL人臉數(shù)據(jù)庫(kù)和稀疏性O(shè)RL人臉數(shù)據(jù)庫(kù)圖像分別進(jìn)行實(shí)驗(yàn)，結(jié)果如下。

2.1 密集性CBCL數(shù)據(jù)庫(kù)仿真

CBCL人臉數(shù)據(jù)庫(kù)是361×2 429密集型數(shù)據(jù)庫(kù)。圖1顯示目標(biāo)函數(shù)‖V－WH‖對(duì)于CBCL密集數(shù)據(jù)在 a為 0、0.5、1、2和 4時(shí)，在 10 s內(nèi)收斂情況。

圖1 誤差曲線(xiàn)Fig.1 Error curve

觀(guān)察到原LSMU-EUC算法(a=0)除了在初始階段比 IILSMU-EUC算法誤差小外，經(jīng)過(guò)1 s后IILSMU-EUC算法目標(biāo)函數(shù)明顯小于原LSMU-EUC算法，說(shuō)明文中第1部分中分析Lee和Seung的原算法不能夠使目標(biāo)函數(shù)足夠降低是正確的。IILSMU-EUC算法收斂性顯著好于原LSMU-EUC算法，也證明IILSMU-EUC算法的實(shí)用價(jià)值。另外，對(duì)于IILSMU-EUC 算法實(shí)驗(yàn)分別取 a=0.5，a=1，a=2，a=4四種情況，可發(fā)現(xiàn)只是在最初的2 s內(nèi)對(duì)誤差影響較大，但是經(jīng)過(guò)2 s以后IILSMU-EUC算法產(chǎn)生的誤差都趨近相同，綜合整體性能a=1時(shí)加速算法從收斂特性和誤差性能都是最好的。

圖2是對(duì)于CBCL數(shù)據(jù)庫(kù)分別采用Lee和Seung的LSMU-EUC算法、IILSMU-EUC算法、Lin提出的PG算法、Cichocki和Zdunek提出的HALS算法進(jìn)行實(shí)驗(yàn)，得出的誤差曲線(xiàn)。

圖2 誤差曲線(xiàn)Fig.2 Error curve

從圖2所示的曲線(xiàn)可以觀(guān)察到，收斂性順序如下:LSMU-EUC算法最差，其次PG算法，稍差I(lǐng)ILSMU-EUC算法，最好的是HALS;經(jīng)過(guò)10 s后誤差也呈現(xiàn)上述排序。但是從算法的復(fù)雜性上考慮，PG和HALS算法都遠(yuǎn)遠(yuǎn)比MU算法復(fù)雜，因此IILSMU-EUC算法比原LSMU-EUC算法性能上有很大提高，又具有算法簡(jiǎn)單、收斂速度快及易于使用特點(diǎn)。

圖3～5分別是CBCL數(shù)據(jù)庫(kù)采樣、LSMU-EUC算法處理和 IILSMU-EUC算法(a=1)處理后的圖像。

圖3 CBCL數(shù)據(jù)庫(kù)人臉采樣Fig.3 Face samples of CBCL

圖4 MU算法處理結(jié)果Fig.4 MU consequence

圖5 IILSMU-EUC算法處理結(jié)果Fig.5 IILSMU-EUC consequence

由圖5可見(jiàn)，運(yùn)用IILSMU-EUC算法分解的圖像噪聲小，主要特征提取明顯，其性能上比原算法分解處理的圖像(圖4)要好。這是由于10 s處理后IILSMU-EUC算法的誤差小于原LSMU-EUC算法處理的誤差，因此直觀(guān)的得出文中算法從精密性上好于原算法。

2.2 稀疏性O(shè)RL人臉圖像仿真

ORL人臉數(shù)據(jù)庫(kù)是一個(gè)2 576×400稀疏數(shù)據(jù)矩陣，從圖6可以看出，與處理密集數(shù)據(jù)CBCL時(shí)有同樣的收斂特性，IILSMU-EUC算法處理的誤差曲線(xiàn)在2 s時(shí)刻已經(jīng)變得很小，已經(jīng)低于LSMU-EUC算法10 s后的誤差值，說(shuō)明IILSMU-EUC算法收斂性明顯好于LSMU-EUC算法，IILSMU-EUC算法在a=1時(shí)整體性能最佳。

圖6 誤差曲線(xiàn)Fig.6 Error curve

圖7是LSMU-EUC算法、IILSMU-EUC算法、Lin提出的PG算法和HALS算法處理稀疏矩陣ORL時(shí)誤差曲線(xiàn)比較，從圖7中可發(fā)現(xiàn)小于4 s時(shí)PG誤差值大于LSMU-EUC，以后誤差值小于LSMU-EUC算法，收斂特性曲線(xiàn)發(fā)現(xiàn)HALS最優(yōu)，其次文中算法，4 s后PG優(yōu)于LSMU-EUC算法。

圖7 誤差曲線(xiàn)Fig.7 Error curve

值得一提的是比較圖2與圖7，可以看到前者誤差大，后者誤差小。這是由于對(duì)于CBCL降維幅度大，特征提取時(shí)誤差大。而對(duì)于ORL來(lái)講降維幅度小，所以誤差小，但這并不影響分析結(jié)果。

通過(guò)實(shí)驗(yàn)證明，無(wú)論是對(duì)密集數(shù)據(jù)還是稀疏數(shù)據(jù)進(jìn)行NMF處理時(shí)，IILSMU-EUC算法(a≠0)的誤差曲線(xiàn)上明顯好于原LSMU-EUC算法(a=0)，并且在a=1時(shí)錯(cuò)誤率最少。從圖2和圖7可得:PG比LSMU-EUC算法優(yōu)越，一般情況下IILSMU-EUC算法都優(yōu)于LSMU-EUC算法;IILSMU-EUC算法性能上優(yōu)于PG算法，略次于HALS算法。

3 結(jié)論

在對(duì)Lee和Seung的LSMU-EUC算法進(jìn)行性能分析中找出該算法的缺陷，針對(duì)該改進(jìn)算法的矩陣運(yùn)算的計(jì)算量進(jìn)行分析，優(yōu)化計(jì)算順序和考慮迭代方法，提出IILSMU-EUC算法。IILSMU-EUC算法中花費(fèi)時(shí)間多的運(yùn)算步驟應(yīng)盡可能少的執(zhí)行，也就是，充分利用已經(jīng)計(jì)算出的MVT和VVT矩陣積，在下一次更新V之前，完成更新U多次;同樣對(duì)V更新也同樣操作，結(jié)果加快迭代收斂速度，提高了分解性能。仿真結(jié)果顯示在a=1時(shí)，IILSMU-EUC算法目標(biāo)函數(shù)明顯小于原算法，即新算法收斂速度快;誤差曲線(xiàn)分析結(jié)果顯示，IILSMU-EUC算法優(yōu)于PG算法，略次于HALS算法，但新算法簡(jiǎn)單且易于使用;圖像分解的對(duì)比實(shí)驗(yàn)結(jié)果說(shuō)明，新算法對(duì)圖象主要特征提取明顯，噪聲小及精密性好。所以IILSMUEUC算法具有算法簡(jiǎn)單、易于使用和收斂速度快，因此更具有推廣價(jià)值。

[1]LEE D D，SEUNG H S.Learning the parts of objects by nonnegative matrix factorization[J].Nature，1999，401(4):788 －791.

[2]HOYER P O.Non negative matrix factorization with sparseness constraints[J].Mach Learning，2004，5(9):1457 －1469.

[3]CHOI S.Algorithms for orthogonal nonnegative matrix factorization[C]//In Proc of IEEE International Joint Conference on Neural，USA:Networks，2008:1828 －1832.

[4]LI S，HOU X，ZHANG H，CHENG Q.Learning spatially localized，parts－based representation[C]//In Proc of IEEE Computer Society Conference on Computer Vision and Pattern Recognition，USA:Hawaii，2001:207 －212.

[5]VIRTANEN T.Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria[J].IEEE Transactions on Audio，Speech and Language Processing，2007，15(3):1066－1074.

[6]ZHANG Y，F(xiàn)ANG Y.A NMF algorithm for blind separation of uncorrelated signals[C]//In Proc of IEEE International Conference on Wavelet Analysis and Pattern Recognition，China:Beijing，2007:999－1003.

[7]GONZALES E F，ZHANG Y.Accelerating the Lee-Seung algorithm for nonnegative matrix factorization[J].Comput Math Rice Univ Houston，2005，26(3):47－51.

[8]LEE D D，SEUNG H S.Algorithms for non-negative matrix factorization[J].In Advances in Neural Information Processing，2001，13(5):456－460.

[9]LIN C J.On the convergence of multiplicative update algorithms for nonnegative Matrix factorization[J].In IEEE Trans on Neural Networks，2007，8(6):117 －122.

[10]BADEAU R，BERTIN N，VINCENT E.Stability analysis of multiplica-tive update algorithms and application to non-negative matrix factorization[J].IEEE Trans Neural Network，2010，21(12):1869－1881.

[11]LIN C J.Projected gradient methods for nonnegative matrix factorization[J].Neural Computation，2007，19(6):2756 －2779.

[12]CICHOCKI A，ZDUNEK R，AMARI S.Hierarchical ALS algorithms for nonnegative matrix and 3D tensor factorization[J].In Lecture Notes in Computer Science，2007，466(5):169－176.

[13]李 樂(lè)，章毓晉.非負(fù)矩陣分解算法綜述[J].電子學(xué)報(bào)，2008，36(4):737－743.

[14]張永鵬，鄭文超，張曉輝.非負(fù)矩陣分解及其在圖像壓縮中的應(yīng)用[J].西安郵電學(xué)院學(xué)報(bào)，2008，13(3):58－61.

[15]劉開(kāi)第，金 斕，龐彥軍，等.與判斷矩陣一致性無(wú)關(guān)的單準(zhǔn)則排序法[J].河北工程大學(xué)學(xué)報(bào):自然科學(xué)版，2013，30(1):107－112.