劉偉強(qiáng),林 鷺

(廈門(mén)大學(xué)數(shù)學(xué)科學(xué)學(xué)院,福建 廈門(mén) 361005)

非負(fù)矩陣分解(nonnegative matrix factorization,NMF)[1-2]是近年來(lái)提出的一種對(duì)大規(guī)模非負(fù)數(shù)據(jù)降維的方法,通過(guò)非負(fù)的限制,讓分解得到的特征矩陣以及系數(shù)矩陣具有明顯的可解釋性,使得一些潛在的信息被有效地挖掘出來(lái).NMF已經(jīng)廣泛地應(yīng)用在人臉識(shí)別[3]、圖像分析[4]、文本聚類[5]、盲源信號(hào)分類[6]、生物基因探測(cè)[7],以及生物數(shù)據(jù)挖掘等領(lǐng)域[8].

近年來(lái),學(xué)者從不同角度出發(fā),提出相應(yīng)的NMF算法.Lee等[2]于2001年首次提出交替方向的乘性迭代算法(MU).Hoyer[9-10]從罰函數(shù)的角度出發(fā),提出了帶稀疏約束條件的NMF算法.Lin[11-12]則提出了基于梯度投影的NMF算法.但對(duì)NMF問(wèn)題算法分析較多,理論分析較少.本文中從線性互補(bǔ)問(wèn)題出發(fā),得到基于不動(dòng)點(diǎn)方程的NMF(fixed point equation NMF,FP-NMF)算法,并從理論上證明算法的收斂性.根據(jù)下降方向的不同,分別提出了不動(dòng)點(diǎn)方程的最速下降(FP-steepest descent,FP-SD)算法與不動(dòng)點(diǎn)方程的最小梯度(FP-minimal gradient,FP-MG)算法,證明了這兩種算法的收斂性.

1 NMF問(wèn)題與不動(dòng)點(diǎn)方程的關(guān)系

s.t.W≥0,H≥0.

(1)

若將目標(biāo)函數(shù)看成是W與H兩個(gè)變量的函數(shù),則式(1)變?yōu)橐粋€(gè)非凸優(yōu)化問(wèn)題,但目標(biāo)函數(shù)對(duì)于單變量W或者H而言均為凸函數(shù).因此,常用交替迭代法求解.

固定W,此時(shí)的目標(biāo)函數(shù)可改寫(xiě)為

其中Vj代表矩陣V的第j列,下同.式(1)等價(jià)于求解n個(gè)非負(fù)最小二乘問(wèn)題[13],即

(2)

同理,對(duì)式(1)固定H,可轉(zhuǎn)化為求解:

(3)

(4)

由定理1可知,利用交替方向迭代求解NMF問(wèn)題時(shí),其KKT條件正好為一個(gè)線性互補(bǔ)問(wèn)題;又因問(wèn)題(2)關(guān)于H是凸函數(shù),故求解問(wèn)題(2)等價(jià)于求解一個(gè)線性互補(bǔ)問(wèn)題(4).

定理2[14]問(wèn)題(4)與下列不動(dòng)點(diǎn)方程同解:

Hj=(Hj-ηjgrad(Hj))+,ηj>0,

(5)

其中H+=max(H,0),下同.

類似地,對(duì)于問(wèn)題(3),有以下結(jié)論:

(6)

定理4[14]式(6)與下列不動(dòng)點(diǎn)方程同解:

(7)

2 FP-SD與FP-MG算法

用梯度下降法求解問(wèn)題(2)與問(wèn)題(3).固定W,按列計(jì)算新的H:

Hj∶=Hj-ηj(grad(H))j,

其中ηj是負(fù)梯度方向上的步長(zhǎng),j=1,2,…,n,grad(H)=WT(WH-V).

同理,固定H,按行計(jì)算新的W:

其中ξi是負(fù)梯度方向上的步長(zhǎng),i=1,2,…,m,grad(W)=(WH-V)HT.

分別記問(wèn)題(2)與問(wèn)題(3)的Hesse矩陣為G(H)與G(W),則

G(H)=WTW,G(W)=HHT.

現(xiàn)在考慮如何得到最優(yōu)的步長(zhǎng)ηj與ξi,使得算法具有更快的收斂速度.為此,用SD[15]優(yōu)化以下式子:

j=1,2,…,n.

(8)

同理,固定H,定義:

(9)

算法1FP-SD

對(duì)k=0,1,…，

1) 計(jì)算(Wk)TWk的譜半徑ρ((Wk)TWk);

2) 計(jì)算grad(Hk)=(Wk)T(WkHk-V),

G(Hk)=(Wk)TWk;

3) 對(duì)j=1,2,…,n,

4) 計(jì)算Hk(Hk)T的譜半徑ρ(Hk(Hk)T);

5) 計(jì)算grad(Wk)=(WkHk-V)(Hk)T,

G(Wk)=Hk(Hk)T;

6) 對(duì)i=1,2,…,m,

注意到式(8)與(9)在選取步長(zhǎng)時(shí)所用的方法為SD方法.當(dāng)然,也可以考慮使得梯度值最小所對(duì)應(yīng)的步長(zhǎng)的MG方法[16],為此,考慮以下優(yōu)化式子:

(10)

式(10)是關(guān)于步長(zhǎng)ηj的一個(gè)非負(fù)最小二乘問(wèn)題,解之可得:

同理,固定H,定義:

i=1,2,…,m.

(11)

式(11)是關(guān)于步長(zhǎng)ξi的一個(gè)非負(fù)最小二乘問(wèn)題,解之可得:

類似于Lin的方法[18-19]將式(10)與(11)應(yīng)用于問(wèn)題(1)中,得到如下基于投影梯度方法的MG算法.

算法2FP-MG

對(duì)k=0,1,…，

1) 計(jì)算(Wk)TWk的譜半徑ρ((Wk)TWk);

2) 計(jì)算grad(Hk)=(Wk)T(WkHk-V),

G(Hk)=(Wk)TWk;

3) 對(duì)j=1,2,…,n,

5) 計(jì)算Hk(Hk)T的譜半徑ρ(Hk(Hk)T);

6) 計(jì)算grad(Wk)=(WkHk-V)(Hk)T,

G(Wk)=Hk(Hk)T;

7) 對(duì)i=1,2,…,m,

3 FP-SD與FP-MG算法的收斂性

證明由已知可設(shè)WTW的特征值為0<λ1≤λ2≤…≤λr=ρ(WTW),則

又有

故有

再由引理1知:

表1 MU、FP-SD、FP-MG算法在不同秩上關(guān)于迭代次數(shù)、相對(duì)誤差和運(yùn)行時(shí)間的比較

Tab.1 Comparison of MU,FP-SD and FP-MG on the number of iterations,reletive error,and running time in different ranks

r迭代次數(shù)相對(duì)誤差運(yùn)行時(shí)間/s MUFP-SDFP-MGMUFP-SDFP-MGMUFP-SDFP-MG 53 0326979350.029 30.028 00.028 2671203302 104 2481 2462 6670.024 20.022 10.022 41 031375907 206 9891 1741 5790.020 50.018 00.018 71 905423613 308 5611 0651 4210.019 00.016 10.017 22 781385597 409181 5220.015 10.016 55751 023 605 2208411 4310.016 20.014 00.016 02 9925701 096 804 2767861 3550.015 20.013 40.015 72 643406840

4 數(shù)值實(shí)驗(yàn)

4.1 實(shí)驗(yàn)數(shù)據(jù)

選取劍橋大學(xué)的ORL人臉圖庫(kù)[20]做數(shù)值實(shí)驗(yàn).ORL圖庫(kù)包含40個(gè)人臉的圖片,每個(gè)人臉從不同角度拍攝10張相片,共400張圖片,每張圖片的像素為112×92,其元素值介于0到255之間.因此V為10 304×400矩陣.

4.2 KKT剩余量

W≥0,

grad(W)≥0,

W.*grad(W)=0,

和

H≥0,

grad(H)≥0,

H.*grad(H)=0.

上述條件可以等價(jià)簡(jiǎn)化為

min(W,grad(W))=0,

min(H,grad(H))=0，

(12)

其中式(12)左邊為按元素取最小所得的矩陣.令

則式(12)成立的充分必要條件是δ=0.顯然δ與W及H的階數(shù)有關(guān).為排除階數(shù)的影響,實(shí)際計(jì)算中,取如下的標(biāo)準(zhǔn)KKT剩余量

作為度量,其中

δW=#(min(W,grad(W))≠0),

δH=#(min(H,grad(H))≠0).

4.3 實(shí)驗(yàn)結(jié)果與分析

表1是以KKT剩余量Δ≤0.01或者迭代步長(zhǎng)上限10 000作為停機(jī)準(zhǔn)則.對(duì)每一組r,所有的算法均取相同的初始矩陣對(duì),迭代次數(shù)、相對(duì)誤差及運(yùn)行時(shí)間均為50組初始矩陣對(duì)計(jì)算后的平均值.

表中空白部分代表迭代次數(shù)超過(guò)上限10 000.觀察表1可知:

1) 從迭代次數(shù)上看,FP-SD以及FP-MG算法均優(yōu)于MU算法,且FP-SD算法表現(xiàn)最佳.當(dāng)r取20,30,40,60時(shí),MU算法未能在5 000次迭代內(nèi)達(dá)到停機(jī)準(zhǔn)則,FP-SD和FP-MG均在1 000步左右達(dá)到停機(jī)準(zhǔn)則,并且FP-SD算法幾乎全在1 000步內(nèi)達(dá)到.當(dāng)r為80時(shí),FP-SD以及FP-MG算法的迭代次數(shù)均不足MU算法的1/3,其中FP-SD不足MU的1/5.除MU算法外,我們發(fā)現(xiàn)一個(gè)有趣的現(xiàn)象,即隨著r的增長(zhǎng),迭代步數(shù)呈先增后減的倒U型趨勢(shì).

2) 從相對(duì)誤差上看,隨著r的增加,3種算法均呈現(xiàn)下降的趨勢(shì).且FP-SD以及FP-MG算法均較MU算法有較好的表現(xiàn),其中相對(duì)誤差平均下降了7%以上,有的下降了15%.

3) 從運(yùn)行時(shí)間上看,隨著r的增加,3種算法均呈現(xiàn)先增長(zhǎng)后下降的趨勢(shì).且FP-SD以及FP-MG算法較MU算法都有明顯的提高,運(yùn)行效率大多提高了2倍以上,最多的提高了7倍.

下面固定r,對(duì)照MU,FP-SD及FP-MG算法的平均相對(duì)誤差以及平均KKT剩余量.由于ORL圖庫(kù)包含40個(gè)人的人臉圖片,所以取r=40,對(duì)比結(jié)果如圖1.

圖1 三種NMF算法的平均相對(duì)誤差對(duì)比Fig.1 Comparison of the average relative error of the three NMF algorithms

圖1對(duì)比了3種算法的平均相對(duì)誤差.從圖中可以看到,FP-SD與FP-MG算法的相對(duì)誤差下降較快,且當(dāng)?shù)螖?shù)在20到40之間時(shí),這2種算法較MU算法的優(yōu)勢(shì)更為顯著.MU算法迭代近100次后,相對(duì)誤差達(dá)到0.02,而FP-SD與FP-MG算法只需50次迭代,相對(duì)誤差即可達(dá)到0.02,步長(zhǎng)減少了一半.圖1從側(cè)面驗(yàn)證了所有算法的相對(duì)誤差是單調(diào)下降的,這與前面的理論分析相吻合.平均KKT剩余量對(duì)比結(jié)果如圖2.

圖2 三種NMF算法的平均KKT剩余量對(duì)照Fig.2 Comparison of the average KKT residual of the three NMF algorithms

由圖2可以看出,迭代15次后,FP-SD以及FP-MG算法的KKT剩余量均小于MU算法.圖中還看到,要使KKT剩余量達(dá)到0.2,FP-SD以及FP-MG算法大約需要迭代40次,而MU迭代100次后仍未達(dá)到.進(jìn)一步,FP-SD算法迭代70次時(shí),KKT剩余量已低于0.1,而FP-MG需要迭代100次.總之,從KKT剩余量看,這3種算法中,FP-SD以及FP-MG算法效果較好.

5 結(jié) 論

本研究將運(yùn)籌學(xué)中的線性互補(bǔ)問(wèn)題與NMF問(wèn)題聯(lián)系起來(lái),利用NMF問(wèn)題與線性互補(bǔ)問(wèn)題的關(guān)系,將NMF問(wèn)題轉(zhuǎn)化為求解不動(dòng)點(diǎn)方程問(wèn)題,從而設(shè)計(jì)出了基于不動(dòng)點(diǎn)方程N(yùn)MF算法,討論了算法的收斂性,并從理論上給予證明.借助這個(gè)思路,利用最速下降法以及最小梯度法選取下降步長(zhǎng),分別得到基于不動(dòng)點(diǎn)方程FP-SD與FP-MG算法.數(shù)值實(shí)驗(yàn)的結(jié)果表明,本文中提出的兩種算法較MU算法在收斂速度上有明顯的改進(jìn).